Perform Kriging interpolatation using ArcPy

4167
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
1 Solution

Accepted Solutions
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'

View solution in original post

25 Replies
JoshuaBixby
MVP Esteemed Contributor

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

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

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor
ShouvikJha
Occasional Contributor III

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. 

0 Kudos
ShouvikJha
Occasional Contributor III

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

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

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

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". 

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

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

0 Kudos
ShouvikJha
Occasional Contributor III

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:]
0 Kudos