roy_hewitt

IDW in python script

Discussion created by roy_hewitt on Jan 19, 2012
Latest reply on Jan 25, 2012 by TPagani-esristaff
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)

Attachments

Outcomes