Perform Kriging interpolatation using ArcPy

[ShouvikJha]

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]

... that because on line 11 you are still producing an empty list ("points") of point featureclasses. See me comment before and change the wildcard into "r*_mean.shp"

[XanderBakker]

Also on line 20 the maximum distance of 15000 map units (I assume they are degrees) is "rather" big. 

[ShouvikJha]

Xander Bakker‌, Thank you. I ran the script by replacing ("r*_mean.shp"). Still i am getting error. 

Traceback (most recent call last):
  File "<module1>", line 27, in <module>
  File "C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 2467, in Idw
    in_barrier_polyline_features)
  File "C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py", line 2459, in Wrapper
    in_barrier_polyline_features)
  File "C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\geoprocessing\_base.py", line 504, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
RuntimeError: Object: Error in executing tool‍‍‍‍‍‍‍‍‍‍‍

Line no 20, i have followed the Arc GIS default value. 

searchRadius = RadiusVariable(10, 150000)

[XanderBakker · Accepted solution]

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'‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[ShouvikJha]

Xander Bakker‌, Thank you very much for your kind cooperation. your script working absolutely correct. 

 I will readjust the search RadiusVariable accordingly. Really i learned a totally new thing from you about the importance of RadiusVariable. Earlier i just kept RadiusVariable as a default value. 

[XanderBakker]

You're welcome. I'm glad it works!