Hello ESRI community,
I have a large raster that covers many farm fields. I would like to iterate through a feature class and clip the raster by each farm field polygon in the feature class. Originally I tried to do this in the Python window of ArcPro, but couldnt get ListRasters to work with the listMaps class. So, I am on to a stand-alone python script.
I keep getting an error when I run the code below. Clip_management wants four points to define the clipping rectangle but I specify a polygon geometry ('SHAPE@'). Im assuming that is the problem, but Im not sure how to fix it.
arcpy.env.workspace = r'C:\GIS\Data\10000_VGF\VGF_Satellite_Imagery\VGF_Satellite_Imagery.gdb'
out_loc = r'C:\GIS\Data\10000_VGF\VGF_Satellite_Imagery\VGF_Imagery_NDVI.gdb'
raster = arcpy.ListRasters('Imagery_merged_2021')[0]
farm = r'C:\GIS\Data\10000_VGF\VGF_Field_Boundaries\VGF_Field_Boundaries.gdb\VGF_Field_Boundaries'
fields = ['SHAPE@', 'Field_Name2', 'Farm']
where_clause = "Farm = 'VGF'"
i = 0
with arcpy.da.SearchCursor(farm, fields, where_clause) as cursor:
for row in cursor:
i = +1
arcpy.Clip_management(raster, row[1], out_loc + '\\' + row[2] + "_"+ str(i) + "_2021_clip", in_template_dataset = row[0], clipping_geometry = "ClippingGeometry") #clip based on layer, clipping geometry will use the polygon extent only
Solved! Go to Solution.
I got it to work with both of your help and much fiddling. Here are my notes of what I changed:
1. Changed the names of some farm fields because they started with numbers, and gdb's hate files that start with a number.
2. Reorganized the fields list to have the field names first, then geometry (SHAPE@), as seen in the SearchCursor help page. I have learned in python that in many instances, the particular order in which things are listed really does matter.
import arcpy
arcpy.env.workspace = r"C:\GIS\VGF_Satellite_Imagery.gdb"
out_loc = r"C:\GIS\VGF_Imagery_NDVI.gdb"
raster = arcpy.ListRasters('Imagery_merged_2021')[0]
farm = r"C:\GIS\VGF_Field_Boundaries.gdb\VGF_Field_Boundaries"
fields = ['Field_Name2', 'Farm', 'SHAPE@']
where_clause = "Farm = 'VGF'"
with arcpy.da.SearchCursor(farm, fields, where_clause) as cursor:
for row in cursor:
arcpy.Clip_management(raster, "", out_loc + '\\' + row[0] + "_2021", row
[2], "0", "ClippingGeometry", "MAINTAIN_EXTENT")
import arcpy
shps = []
with arcpy.da.SearchCursor(fc0, "SHAPE@") as cursor:
for row in cursor:
shp = row[0]
shps.append(shp)
s0 = shps[0]
ext = str(s0.extent)
" ".join([i for i in ext.split(" ")[:4]])
'300000 5000000 300010 5000010'
You are setting the rectangle parameter of the function to what is in your Fieldname2 by assigning row[1] in that spot. Is that field a X-min, Y-min, X-Max, Y-Max rectangle?
I remember something similar happening a few months ago with a student project and I think the solution we found was to '' or None the rectangle parameter if the ClippingGeometry is set. I think we discussed it a zoom call and I can't find any emails talking about it to verify that is what we did so it may or may not be a solution. I remember submitting documentation feedback for it though, but I don't see any mention of it yet.
Thanks Jeff. I tried your suggestion of using "" for the rectangle parameter since I set the clipping geometry. It worked for the first field, but did not update to use the next field boundary row as the clipping geometry for the next iteration. The program ran, but the output was 37 clipped rasters that were exactly the same.
I wonder if you have to select the polygon feature with a selectByAttribute on the featureclass, instead of passing the geometry. So maybe get a list of objectids from the search cursor and sql query and then iterate over them selecting the objectid, using the FC with selected polygon to clip, repeating.
I know you can pass the geometry to some functions, but I cant remember if clip management is one of them.
the code example shows how to get the extent of input geometry through a cursor. Rasters are clipped by rectangles in any event, If you want the outside of the polygon shape (not its extent) to be nodata, then extra steps are needed
I got it to work with both of your help and much fiddling. Here are my notes of what I changed:
1. Changed the names of some farm fields because they started with numbers, and gdb's hate files that start with a number.
2. Reorganized the fields list to have the field names first, then geometry (SHAPE@), as seen in the SearchCursor help page. I have learned in python that in many instances, the particular order in which things are listed really does matter.
import arcpy
arcpy.env.workspace = r"C:\GIS\VGF_Satellite_Imagery.gdb"
out_loc = r"C:\GIS\VGF_Imagery_NDVI.gdb"
raster = arcpy.ListRasters('Imagery_merged_2021')[0]
farm = r"C:\GIS\VGF_Field_Boundaries.gdb\VGF_Field_Boundaries"
fields = ['Field_Name2', 'Farm', 'SHAPE@']
where_clause = "Farm = 'VGF'"
with arcpy.da.SearchCursor(farm, fields, where_clause) as cursor:
for row in cursor:
arcpy.Clip_management(raster, "", out_loc + '\\' + row[0] + "_2021", row
[2], "0", "ClippingGeometry", "MAINTAIN_EXTENT")