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)
The error?
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.
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.
lat/long in decimal degrees are the coordinates. you should be using projected data
Thank you very much