Kriging using ArcPy for a bunch of ponit shapefiles

540
6
05-28-2022 07:51 AM
Ehsanfrmd
New Contributor

Hi,

I am trying to perform spatial interpolation for a few set of points. I have 17 shapefiles and I need to do kriging for each of them. I tried to do it using the following code. The code is working and it saves the raster files but it looks like it has a problem because when I try to open the output using Arcmap, it gives me an error. I also have used numpy to read the raster output and it was just one number. Would you please help me with that?

 

import arcpy
from arcpy import env
from arcpy.sa import *
import os

# Specify your environment
arcpy.env.workspace = r"C:\Users\...."


filess = 'C:/Users/Desktop/python code/points'

point_file_list = [ ]
for file in os.listdir(filess):
    if ".shp" in file:
        point_file_list.append(file)

path_list = [ ]
for points in point_file_list:
    inFeatures = os.path.join(filess,points)
    field = 'dro'
    cellSize = 12500
    outRaster = os.path.join("C:/Users/Desktop/python code/rasters")
    lagSize = 2000
    majorRange = 2.6
    partialSill = 542
    nugget = 0
    kModelOrdinary = KrigingModelOrdinary("CIRCULAR", lagSize, majorRange, partialSill, nugget)
    kRadius = RadiusFixed(20000, 1)
    outKriging = Kriging(inFeatures, field, kModelOrdinary, cellSize, kRadius)
    path_list.append(outRaster)
    tifname = points[-10:]
    kriging_OUT= os.path.join(outRaster, 'r{0}.tif'.format(tifname))
    outKriging.save(kriging_OUT)

0 Kudos
6 Replies
DanPatterson
MVP Esteemed Contributor

The error?


... sort of retired...
0 Kudos
Ehsanfrmd
New Contributor

Here is the error.

0 Kudos
DanPatterson
MVP Esteemed Contributor

You should be using data with a projected coordinate system and all the data you are trying to krige should have the same coordinate system.  The error suggests that you have a mix of data coordinates.  Also. set the analysis extent and cell size within the script for each shapefile that you are using, or you will have no control over the output.


... sort of retired...
0 Kudos
Ehsanfrmd
New Contributor

But I am using the same latitudes and longitudes for all of my data. Also, I am importing the  latitudes and longitudes from a csv file. They are just numbers and there is no coordinates.

0 Kudos
DanPatterson
MVP Esteemed Contributor

lat/long in decimal degrees are the coordinates.  you should be using projected data


... sort of retired...
0 Kudos
Ehsanfrmd
New Contributor

Thank you very much

0 Kudos