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

nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

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"

<SPAN class="keyword token">import</SPAN> arcpy
<SPAN class="keyword token">from</SPAN> arcpy <SPAN class="keyword token">import</SPAN> env
<SPAN class="keyword token">from</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>sa <SPAN class="keyword token">import</SPAN> <SPAN class="operator token">*</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>overwriteOutput <SPAN class="operator token">=</SPAN> <SPAN class="token boolean">True</SPAN>

<SPAN class="comment token"># Check out the ArcGIS Spatial Analyst extension license</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>CheckOutExtension<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"Spatial"</SPAN><SPAN class="punctuation token">)</SPAN>

In_Point <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\MERRA_TEMP_2012C"</SPAN> <SPAN class="comment token">#(Point feature name:r001_mean, r002_mean.....r012_mean )</SPAN>
Out_Raster <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\MERRA_TEMP_2012D"</SPAN>

points <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListFeatureClasses<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">)</SPAN>


zFields <SPAN class="operator token">=</SPAN> <SPAN class="string token">"GRID_CODE"</SPAN>

<SPAN class="comment token">#Kriging Veriable</SPAN>
cellSize <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.05</SPAN>
lagSize <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.5780481172534</SPAN>
majorRange <SPAN class="operator token">=</SPAN> <SPAN class="number token">6</SPAN>
partialSill <SPAN class="operator token">=</SPAN> <SPAN class="number token">3.304292110</SPAN>
nugget <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.002701348</SPAN>
kRadius <SPAN class="operator token">=</SPAN> RadiusFixed<SPAN class="punctuation token">(</SPAN><SPAN class="number token">20000</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">1</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token">#Mask region of interest</SPAN>
mask <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\Gujarta Shape file\GUJARATSTATE.shp"</SPAN>

<SPAN class="keyword token">for</SPAN> zField <SPAN class="keyword token">in</SPAN> zFields<SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; Point <SPAN class="operator token">=</SPAN> Point_Num<SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="number token">3</SPAN><SPAN class="punctuation token">]</SPAN>
&nbsp;&nbsp;&nbsp; kModelUniversalObj <SPAN class="operator token">=</SPAN> KrigingModelUniversal<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"LINEARDRIFT"</SPAN><SPAN class="punctuation token">,</SPAN> lagSize<SPAN class="punctuation token">,</SPAN> majorRange<SPAN class="punctuation token">,</SPAN> partialSill<SPAN class="punctuation token">,</SPAN> nugget<SPAN class="punctuation token">)</SPAN>

&nbsp;&nbsp;&nbsp; OutKriging <SPAN class="operator token">=</SPAN> Kriging<SPAN class="punctuation token">(</SPAN>inPointFeatures<SPAN class="punctuation token">,</SPAN> zField<SPAN class="punctuation token">,</SPAN> kModelUniversalObj<SPAN class="punctuation token">,</SPAN> cellSize<SPAN class="punctuation token">,</SPAN> kRadius<SPAN class="punctuation token">)</SPAN>


&nbsp;&nbsp;&nbsp; <SPAN class="comment token">#IDWMASk = ExtractByMask(outIDW, mask)</SPAN>
&nbsp;&nbsp;&nbsp; KrigMask <SPAN class="operator token">=</SPAN> ExtractByMask<SPAN class="punctuation token">(</SPAN>OutKriging<SPAN class="punctuation token">,</SPAN> mask<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="comment token">#Save outraster as the same name of input</SPAN>
&nbsp;&nbsp;&nbsp; KrigMask<SPAN class="punctuation token">.</SPAN>save<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"r{}.tif"</SPAN><SPAN class="punctuation token">.</SPAN>format<SPAN class="punctuation token">(</SPAN>Point_Num<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

Vind meer berichten getagd met

Accepted answers

XanderBakker

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:

<SPAN class="keyword token">import</SPAN> arcpy
<SPAN class="keyword token">from</SPAN> arcpy <SPAN class="keyword token">import</SPAN> env
<SPAN class="keyword token">from</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>sa <SPAN class="keyword token">import</SPAN> <SPAN class="operator token">*</SPAN>
<SPAN class="keyword token">import</SPAN> os
arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>parallelProcessingFactor <SPAN class="operator token">=</SPAN> <SPAN class="string token">"100%"</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>overwriteOutput <SPAN class="operator token">=</SPAN> <SPAN class="token boolean">True</SPAN>

<SPAN class="comment token"># settings</SPAN>
env<SPAN class="punctuation token">.</SPAN>workspace <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\Xander\GeoNet\KrigingMultiFC\KRIGING\Test_Points"</SPAN>
mask <SPAN class="operator token">=</SPAN> <SPAN class="string token">"D:\Xander\GeoNet\KrigingMultiFC\KRIGING\MASK_FILE\GUJARAT-STATE_NEW.shp"</SPAN>
out_folder <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\Xander\GeoNet\KrigingMultiFC\KRIGING\OUTPUT_RAS"</SPAN>
zField <SPAN class="operator token">=</SPAN> <SPAN class="string token">"GRID_CODE"</SPAN> <SPAN class="comment token"># has to be string, not list</SPAN>

points <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListFeatureClasses<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"r*_mean.shp"</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">"POINT"</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">print</SPAN> points

cellSize <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.025</SPAN> <SPAN class="comment token"># 0.002298707671</SPAN>
power <SPAN class="operator token">=</SPAN> <SPAN class="number token">2</SPAN>
searchRadius <SPAN class="operator token">=</SPAN> RadiusVariable<SPAN class="punctuation token">(</SPAN><SPAN class="number token">10</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">5</SPAN><SPAN class="punctuation token">)</SPAN> <SPAN class="comment token"># smaller max distance</SPAN>

<SPAN class="comment token"># Check out the ArcGIS Spatial Analyst extension license</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>CheckOutExtension<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"Spatial"</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token"># Execute IDW</SPAN>
<SPAN class="keyword token">for</SPAN> fc <SPAN class="keyword token">in</SPAN> points<SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="keyword token">print</SPAN> fc
&nbsp;&nbsp;&nbsp; outIDW <SPAN class="operator token">=</SPAN> Idw<SPAN class="punctuation token">(</SPAN>fc<SPAN class="punctuation token">,</SPAN> zField<SPAN class="punctuation token">,</SPAN> cellSize<SPAN class="punctuation token">,</SPAN> power<SPAN class="punctuation token">,</SPAN> searchRadius<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># Execute Mask</SPAN>
&nbsp;&nbsp;&nbsp; IDWMASk <SPAN class="operator token">=</SPAN> ExtractByMask<SPAN class="punctuation token">(</SPAN>outIDW<SPAN class="punctuation token">,</SPAN> mask<SPAN class="punctuation token">)</SPAN>

&nbsp;&nbsp;&nbsp; <SPAN class="comment token">#Save output raster to output workspace</SPAN>
&nbsp;&nbsp;&nbsp; tifname <SPAN class="operator token">=</SPAN> fc<SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="number token">4</SPAN><SPAN class="punctuation token">]</SPAN>&nbsp; <SPAN class="comment token"># captures "r001" for output tiff name</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="keyword token">print</SPAN> tifname

&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># Save output, except Solar Radiation raster</SPAN>
&nbsp;&nbsp;&nbsp; IDW_OUT<SPAN class="operator token">=</SPAN> os<SPAN class="punctuation token">.</SPAN>path<SPAN class="punctuation token">.</SPAN>join<SPAN class="punctuation token">(</SPAN>out_folder<SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'r{0}.tif'</SPAN><SPAN class="punctuation token">.</SPAN>format<SPAN class="punctuation token">(</SPAN>tifname<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; IDWMASk<SPAN class="punctuation token">.</SPAN>save<SPAN class="punctuation token">(</SPAN>IDW_OUT<SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token">#print done</SPAN>
<SPAN class="keyword token">print</SPAN> <SPAN class="string token">'done'</SPAN><SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

All comments

XanderBakker

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

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"

ShouvikJha

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

<SPAN class="keyword token">import</SPAN> arcpy
<SPAN class="keyword token">from</SPAN> arcpy <SPAN class="keyword token">import</SPAN> env
<SPAN class="keyword token">from</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>sa <SPAN class="keyword token">import</SPAN> <SPAN class="operator token">*</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>parallelProcessingFactor <SPAN class="operator token">=</SPAN> <SPAN class="string token">"100%"</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>overwriteOutput <SPAN class="operator token">=</SPAN> <SPAN class="token boolean">True</SPAN>

<SPAN class="comment token"># Check out the ArcGIS Spatial Analyst extension license</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>CheckOutExtension<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"Spatial"</SPAN><SPAN class="punctuation token">)</SPAN>

env<SPAN class="punctuation token">.</SPAN>workspace <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\KRIGING\Test_Points"</SPAN>
points <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListFeatureClasses<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"r*_mean"</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">"POINT"</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">print</SPAN> points

out_folder <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\KRIGING\OUTPUT_RAS"</SPAN>

zFields <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="string token">"GRID_CODE"</SPAN><SPAN class="punctuation token">]</SPAN>

cellSize <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.002298707671</SPAN>
power <SPAN class="operator token">=</SPAN> <SPAN class="number token">2</SPAN>
searchRadius <SPAN class="operator token">=</SPAN> RadiusVariable<SPAN class="punctuation token">(</SPAN><SPAN class="number token">10</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">150000</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="comment token">#Mask region of interest</SPAN>
mask<SPAN class="operator token">=</SPAN><SPAN class="string token">"D:\NPP_Test\KRIGING\MASK_FILE\GUJARAT-STATE_NEW.shp"</SPAN>

<SPAN class="comment token"># Execute IDW</SPAN>
<SPAN class="keyword token">for</SPAN> fc <SPAN class="keyword token">in</SPAN> points<SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="keyword token">print</SPAN> fc
&nbsp;&nbsp;&nbsp; outIDW <SPAN class="operator token">=</SPAN> Idw<SPAN class="punctuation token">(</SPAN>fc<SPAN class="punctuation token">,</SPAN> zField<SPAN class="punctuation token">,</SPAN> cellSize<SPAN class="punctuation token">,</SPAN> power<SPAN class="punctuation token">,</SPAN> searchRadius<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># Execute Mask</SPAN>
&nbsp;&nbsp;&nbsp; IDWMASk <SPAN class="operator token">=</SPAN> ExtractByMask<SPAN class="punctuation token">(</SPAN>outIDW<SPAN class="punctuation token">,</SPAN> mask<SPAN class="punctuation token">)</SPAN>

&nbsp;&nbsp;&nbsp; <SPAN class="comment token">#Save output raster to output workspace</SPAN>
&nbsp;&nbsp;&nbsp; tifname <SPAN class="operator token">=</SPAN> fc<SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="number token">4</SPAN><SPAN class="punctuation token">]</SPAN>&nbsp; <SPAN class="comment token"># captures "r001" for output tiff name</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="keyword token">print</SPAN> tifname

&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># Save output, except Solar Radiation raster</SPAN>
&nbsp;&nbsp;&nbsp; IDW_OUT<SPAN class="operator token">=</SPAN> os<SPAN class="punctuation token">.</SPAN>path<SPAN class="punctuation token">.</SPAN>join<SPAN class="punctuation token">(</SPAN>out_folder<SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'r{0}.tif'</SPAN><SPAN class="punctuation token">.</SPAN>format<SPAN class="punctuation token">(</SPAN>tifname<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; IDWMASk<SPAN class="punctuation token">.</SPAN>save<SPAN class="punctuation token">(</SPAN>IDW_OUT<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; 

<SPAN class="comment token">#print done</SPAN>
<SPAN class="keyword token">print</SPAN> <SPAN class="string token">'done'</SPAN><SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

XanderBakker

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?

XanderBakker

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 <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListFeatureClasses<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"r*_mean.shp"</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">"POINT"</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="line-numbers-rows"><SPAN>‍</SPAN></SPAN>

ShouvikJha

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

DanPatterson_Retired

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

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

<SPAN class="keyword token">import</SPAN> arcpy
<SPAN class="keyword token">import</SPAN> os
<SPAN class="keyword token">from</SPAN> arcpy <SPAN class="keyword token">import</SPAN> env
<SPAN class="keyword token">from</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>sa <SPAN class="keyword token">import</SPAN> <SPAN class="operator token">*</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>overwriteOutput <SPAN class="operator token">=</SPAN> <SPAN class="token boolean">True</SPAN>

<SPAN class="comment token"># Check out the ArcGIS Spatial Analyst extension license</SPAN>
arcpy<SPAN class="punctuation token">.</SPAN>CheckOutExtension<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"Spatial"</SPAN><SPAN class="punctuation token">)</SPAN>

zField <SPAN class="operator token">=</SPAN> <SPAN class="string token">"GRID_CODE"</SPAN>&nbsp; <SPAN class="comment token"># field in each input point fc</SPAN>
<SPAN class="comment token">#Kriging Veriable</SPAN>
cellSize <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.05</SPAN>
lagSize <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.5780481172534</SPAN>
majorRange <SPAN class="operator token">=</SPAN> <SPAN class="number token">6</SPAN>
partialSill <SPAN class="operator token">=</SPAN> <SPAN class="number token">3.304292110</SPAN>
nugget <SPAN class="operator token">=</SPAN> <SPAN class="number token">0.002701348</SPAN>
kRadius <SPAN class="operator token">=</SPAN> RadiusFixed<SPAN class="punctuation token">(</SPAN><SPAN class="number token">20000</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">1</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token"># get list of point fcs</SPAN>
env<SPAN class="punctuation token">.</SPAN>workspace <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\KRIGING\Test_Points"</SPAN>
points <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListFeatureClasses<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"r*mean"</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">"POINT"</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token"># where to save output</SPAN>
out_folder <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\KRIGING\OUTPUT_RAS"</SPAN>

<SPAN class="comment token">#Mask region of interest</SPAN>
mask <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\KRIGING\MASK_FILE\GUJARAT-STATE_NEW.shp"</SPAN>


<SPAN class="keyword token">for</SPAN> fc <SPAN class="keyword token">in</SPAN> points<SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; kModelUniversalObj <SPAN class="operator token">=</SPAN> KrigingModelUniversal<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"LINEARDRIFT"</SPAN><SPAN class="punctuation token">,</SPAN> lagSize<SPAN class="punctuation token">,</SPAN> majorRange<SPAN class="punctuation token">,</SPAN> partialSill<SPAN class="punctuation token">,</SPAN> nugget<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; OutKriging <SPAN class="operator token">=</SPAN> Kriging<SPAN class="punctuation token">(</SPAN>fc<SPAN class="punctuation token">,</SPAN> zField<SPAN class="punctuation token">,</SPAN> kModelUniversalObj<SPAN class="punctuation token">,</SPAN> cellSize<SPAN class="punctuation token">,</SPAN> kRadius<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; KrigMask <SPAN class="operator token">=</SPAN> ExtractByMask<SPAN class="punctuation token">(</SPAN>OutKriging<SPAN class="punctuation token">,</SPAN> mask<SPAN class="punctuation token">)</SPAN>

&nbsp;&nbsp;&nbsp; <SPAN class="comment token">#Save output raster to output workspace</SPAN>
&nbsp;&nbsp;&nbsp; tifname <SPAN class="operator token">=</SPAN> fc<SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="number token">4</SPAN><SPAN class="punctuation token">]</SPAN>&nbsp; <SPAN class="comment token"># captures "r001" for output tiff name</SPAN>
&nbsp;&nbsp;&nbsp; Krigout <SPAN class="operator token">=</SPAN> os<SPAN class="punctuation token">.</SPAN>path<SPAN class="punctuation token">.</SPAN>join<SPAN class="punctuation token">(</SPAN>out_folder<SPAN class="punctuation token">,</SPAN> <SPAN class="string token">'r{0}.tif'</SPAN><SPAN class="punctuation token">.</SPAN>format<SPAN class="punctuation token">(</SPAN>tifname<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; KrigMask<SPAN class="punctuation token">.</SPAN>save<SPAN class="punctuation token">(</SPAN>Krigout<SPAN class="punctuation token">)</SPAN>

‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

KRIGING.rar

curtvprice

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 <SPAN class="operator token">=</SPAN> <SPAN class="string token">"GRID_CODE"</SPAN>&nbsp; <SPAN class="comment token"># field in each input point fc</SPAN>

<SPAN class="comment token"># get list of point fcs</SPAN>
env<SPAN class="punctuation token">.</SPAN>workspace <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\MERRA_TEMP_2012C"</SPAN>
points <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListFeatureClasses<SPAN class="punctuation token">(</SPAN><SPAN class="string token">"r*mean"</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="string token">"POINT"</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="comment token"># where to save output</SPAN>
out_folder <SPAN class="operator token">=</SPAN> r<SPAN class="string token">"D:\NPP_Test\MERRA_TEMP_2012D"</SPAN>

<SPAN class="keyword token">for</SPAN> fc <SPAN class="keyword token">in</SPAN> points<SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; kModelUniversalObj <SPAN class="operator token">=</SPAN> KrigingModelUniversal<SPAN class="punctuation token">(</SPAN>
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; <SPAN class="string token">"LINEARDRIFT"</SPAN><SPAN class="punctuation token">,</SPAN> lagSize<SPAN class="punctuation token">,</SPAN> 
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; majorRange<SPAN class="punctuation token">,</SPAN> partialSill<SPAN class="punctuation token">,</SPAN> nugget<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; OutKriging <SPAN class="operator token">=</SPAN> Kriging<SPAN class="punctuation token">(</SPAN>fc<SPAN class="punctuation token">,</SPAN> zField<SPAN class="punctuation token">,</SPAN> kModelUniversalObj<SPAN class="punctuation token">,</SPAN> cellSize<SPAN class="punctuation token">,</SPAN> kRadius<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; KrigMask <SPAN class="operator token">=</SPAN> ExtractByMask<SPAN class="punctuation token">(</SPAN>OutKriging<SPAN class="punctuation token">,</SPAN> mask<SPAN class="punctuation token">)</SPAN>

&nbsp;&nbsp;&nbsp; <SPAN class="comment token">#Save output raster to output workspace</SPAN>
&nbsp;&nbsp;&nbsp; tifname <SPAN class="operator token">=</SPAN> fc<SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="number token">4</SPAN><SPAN class="punctuation token">]</SPAN>&nbsp; <SPAN class="comment token"># captures "r001" for output tiff name</SPAN>
&nbsp;&nbsp;&nbsp; <SPAN class="comment token"># construct output path - "r001.tif" in out_folder</SPAN>
&nbsp;&nbsp;&nbsp; tif_out <SPAN class="operator token">=</SPAN> os<SPAN class="punctuation token">.</SPAN>path<SPAN class="punctuation token">.</SPAN>join<SPAN class="punctuation token">(</SPAN>out_folder<SPAN class="punctuation token">,</SPAN> <SPAN class="string token">"{}.tif"</SPAN><SPAN class="punctuation token">.</SPAN>format<SPAN class="punctuation token">(</SPAN>tifname<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; KrigMask<SPAN class="punctuation token">.</SPAN>save<SPAN class="punctuation token">(</SPAN>tif_out<SPAN class="punctuation token">)</SPAN><SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

ShouvikJha

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 <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListFeatureClasses<SPAN class="punctuation token">(</SPAN><SPAN class="string token">'*'</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="string token">'POINT'</SPAN><SPAN class="punctuation token">)</SPAN>
<SPAN class="keyword token">for</SPAN> fc <SPAN class="keyword token">in</SPAN> points<SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; path <SPAN class="operator token">=</SPAN> os<SPAN class="punctuation token">.</SPAN>path<SPAN class="punctuation token">.</SPAN>join<SPAN class="punctuation token">(</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>env<SPAN class="punctuation token">.</SPAN>workspace<SPAN class="punctuation token">,</SPAN> fc<SPAN class="punctuation token">)</SPAN>
&nbsp;&nbsp;&nbsp; Point <SPAN class="operator token">=</SPAN> Point_Num<SPAN class="punctuation token">[</SPAN><SPAN class="number token">3</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

DanPatterson_Retired

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

ShouvikJha

Dan , i have made changes in line no 28, as i changed it ,

<SPAN class="keyword token">for</SPAN> fc <SPAN class="keyword token">in</SPAN> points<SPAN class="punctuation token">:</SPAN><SPAN class="line-numbers-rows"><SPAN>‍</SPAN></SPAN>

but still error i am getting.

and error message is "object is not iterable".

DanPatterson_Retired

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

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

ShouvikJha

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.

DanPatterson_Retired

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

ShouvikJha

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

JoshuaBixby

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?

Doe mee aan het gesprek

Aangemelde leden kunnen berichten plaatsen, updates volgen en meer. Nieuw hier? Registreer een gratis account.

Snelle Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.