topic Re: Perform Kriging interpolatation using ArcPy in Python Questions

Perform Kriging interpolatation using ArcPy

ShouvikJha — Sat, 11 Dec 2021 08:24:08 GMT

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

Re: Perform Kriging interpolatation using ArcPy

JoshuaBixby — Tue, 22 Nov 2016 15:49:51 GMT

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?

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Tue, 22 Nov 2016 16:00:19 GMT

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

Re: Perform Kriging interpolatation using ArcPy

DanPatterson_Retired — Tue, 22 Nov 2016 16:03:06 GMT

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-attribute-table-using-arcpy

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Tue, 22 Nov 2016 16:11:36 GMT

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.

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Tue, 22 Nov 2016 17:45:20 GMT

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

Re: Perform Kriging interpolatation using ArcPy

DanPatterson_Retired — Tue, 22 Nov 2016 17:55:40 GMT

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.

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Tue, 22 Nov 2016 19:22:11 GMT

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

Re: Perform Kriging interpolatation using ArcPy

DanPatterson_Retired — Tue, 22 Nov 2016 22:28:33 GMT

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

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Sat, 11 Dec 2021 08:24:11 GMT

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:]‍‍‍‍

Re: Perform Kriging interpolatation using ArcPy

curtvprice — Sat, 11 Dec 2021 08:24:13 GMT

In the script example you gave you aren't setting the workspace. If you don't do that you will not get a list of fc's, as ListFeatureClasses looks in the current env.workspace.

Once you are looping, you don't need to provide the full path, just the fc name from the list as the tools will look for the points fc in env.workspace.

Try this:

zField = "GRID_CODE"  # field in each input point fc

# get list of point fcs
env.workspace = r"D:\NPP_Test\MERRA_TEMP_2012C"
points = arcpy.ListFeatureClasses("r*mean", "POINT")

# where to save output
out_folder = r"D:\NPP_Test\MERRA_TEMP_2012D"

for fc in points:
    kModelUniversalObj = KrigingModelUniversal(
                             "LINEARDRIFT", lagSize, 
                             majorRange, partialSill, nugget)
    OutKriging = Kriging(fc, zField, kModelUniversalObj, cellSize, kRadius)
    KrigMask = ExtractByMask(OutKriging, mask)

    #Save output raster to output workspace
    tifname = fc[:4]  # captures "r001" for output tiff name
    # construct output path - "r001.tif" in out_folder
    tif_out = os.path.join(out_folder, "{}.tif".format(tifname))
    KrigMask.save(tif_out)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Sat, 11 Dec 2021 08:24:16 GMT

Hi, Curis, Thank you. I have updated my code according to you. and i did bit changes in line no 37.

Below is my updated script, this script neither producing any error nor output.

Also I have attached sample input file, which i trying to interpolate. Plase find the sample attached file for reference

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

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

zField = "GRID_CODE"  # field in each input point fc
#Kriging Veriable
cellSize = 0.05
lagSize = 0.5780481172534
majorRange = 6
partialSill = 3.304292110
nugget = 0.002701348
kRadius = RadiusFixed(20000, 1)

# get list of point fcs
env.workspace = r"D:\NPP_Test\KRIGING\Test_Points"
points = arcpy.ListFeatureClasses("r*mean", "POINT")

# where to save output
out_folder = r"D:\NPP_Test\KRIGING\OUTPUT_RAS"

#Mask region of interest
mask = r"D:\NPP_Test\KRIGING\MASK_FILE\GUJARAT-STATE_NEW.shp"


for fc in points:
    kModelUniversalObj = KrigingModelUniversal("LINEARDRIFT", lagSize, majorRange, partialSill, nugget)
    OutKriging = Kriging(fc, zField, kModelUniversalObj, cellSize, kRadius)
    KrigMask = ExtractByMask(OutKriging, mask)

    #Save output raster to output workspace
    tifname = fc[:4]  # captures "r001" for output tiff name
    Krigout = os.path.join(out_folder, 'r{0}.tif'.format(tifname))
    KrigMask.save(Krigout)

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Perform Kriging interpolatation using ArcPy

DanPatterson_Retired — Wed, 23 Nov 2016 08:27:04 GMT

When you get no output, no messages and no errors, it usually means it worked...but not the way you wanted to. It is then time for the print statements... like

line 22 print(points)

line 30.5 print(fc)

line 36.5 print(tifname)

line 37.5 print(krigout)

what were the results?

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Wed, 23 Nov 2016 08:42:40 GMT

Dan, I have added the print statements, but not any files shows while running the script.

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Wed, 23 Nov 2016 09:41:59 GMT

Xander Bakker‌, I have modified this code from your early code which was multiple Z fields and one feature class. I have tried hard to fix the issue suggested by many of people but unfortunately I am still stuck yet. Kindly cooperate with us to solve the issue. Thank you

Re: Perform Kriging interpolatation using ArcPy

DanPatterson_Retired — Wed, 23 Nov 2016 10:55:02 GMT

could you show the results with the print statements? Did the print(points) line print anything? etcetera

Re: Perform Kriging interpolatation using ArcPy

XanderBakker — Wed, 23 Nov 2016 11:37:21 GMT

It seems you are not getting any results when you list your featureclasses (the list points is empty). That is why you don't get any error and nothing seems to happen. As Dan Patterson‌ mentioned, it is good to include some print statements to verify the process or debug the code and check the content of variables along the way.

Since you list featureclasses, your list points contains:

r001_mean.shp
r002_mean.shp
r003_mean.shp

...

You should use:

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

Re: Perform Kriging interpolatation using ArcPy

XanderBakker — Wed, 23 Nov 2016 12:15:01 GMT

However, when I execute the code it throws an error:

ExecuteError: ERROR 010079: Unable to estimate semivariogram.

Since you are using point featureclasses with a GRID_CODE field, may I assume that the point featureclasses where created from rasters? If so, why use Kriging to turn them back into rasters?

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Sat, 11 Dec 2021 08:24:19 GMT

Xander Bakker‌ , Thank you very much.

Yes, you are absolutely correct. i have generated that point from raster data, actually i want to downscale the raster data to meet the spatial resolution with base map . Thats why i am trying to apply interpolation method.

But from Kriging method we getting error, possibly error due to downscaling cell size, I am thinking to apply linear interpolation IDW method, though i have modified that code according to IDW interpolation, but the problem is same with the script, neither error nor output.

Below is updated script for IDW

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

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

env.workspace = r"D:\NPP_Test\KRIGING\Test_Points"
points = arcpy.ListFeatureClasses("r*_mean", "POINT")
print points

out_folder = r"D:\NPP_Test\KRIGING\OUTPUT_RAS"

zFields = ["GRID_CODE"]

cellSize = 0.002298707671
power = 2
searchRadius = RadiusVariable(10, 150000)
#Mask region of interest
mask="D:\NPP_Test\KRIGING\MASK_FILE\GUJARAT-STATE_NEW.shp"

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

Re: Perform Kriging interpolatation using ArcPy

ShouvikJha — Thu, 24 Nov 2016 07:11:05 GMT

Dan , i have added the print statement but nothing file showing during program execution, i have posted below my full code with print statement. thank you