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

curtvprice · ‎11-22-2016

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

ShouvikJha · ‎11-22-2016

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)

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

DanPatterson_Retired · ‎11-23-2016

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?

ShouvikJha · ‎11-23-2016

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

DanPatterson_Retired · ‎11-23-2016

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

ShouvikJha · ‎11-23-2016

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

ShouvikJha · ‎11-23-2016

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

XanderBakker · ‎11-23-2016

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

XanderBakker · ‎11-23-2016

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?

ShouvikJha · ‎11-23-2016

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