topic Re: raster clipping with overlapping polygons in Python Questions

raster clipping with overlapping polygons

matseymour — Thu, 16 Aug 2012 11:04:47 GMT

Hello,
I am quite new to using python, but...
I have a set of overlapping polygons that I would like to use to clip several raster layers

I have the following code thus far

import arcpy # set environments arcpy.env.workspace = "Path" arcpy.env.overwriteOutput = True #set parameters featureTable="Individual_catchments" baseFeature="elevation" #set cursor rows = arcpy.SearchCursor(featureClass) # loop  for row in rows:      #here I select the geometry       feat = row.Shape      arcpy.Clip_management(baseFeature, "#", "catch"+str(row.Local)+".tif",feat, "0", "ClippingGeometry")

however, when I try to run I get an error stating "ERROR 999999: Error executing function. Raster dataset already exists Failed to execute (Clip). "

thank you for your time and help

cheers

Mat

Re: raster clipping with overlapping polygons

ChristopherThompson — Sun, 12 Dec 2021 08:27:13 GMT

A couple of things that come to mind but i haven't tested:

try building your output raster name outside of the tool, and pass it in as a variable, putting those sorts of expressions in a tool can be confusing to read and lead to syntax errors

If the row.Local field contains duplicate values, you could indeed be trying to overwrite an existing raster that has been created, so you might want to double check that.

you're passing in the geometry of a polygon as the 'feat' variable - if you read the help on this tool i think it wants a feature layer or another raster layer instead of the geometry. So create a feature layer based on the current row you are processing like this:

feat_id = row.FID #or however you want to uniquely identify a given polygon)
qry = """ "FID" = """ + feat_id
feat = arcpy.MakeFeatureLayer(featureClass,'feat_lyr',qry)
out_rast = "catch" + str(row.Local) + ".tif"
arcpy.Clip_management(baseFeature, "#", out_rast,feat, "0", "ClippingGeometry")

The MakeFeatureLayer makes an in-memory layer with a query definition based on the current row's FID and then is passed in as an argument to the Clip tool.

Re: raster clipping with overlapping polygons

MathewCoyle — Sun, 12 Dec 2021 08:27:15 GMT

Further to what Christopher said, you want to extract your feature class data of interest using your cursor, and then access that outside the cursor. Here is an example using a dictionary.

import arcpy
# set environments
arcpy.env.workspace = "Path"
arcpy.env.overwriteOutput = True
#set parameters
featureTable="Individual_catchments"
baseFeature="elevation"
temp_featureTable = r"in_memory\temp1"
temp_baseFeature = r"in_memory\temp2"
temp_featureClass = r"in_memory\temp_fc"
temp_layer = r"in_memory\temp"

arcpy.MakeFeatureLayer_management(featureTable, temp_featureTable)
arcpy.MakeFeatureLayer_management(baseFeature, temp_baseFeature)
arcpy.MakeFeatureLayer_management(featureClass, temp_featureClass)

# Get OID field name
f_list = arcpy.ListFields(temp_featureClass)
for f in f_list:
    if f.type == "OID":
        f_oid = f.name

# Set Dictionary
feat_dict = {}

#set cursor
rows = arcpy.SearchCursor(temp_featureClass)

# loop
for row in rows:
     feat_dict[row.getValue(f_oid)] = row.getValue("Local")

for id in feat_dict:
    arcpy.MakeFeatureLayer_management(temp_featureClass, temp_layer, "%s = %s" % (f_oid, id))
    arcpy.Clip_management(temp_baseFeature, "#", "catch%s.tif" % (feat_dict[id]), temp_layer, "0", "ClippingGeometry")

Re: raster clipping with overlapping polygons

matseymour — Wed, 22 Aug 2012 09:31:31 GMT

Thanks for the input folks. After several instances of trial and error I managed to get the below script to produce the output I needed, and feel much more confident in using python to retrieve data. I am sure the code could be simplified significantly, but it is enough for my simple mind. Of course any suggestions would be welcomed.
cheers
Mat

###this first part creates individual raster files (from as single larger raster file) based on overlapping polygons import math import os import arcpy # set the workspace environment  arcpy.env.workspace = "D:/path.mdb" arcpy.env.overwriteOutput = True # set the cliping regions featureClass="Individual_catchments" # set the base feature baseFeature="elevation_swiss" #create list to store catchment elevations catchmentList = [] #create table to store catchment id and average elevation CreateTable_management ("D:/path.mdb", "Average_elevations")   # need to use a cursor argument since we are looking at a single table # this one assigns a search cursor to the Individual catchments table rows = arcpy.SearchCursor(featureClass)   # loop to create all of my individual files                            for row in rows:      #here we only want the geometry (shape)      feat = row.Shape      output="c"+str(row.LOCALITYID)          arcpy.Clip_management(baseFeature, "#", output,feat, "0", "ClippingGeometry") # delete the cursor del rows  ################################################################################################# ####### second part of script to take individual catchments and calculate some statistics ####### # i created this using a separate script from the first so the environments are retyped #################################################################################################  # make sure to put the module you are working with before the function (or it will spit back 'not defined') - one my biggest pitfalls ###here i want to update the catchemnt feature with additional information based on other layers import math import os import arcpy # set the workspace environment  arcpy.env.workspace = "D:/path.mdb"  # Overwrite pre-existing files arcpy.env.overwriteOutput = True  #parameters path = "D:/path.mdb" arcpy.Delete_management("SumTable") arcpy.CreateTable_management (path,"SumTable") arcpy.AddField_management ("SumTable", "NameID", "TEXT") arcpy.AddField_management ("SumTable", "TotalCount", "LONG") arcpy.AddField_management ("SumTable", "TotalElv", "FLOAT",'',10) arcpy.AddField_management ("SumTable", "AVG_ELV", "FLOAT",'',10) SummaryT= "D:/path.mdb/SumTable"   statsFields = [["COUNT", "SUM"],["ELEVC","SUM"]]   #create list to point to catchment elevations catchmentList = arcpy.ListRasters ("c*")   for i in catchmentList:     #make sure each raster has its own attribute table     arcpy.BuildRasterAttributeTable_management(i,"Overwrite")     #create a new field with elevation times number of occurances (COUNT * VALUE)     arcpy.AddField_management (i, "ELEVC", "LONG")     arcpy.CalculateField_management(i,"ELEVC","!Value! * !Count!","PYTHON_9.3","#")     temp_outtable = "outtable"     arcpy.Statistics_analysis (i, temp_outtable , statsFields)      #python's tedious way to get data from a table     rows=arcpy.SearchCursor(temp_outtable) #this 'points' to the first row      row=rows.next() #this 'points' to the field     print row.SUM_COUNT,"_",row.SUM_ELEVC  # Create the insert cursor and a new empty row     rowInserter = arcpy.InsertCursor(SummaryT)     newIncident = rowInserter.newRow()  # Populate attributes of new row     newIncident.NameID = str(i[0:7])     newIncident.TotalCount = row.SUM_COUNT     newIncident.TotalElv = row.SUM_ELEVC     newIncident.AVG_ELV = row.SUM_ELEVC / row.SUM_COUNT  # Insert the new row into the shapefile     rowInserter.insertRow(newIncident)     # delete catchment specific search cursor and temp sum file     del row     del rows     arcpy.Delete_management(temp_outtable)       # delete the insert cursor del rowInserter