arcpy problem grabbing SHAPE@ token using index

258
9
Jump to solution
01-20-2022 06:39 AM
Labels (1)
lkline97
New Contributor II

Hi all,

I have written code that selects rows of a buffer shapefile based on that file's intersection with a polygon that represents the extent of a raster. In the photo below, the code selects all of the red dots for this given polygon.

image.png

Using a cursor, I am going through those selected buffers and clipping an associated raster by each buffer. The code below was working when I wrote it three months ago; however, since updating, I am having a problem grabbing the SHAPE@ token as a polygon to clip the raster. EDIT: I would have been using Version 2.8.3 when I wrote and ran the code successfully, and I have Version 2.9.1 now.

I believe I've isolated the problem to the row[1] command. Previously, I'm under the impression that I was grabbing a polygon object. Now, however, my output for the lines,

 

 

print(' row: ', row)
print(' row[1]: ', row[1])

 

 

is showing row[1] as two different outputs, a Polygon object that I want and the geoprocessing geometry object that it now gives me:

 

 

row: ['LCRas_0', <Polygon object at 0x1f56ab62108[0x1f56db53300]>]
row[1]: <geoprocessing describe geometry object object at 0x000001F56DB53300>

 

 

The clip code that once worked is now giving me a blank output. I've also tried Extract by Mask, but receive an error (good ole 99999) when I try to specify that the output extent should be that of the buffer.

Could this be a bug in the versions, or am I missing something simple? Thank you for any input!

 

 

env.workspace = r'I:\test\Polygons'
arcpy.env.addOutputsToMap = False
fields = ['RasterName', 'SHAPE@']
j=0 #code was made pieacemeal to work on 3 comps; number changes depending where in the shapefile you're at
buff_file = r"I:\test\test_buffs.shp"

for poly in listp:
#record lc file that corresponds
     lc_file = poly[:-9].replace(" ", "")
     path = r"I:\test\Classified"
     cor_rast = os.path.join(path, str(lc_file+"_classified.tif"))
     raster_xmin = arcpy.GetRasterProperties_management(cor_rast, "LEFT")
     raster_xmax = arcpy.GetRasterProperties_management(cor_rast, "RIGHT")
     raster_ymin = arcpy.GetRasterProperties_management(cor_rast, "BOTTOM")
     raster_ymax = arcpy.GetRasterProperties_management(cor_rast, "TOP")
     extent = raster_xmin[0] + " " + raster_ymin[0] + " " + raster_xmax[0] + " " + raster_ymax[0]
     print('Starting arcpy.da.UpdateCursor')

     with arcpy.da.UpdateCursor(arcpy.management.SelectLayerByLocation(buff_file, "INTERSECT", poly),fields) as cursor: #select by polygon
          for row in cursor:
               rname = str("LCRas_" + str(j)) #make raster name
               path = r"I:\test\Rasters" #set path where raster output will happen
               rpath = os.path.join(path, str(rname+".tif")) #create name of raster file
               print(' Starting arcpy.Clip_management for ', rpath)
               print(' cor_rast: ', cor_rast)
               print(' raster_xmin: ', raster_xmin)
               print(' raster_xmax: ', raster_xmax)
               print(' raster_ymin: ', raster_ymin)
               print(' raster_ymax: ', raster_ymax)
               print(' extent: ', extent)
               print(' row: ', row)
               print(' row[1]: ', row[1])
               arcpy.Clip_management(cor_rast, extent, rpath, row[1], "0", "ClippingGeometry", 'MAINTAIN_EXTENT') #extract raster from land cover
               print(' Ended arcpy.Clip_management for ', rpath)
               row[0] = rname #set row to equal raster name
               cursor.updateRow(row) #updates row to have raster column
               j=j+1
print('Ended arcpy.da.UpdateCursor')

 

 

 

P.S. I am cross-posting this on GIS stack exchange and will share any updates if I get them there as well.

0 Kudos
1 Solution

Accepted Solutions
lkline97
New Contributor II

Hi everyone,

It looks like I may have sent us on a wild goose chase. There seems to have been an issue with the extent rather than the polygon. By chance I realized that the output from the code above was not a blank raster, but instead one at such a large scale that I couldn't see the small 0.6 radius raster. 

When I ran the Clip Raster tool in ArcGIS Pro, I got the following Python output for clipping the raster by the first buffer successfully:

arcpy.management.Clip(cor_rast, "464910.1475 4860275.062 535365.525099999 4890685.978", path, "test_buffs", "15", "ClippingGeometry", "NO_MAINTAIN_EXTENT")

 This produces a raster with the appropriate final extent of "465252.920 465254.140 4860631.900 4860633.120".

When I used this code:

for poly in listp:
#record lc file that corresponds
     lc_file = poly[:-9].replace(" ", "")
     path = r"I:\test\Classified"
     cor_rast = os.path.join(path, str(lc_file+"_classified.tif"))
     raster_xmin = arcpy.GetRasterProperties_management(cor_rast, "LEFT")
     raster_xmax = arcpy.GetRasterProperties_management(cor_rast, "RIGHT")
     raster_ymin = arcpy.GetRasterProperties_management(cor_rast, "BOTTOM")
     raster_ymax = arcpy.GetRasterProperties_management(cor_rast, "TOP")
     extent = raster_xmin[0] + " " + raster_ymin[0] + " " + raster_xmax[0] + " " + raster_ymax[0]

to get the extent for the  classified raster, I instead got that the extent of the classified raster is

464479.22 4859161.04 465380.34 4860963.28

Instead of setting extent in line 33 as the raster extent, I changed it to be the polygon extent. This was successful. The fixed code is below. I have no idea why the code worked when I wrote it three months ago, but it seems to be in order now! Thank you to everyone for the help; it made me sit and consider my code more critically.

env.workspace = r'I:\test\Polygons'
arcpy.env.addOutputsToMap = False
fields = ['RasterName', 'SHAPE@']
j=0 #code was made pieacemeal to work on 3 comps; number changes depending where in the shapefile you're at
buff_file = r"I:\test\test_buffs.shp"

for poly in listp:
    #record lc file that corresponds
    lc_file = poly[:-9].replace(" ", "")
    path = r"I:\test\Classified"
    cor_rast = os.path.join(path, str(lc_file+"_classified.tif"))
    print('Starting arcpy.da.UpdateCursor')   
    
    with arcpy.da.UpdateCursor(arcpy.management.SelectLayerByLocation(buff_file, "INTERSECT", poly),fields) as cursor: #select by polygon
        for row in cursor:
            rname = str("LCRas_" + str(j)) #make raster name
            path = r"I:\test\Rasters" #set path where raster output will happen
            rpath = os.path.join(path, str(rname+".tif")) #create name of raster file
            extent = row[1].extent
            extentBuff = str(extent.XMin) + " " + str(extent.YMin) + " " + str(extent.XMax) + " " + str(extent.YMax)            
            arcpy.management.Clip(cor_rast, extentBuff, rpath, row[1], "15", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
            print('   Ended arcpy.Clip_management for ', rpath)
            row[0] = rname #set row to equal raster name
            cursor.updateRow(row) #updates row to have raster column 
            j=j+1  
    print('Ended arcpy.da.UpdateCursor')

View solution in original post

9 Replies
JoshuaBixby
MVP Esteemed Contributor

I haven't had a chance to dive in, but you should always share what version it last worked for you and what version you are using now.  Saying "since updating, I am having a problem" doesn't give folks much to go on, especially if you think the issue is related to a change in Pro versions.

lkline97
New Contributor II

Thank you for the feedback! I will edit this into the post, but I would have been using Version 2.8.3 when I wrote and ran the code successfully, and I have Version 2.9.1 now.

0 Kudos
DanPatterson
MVP Esteemed Contributor

Your indentation doesn't appear to be correct... perhaps a copy paste issue

Code formatting ... the Community Version - Esri Community


... sort of retired...
0 Kudos
lkline97
New Contributor II

Copy/paste error indeed. I've editted the original post so the indentation should be correct now. Thanks!

0 Kudos
JoshuaBixby
MVP Esteemed Contributor

I just applied the 2.9.1 patch, and I ran some basic tests against polygon feature classes, and I am not seeing what you are seeing.  I get <Polygon object at ...> when printing row and row[1].

Can you share a sample of the data set?

0 Kudos
lkline97
New Contributor II

Hmm, interesting. I hope it's okay to share through a Google Drive link as the files are too large to share here. You can find those here. line 7 references a variable listp that I didn't specify above; that should be set to the shapefile in the Polygons folder. Please let me know if I can clarify anything, and I appreciate the help!

0 Kudos
and_viceversa
Occasional Contributor

It looks like that's working as intended.

with arcpy.da.UpdateCursor(r'path/to/some_polygon.shp', 'SHAPE@') as cursor:
    for row in cursor:
        print(row)
        print(type(row))
        print(row[0])
        print(type(row[0]))

[<Polygon object at 0x286a858f748[0x286a80f2a20]>]
<class 'list'>
<geoprocessing describe geometry object object at 0x00000286A80F2A20>
<class 'arcpy.arcobjects.geometries.Polygon'>

 

I think you're still getting the Polygon? The output you showed is "printing" that object, which does not necessarily show the type, but shows the string conversion such as a Describe object.

Perhaps the problem is elsewhere? 

.
0 Kudos
DanPatterson
MVP Esteemed Contributor

Once you get a polygon there is loads you can get

p = array_poly(sq2, "POLYGON")
p
[(<Polygon object at 0x190d54416c8[0x190d53f25a0]>,), ... snip ...,
 (<Polygon object at 0x190d5441888[0x190d5455510]>,)]

p0 = p[0][0]
dir(p0)
['JSON', 'WKB', 'WKT', '__add__', '__class__', '__cmp__', '__delattr__',
 '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__from_scripting_arc_object__', '__ge__', '__geo_interface__',
 '__getSVG__', '__getattribute__', '__getitem__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__iter__', '__le__', '__len__', '__lt__',
 '__module__', '__ne__', '__new__', '__or__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__str__', '__sub__',
 '__subclasshook__', '__type_mapping__', '__type_string__', '__weakref__', '__xor__', '_arc_object', '_fromGeoJson', '_go', '_passthrough',
 '_repr_svg_', 'angleAndDistanceTo', 'area', 'boundary', 'buffer', 
 'centroid', 'clip', 'contains', 'convexHull', 'crosses', 'cut',
 'densify', 'difference', 'disjoint', 'distanceTo', 'equals', 'extent', 
 'firstPoint', 'generalize', 'getArea', 'getGeohash', 'getLength', 'getPart',
 'hasCurves', 'hullRectangle', 'intersect', 'isMultipart', 'labelPoint', 
 'lastPoint', 'length', 'length3D', 'measureOnLine', 'overlaps', 'partCount',
 'pointCount', 'pointFromAngleAndDistance', 'positionAlongLine', 'projectAs',
 'queryPointAndDistance', 'segmentAlongLine', 'snapToLine',
 'spatialReference', 'symmetricDifference', 'touches', 'trueCentroid', 'type', 'union', 'within']

until such time, you have just the polygon

polygon.png

But there is lots of other stuff you can explore on your own from the above

p0.isMultipart
False

p0.hullRectangle
'2.94135809733738 11.7648938434666 11.1766832588924 9.70596830869874 8.23532516155499 -2.05892553476782 0 0'

Or if you want to convert shapes or get things like the convex hull etc

ch.png

explore, there is loads of stuff that can be obtained directly from the object itself


... sort of retired...
Tags (3)
lkline97
New Contributor II

Hi everyone,

It looks like I may have sent us on a wild goose chase. There seems to have been an issue with the extent rather than the polygon. By chance I realized that the output from the code above was not a blank raster, but instead one at such a large scale that I couldn't see the small 0.6 radius raster. 

When I ran the Clip Raster tool in ArcGIS Pro, I got the following Python output for clipping the raster by the first buffer successfully:

arcpy.management.Clip(cor_rast, "464910.1475 4860275.062 535365.525099999 4890685.978", path, "test_buffs", "15", "ClippingGeometry", "NO_MAINTAIN_EXTENT")

 This produces a raster with the appropriate final extent of "465252.920 465254.140 4860631.900 4860633.120".

When I used this code:

for poly in listp:
#record lc file that corresponds
     lc_file = poly[:-9].replace(" ", "")
     path = r"I:\test\Classified"
     cor_rast = os.path.join(path, str(lc_file+"_classified.tif"))
     raster_xmin = arcpy.GetRasterProperties_management(cor_rast, "LEFT")
     raster_xmax = arcpy.GetRasterProperties_management(cor_rast, "RIGHT")
     raster_ymin = arcpy.GetRasterProperties_management(cor_rast, "BOTTOM")
     raster_ymax = arcpy.GetRasterProperties_management(cor_rast, "TOP")
     extent = raster_xmin[0] + " " + raster_ymin[0] + " " + raster_xmax[0] + " " + raster_ymax[0]

to get the extent for the  classified raster, I instead got that the extent of the classified raster is

464479.22 4859161.04 465380.34 4860963.28

Instead of setting extent in line 33 as the raster extent, I changed it to be the polygon extent. This was successful. The fixed code is below. I have no idea why the code worked when I wrote it three months ago, but it seems to be in order now! Thank you to everyone for the help; it made me sit and consider my code more critically.

env.workspace = r'I:\test\Polygons'
arcpy.env.addOutputsToMap = False
fields = ['RasterName', 'SHAPE@']
j=0 #code was made pieacemeal to work on 3 comps; number changes depending where in the shapefile you're at
buff_file = r"I:\test\test_buffs.shp"

for poly in listp:
    #record lc file that corresponds
    lc_file = poly[:-9].replace(" ", "")
    path = r"I:\test\Classified"
    cor_rast = os.path.join(path, str(lc_file+"_classified.tif"))
    print('Starting arcpy.da.UpdateCursor')   
    
    with arcpy.da.UpdateCursor(arcpy.management.SelectLayerByLocation(buff_file, "INTERSECT", poly),fields) as cursor: #select by polygon
        for row in cursor:
            rname = str("LCRas_" + str(j)) #make raster name
            path = r"I:\test\Rasters" #set path where raster output will happen
            rpath = os.path.join(path, str(rname+".tif")) #create name of raster file
            extent = row[1].extent
            extentBuff = str(extent.XMin) + " " + str(extent.YMin) + " " + str(extent.XMax) + " " + str(extent.YMax)            
            arcpy.management.Clip(cor_rast, extentBuff, rpath, row[1], "15", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
            print('   Ended arcpy.Clip_management for ', rpath)
            row[0] = rname #set row to equal raster name
            cursor.updateRow(row) #updates row to have raster column 
            j=j+1  
    print('Ended arcpy.da.UpdateCursor')