Perform Kriging interpolatation using ArcPy

ShouvikJha · ‎11-22-2016

I have multiple point feature class (r001_mean, r002_mean.....012_mean) into one folder. I am trying to write a python script to perform Kriging interpolation which will loop all the point feature class and outraster will save as it name of input raster.

Each Point feature class hold one Z filed only (GRID_CODE)

Below is my script, what i have so far.....its giving error massage of "not defined point_num"

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

In_Point = r"D:\NPP_Test\MERRA_TEMP_2012C" #(Point feature name:r001_mean, r002_mean.....r012_mean )
Out_Raster = r"D:\NPP_Test\MERRA_TEMP_2012D"

points = arcpy.ListFeatureClasses()


zFields = "GRID_CODE"

#Kriging Veriable
cellSize = 0.05
lagSize = 0.5780481172534
majorRange = 6
partialSill = 3.304292110
nugget = 0.002701348
kRadius = RadiusFixed(20000, 1)

#Mask region of interest
mask = r"D:\Gujarta Shape file\GUJARATSTATE.shp"

for zField in zFields:
    Point = Point_Num[:3]
    kModelUniversalObj = KrigingModelUniversal("LINEARDRIFT", lagSize, majorRange, partialSill, nugget)

    OutKriging = Kriging(inPointFeatures, zField, kModelUniversalObj, cellSize, kRadius)


    #IDWMASk = ExtractByMask(outIDW, mask)
    KrigMask = ExtractByMask(OutKriging, mask)
    #Save outraster as the same name of input
    KrigMask.save("r{}.tif".format(Point_Num))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

XanderBakker · ‎11-24-2016

There are a number of things going wrong in your code.

you have to include a "import os" at the start of the script, since you use it later on
change zFields = ["GRID_CODE"] by zField = "GRID_CODE", since IDW will require a string indicating the field, not a list
the "default" search distance in the example script is based on the (projected) data used for the example, this should be adjusted to a valid value for your data. For instance I used 5 degrees.
the cellsize is very small and will generate a very detailed result (which is hard to defend based on the input data). I used 0.025 in my example to run your data

Below the results based on the points in r001_mean.shp:

Code used:

import arcpy
from arcpy import env
from arcpy.sa import *
import os
arcpy.env.parallelProcessingFactor = "100%"
arcpy.env.overwriteOutput = True

# settings
env.workspace = r"D:\Xander\GeoNet\KrigingMultiFC\KRIGING\Test_Points"
mask = "D:\Xander\GeoNet\KrigingMultiFC\KRIGING\MASK_FILE\GUJARAT-STATE_NEW.shp"
out_folder = r"D:\Xander\GeoNet\KrigingMultiFC\KRIGING\OUTPUT_RAS"
zField = "GRID_CODE" # has to be string, not list

points = arcpy.ListFeatureClasses("r*_mean.shp", "POINT")
print points

cellSize = 0.025 # 0.002298707671
power = 2
searchRadius = RadiusVariable(10, 5) # smaller max distance

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Execute IDW
for fc in points:
    print fc
    outIDW = Idw(fc, zField, cellSize, power, searchRadius)
    # Execute Mask
    IDWMASk = ExtractByMask(outIDW, mask)

    #Save output raster to output workspace
    tifname = fc[:4]  # captures "r001" for output tiff name
    print tifname

    # Save output, except Solar Radiation raster
    IDW_OUT= os.path.join(out_folder, 'r{0}.tif'.format(tifname))
    IDWMASk.save(IDW_OUT)

#print done
print 'done'‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

View solution in original post

JoshuaBixby · ‎11-22-2016

My guess, line 29 is generating the error: Point = Point_Num[:3] . You are trying to slice Point_Num but it hasn't been defined anywhere in the script. Is it being inherited?

ShouvikJha · ‎11-22-2016

Thank you. yes, i am getting the error at same line which you mentioned, how to fix the error

DanPatterson_Retired · ‎11-22-2016

Joshua, I think this thread follows on from this one https://community.esri.com/thread/181781-perform-idw-interpolation-to-select-multiple-z-field-from-a...

ShouvikJha · ‎11-22-2016

Dan , thank you. that link for multiple Z fields in one point feature class, but here i have multiple point feature class and each point feature class have one z field.

ShouvikJha · ‎11-22-2016

Hi, All Please provide your valuable input to fix the error. Thank you

DanPatterson_Retired · ‎11-22-2016

I don't see where you are cycling through a list of input files. All I see is your code cycling through a list of fields in a file. Is something missing? or did you just change the interpolator from the associated example? It appears you want someone to write the loop to cycle through a list of featureclasses or shapefiles, then pick out one field from it. You can emulate your previous example by just changing what you are cycling through and what is fixed. Perhaps if we could see that code rather than code based on your previous question, it might help.

ShouvikJha · ‎11-22-2016

Dan , i have made changes in line no 28, as i changed it ,

for fc in points:‍

but still error i am getting.

and error message is "object is not iterable".

DanPatterson_Retired · ‎11-22-2016

Sadly 'points' doesn't provide the necessary information. If you need to do scripting, I suggest you begin examining the code snippets in the Arcpy section What is ArcPy?—Help | ArcGIS for Desktop

beginning with ListFeatureClasses—Help | ArcGIS for Desktop to get some idea how to modify your codebase

ShouvikJha · ‎11-22-2016

I am trying in my best, i am posting here what i have so far understood,

can anyone suggest me how to iterate Friging function to all point data

points = arcpy.ListFeatureClasses('*','POINT')
for fc in points:
    path = os.path.join(arcpy.env.workspace, fc)
    Point = Point_Num[3:]‍‍‍‍