Perform Kriging interpolatation using ArcPy

6263
25
Jump to solution
11-22-2016 07:23 AM
ShouvikJha
Occasional Contributor III

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))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
25 Replies
XanderBakker
Esri Esteemed Contributor

... 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
Esri Esteemed Contributor

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

ShouvikJha
Occasional Contributor III

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)
0 Kudos
XanderBakker
Esri Esteemed Contributor

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
Occasional Contributor III

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. 

0 Kudos
XanderBakker
Esri Esteemed Contributor

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

0 Kudos