IDW in python script

[RoyHewitt]

Hello All,

I'm working on a script that creates an interpolated raster surface of wetland vegetation.  Data points were collected for coastal marshes in boat surveys.  Each point contains fields for each of the 15 different types of vegetation documented (as percent cover).

The rivers are disjoint, and one river's vegetation should not affect the others (because of variables including salinity, tides, etc.).  Each set of points has been separated into "ModelUnits" of which there are 8.

My script has a nested loop: First, it loops through the study area, then loops through each vegetation type to create an interpolated surface (IDW).

Everything seems to work until the script reaches the IDW tool:  There are several IDW tools in different toolboxes (3D, Geostatistical Analysis, and Spatial Analysis).  I have chosen the spatial analysis version, but I'm confused as to syntax.  The help documents have only:  arcpy.IDW(parameters) as opposed to: arcpy.IDW_ga(parameters).

Trying IDW in the ArcMap python window, there is no option for the spatial analyst version of IDW.

I've attached a map image to give you an idea of what the files look like and my python script.

Any and all help is appreciated!

import arcpy
arcpy.env.overwriteOutput = True

# Create file paths
barrier = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/NutriaProj.mdb/WetlandBarrier/"
rhaPoints = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/NutriaProj.mdb/RHA_Waypoints"
habitatData = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/NutriaProj.mdb/habitatDetails"
rasterOutput = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/NutriaProj.mdb/"

# Create field variables
rhaField = "IDENT"
habitatField = "PointName"
rhaUnit = "ModelUnit"
barrierUnit = "ModelUnit"

# Create List Variables
fieldList = []
unitList = [1,2,3,4,5,6,7,8]

# Create join between RHA points and Habitat data table
try:
    arcpy.AddJoin_management(rhaPoints, rhaField, habitatData, habitatField)
except:
    print arcpy.GetMessages(2)

# Create list of vegetation attributes for loop, local variables
for field in arcpy.ListFields(habitatData):
    fieldList.append(field.name)
fList = fieldList[2:17]
    
power = 2  # As distance increases, point has less impact on interpolation
cellSize = 60  # Raster cell size

try:
    for unit in unitList:
        # Create where clause for current study area.
        whereRHA = '['+ rhaUnit + '] = ' + "'" + str(unitList[unit - 1]) + "'"
        whereBarrier = '['+ barrierUnit + '] = ' + "'" + str(unitList[unit - 1]) + "'"

        # Create feature layer for RHA waypoints and Barrier for current study unit
        arcpy.MakeFeatureLayer_management(rhaPoints, "currentRHA", whereRHA)
        arcpy.MakeFeatureLayer_management(barrier, "currentBarrier", whereBarrier)
        
        for feat in fList:
            try:
                # IDW geostat analyst
                #arcpy.Idw_ga("currentRHA", feat, "", rasterOutpath + feat + "_" + str(unit),cellSize, power, "", "currentBarrier")

                # IDW spatial analyst
                outRaster = arcpy.Idw_sa("currentRHA", feat, cellSize, 2, "", "currentBarrier")
                outRaster.save(rasterOutput + feat + "_" + str(unit))
                print "Successfully ran IDW for " + feat + "."
            except:
                print arcpy.GetMessages(2)
        # Delete feature layers
        arcpy.Delete_management("currentRHA")
        arcpy.Delete_management("currentBarrier")
except:
    print arcpy.GetMessages(2)

[JakeSkinner]

Hi Roy,

You'll need to checkout the spatial analyst extension first:

arcpy.CheckOutExtension("Spatial")

Then you can use 'arcpy.sa.Idw'.  Take a look here for further examples.

[RoyHewitt]

I added the check out / in methods to my script, but I don't think that's the problem (I'm loading and running the script in the python window with the extension checked in ArcMap.

I suspect the problem is with the field names.  Since the fields that I'm using for the z-value are from a joined table, the IDW isn't finding them.

For the Z-value field I've tried:
1.  Looping through Python list of fields (ex. SaltmarshBulrushPercentage)
2.  Adding name of joined table (ex. habitatDetails.SaltmarshBulrushPercentage)

Both of these methods return the same error:
ERROR 001000: Z value field: Field SaltmarshBulrushPercentage does not exist


RESOLVED Error 001000
AddJoin does not work with a feature class, only a feature layer.  I rewrote the top of the script to get the join to work, now the script recognizes the field name.

[RoyHewitt]

Updated Script and Error:


Script runs correctly up to the point of IDW.
Error: 010092: Invalid output extent.
 

#  Author: Roy Hewitt
#  US Fish and Wildlife Service
#  December 2011

# Import modules, setup overwrite in environments
import arcpy
from arcpy.sa import *

arcpy.env.overwriteOutput = True

# Create file paths
barrier = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/NutriaProj.mdb/WetlandBarrier/"
rhaPoints = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/NutriaProj.mdb/RHA_Waypoints"
habitatData = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/NutriaProj.mdb/habitatDetails"
outPath = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/"
rasterPath = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/Rasters/"
rhaLayer = "P:/Employee_GIS_Data/Hewitt_GIS_Data/Python/RHA_Waypoints.lyr"

# Create field variables
rhaField = "IDENT"
habitatField = "PointName"
rhaUnit = "ModelUnit"
barrierUnit = "ModelUnit"

# Create List Variables
fieldList = []
unitList = [1,2,3,4,5,6,7,8]

# Create join between RHA points and Habitat data table
try:
    print "Joining waypoints with habitat data..."
    arcpy.MakeFeatureLayer_management(rhaPoints, "rhaLyr")
    arcpy.AddJoin_management("rhaLyr", rhaField, habitatData, habitatField)
    arcpy.SaveToLayerFile_management("rhaLyr", rhaLayer)
except:
    print arcpy.GetMessages(0)

# Create list of vegetation attributes for loop, local variables
print "Creating list of vegetation fields for IDW modeling..."
for field in arcpy.ListFields(habitatData):
    fieldList.append(field.name)
fList = fieldList[2:17]
    
power = 2  # As distance increases, point has less impact on interpolation
cellSize = 60  # Raster cell size

# Check out Spatial Analyst extension
try:
    if arcpy.CheckExtension("spatial")== "Available":
        arcpy.CheckOutExtension("spatial")
        print "Spatial license checked out."
except:
    print "Spatial Analyst extension not available."
    print arcpy.GetMessages(2)
    
try:
    for unit in unitList:
        # Create where clause for current study area.
        print "Creating where clause for study unit " + str(unit) + "..."
        whereRHA = '['+ rhaUnit + '] = ' + "'" + str(unitList[unit - 1]) + "'"
        whereBarrier = '['+ barrierUnit + '] = ' + "'" + str(unitList[unit - 1]) + "'"

        # Create feature layer for RHA waypoints and Barrier for current study unit
        arcpy.MakeFeatureLayer_management("rhaLyr", "currentRHA", whereRHA)
        arcpy.MakeFeatureLayer_management(barrier, "currentBarrier", whereBarrier)
        
        for feat in fList:
            try:
                print "Running IDW model for " + feat + "..."
                # IDW geostat analyst
                #arcpy.Idw_ga("currentRHA", feat, "", rasterOutpath + feat + "_" + str(unit),cellSize, power, "", "currentBarrier")

                # IDW spatial analyst
                outRaster = arcpy.sa.Idw("currentRHA", "habitatDetails." + feat, cellSize, power, "", "currentBarrier")
                outRaster.save(outRaster + feat + "_" + str(unit))
                print "Successfully ran IDW for " + feat + "."
            except:
                print feat + " interpolation failed."
                print arcpy.GetMessages(2)

        # Delete feature layer
        arcpy.Delete_management("currentRHA")
        arcpy.Delete_management("currentBarrier")
except:
    print arcpy.GetMessages(2)
finally:
    arcpy.CheckInExtension("spatial")
    print "Spatial license checked back in."

print "Finished running IDW module."
                
                

[CarrieDavis]

Hi Roy,

Just a couple trouble shooting suggestions from looking real quick at your script.

Try importing the arcpy environments settings at the beginning of your script like below.

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *

Are you able to run the IDW tool from the toolbox?

If you enter your data into the format of the sample IDW stand alone python script given at http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009z0000006m000000.htm.  Does this work?

[TonyPagani]

Hi Roy,

In the script it looks like you are using a a barrier in the idw. Does the barrier feature class have the same projected coordinate system as the point feature class?  If they are different you may get the error you are seeing. You should keep all the data in your analysis in the same projected coordinate system when doing spatial analysis to get accurate results.

Tony