Select to view content in your preferred language

Iterating through a feature to clip raster

3146
6
Jump to solution
12-15-2021 04:24 PM
DryCreekEng
Frequent Contributor

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.

SearchCursor error png.png 

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

 

0 Kudos
1 Solution

Accepted Solutions
DryCreekEng
Frequent Contributor

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

 

View solution in original post

0 Kudos
6 Replies
DanPatterson
MVP Esteemed Contributor
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'

polygon_extent.png


... sort of retired...
by Anonymous User
Not applicable

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.

 

 

 

0 Kudos
DryCreekEng
Frequent Contributor

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. 

0 Kudos
by Anonymous User
Not applicable

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.

0 Kudos
DanPatterson
MVP Esteemed Contributor

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


... sort of retired...
0 Kudos
DryCreekEng
Frequent Contributor

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

 

0 Kudos