nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura

Hi,

How can I iterate through feature layer or feature class by arcpy.da.SerachCursor and arcpy.da.UpdateCursor within a certain geographic extent? Right now I am making feature class copying the content. Is this most efficient method? Basically I have giant feature class, and I want to visit geographic subdivision of the feature class and then do my processing there.

# I have a feature class.  I made it into feature layer
arcpy.MakeFeatureLayer_management(poly_fc, tmp_lyr)

# I then set extent to the part of data that I want to process
arcpy.env.extent = extent1

# if I count the feature from the layer, I got the count of feature in the extent. 
#  if I copy it into feature class I got the featueres inside the extent
n0 = int( arcpy.GetCount_management(tmp_lyr).getOutput(0)) # I got corrent n0
arcpy.CopyFeatures_management(tmp_lyr, selected_fc)  # I got fc of only the feature within extent1

# But iterating the tmp_lyr visits all featuers in tmp_lyr, ignoring the arcpy.env.extent
lst =  [ _[0] for _ in arcpy.da.SearchCursor(tmp_lyr, ['SHAPE@'])]   
#     lst above contains all features, not just the one in extent

# if I access through the copied fc, I get what I want, i.e., features within the extent
lst_small =  [ _[0] for _ in arcpy.da.SearchCursor(selected_fc, ['SHAPE@'])]   
#     lst_small above contains what I want

# so I just process on the copy....
with arcpy.da.SearchCursor(selected_fc, ['SHAPE@']) as cur:
  for row in cur:
    # my processing here...
    pass

Find more posts tagged with

Accepted answers

bixb0012

I am using ArcGIS 10.3, and it isn't fixed!

There is a bug with the arcpy.SelectLayerByLocation_management tool. At first I thought it might be related to using shape files, but I can generate the same issues with file geodatabases and even ArcSDE databases. The following code demonstrates how arcpy.SelectLayerByLocation_management misses the polygon while looping through a search cursor and checking for intersection gets the polygon.

>>> arcpy.SelectLayerByLocation_management("xyz_feat", "INTERSECT","xyz_ext1", "", "NEW_SELECTION")
<Result 'xyz_feat'>
>>> arcpy.GetCount_management("xyz_feat")
<Result '1357'>
>>> arcpy.SelectLayerByAttribute_management("xyz_feat","CLEAR_SELECTION")
<Result 'xyz_feat'>
>>>
>>> i = 0
>>> ext = next(arcpy.da.SearchCursor("xyz_ext1","SHAPE@"))[0]
>>> with arcpy.da.SearchCursor("xyz_feat","SHAPE@") as cur:
...    for shape, in cur:
...        if not shape.disjoint(ext):
...            i += 1
...           
>>> i
1358
>>>

The code above does raise the possibility of another approach. If you are always selecting from the same set of features, you could build a search cursor and loop over the same cursor each time the extent object changes. There are a couple of pluses with this approach. One, creating a new cursor takes more resources than simply resetting one. Two, the .disjoint method is quite quick.

All comments

yosukekimura

Thank you for help!

So my extent is one giant box. my feature having trouble in particular this time is a small rectangle sitting inside the extent rectangle. I understand that INTERSECT, CONTAINS etc has slight difference but I understand that they both select the blue feature in the pic below. pink is the extent polygon

The extent is a lot larger than the feature. Below show the extent, and tiny dot near top left corner is the feature causing trouble.

bixb0012

Regarding spatial reference, yes, it is always a good idea to specify a spatial reference when creating ArcPy geometry objects. There are numerous, known precision issues that crop up with ArcPy geometry objects when spatial references aren't specified. Is lack of a spatial reference the issue here? We will know soon enough when you test it.

When you are talking about a "shared point with the extent rectangle," I am a bit uncertain whether that is what you mean. "Shared point" could be taken to mean something very specific, like a polygon in question must share a vertex with the extent rectangle. Try adding the spatial reference, and then report back. Also, are the differences you are finding only along the periphery of the extent object? Are any polygons more in the interior being handled different.

There are 15+ different types of overlaps one can specify. I picked Intersect in this case because it is the default and most common. It could be the CopyFeatures_management tool uses a different overlap type, or it could be a bug in one tool or the other.

yosukekimura

Actually, I am having trouble with Select Layer By Location approach. I am getting different results with two method. I dumped what's selected into arcmap and found that setting extent and copy method selects fearues I expected, which is any features which has shared point with the extent rectangle. Somehow method 1 failed to choose some of feature. My questions are

- Is "INTERSECT" correct method to mimic what extent does?

EDIT

Looking into further, this occurred for one small polygon near the edge of the extent (but completely contained by the extent rectangle), not being selected by Select Layer By Location method, whereas setting extent and then copying features catches this one feature. I heve everything in same sapatial reference. The extent polygon has dimension of tens of degree (WGS84), and the feature having trouble is an rectangle of half minute in each side, off from the extent boundary by about half minutes.

I am using ArcGIS 10.2.2, by the way.

# prepare extent from precalculated corner coordinates
sr = arcpy.Describe(feat_lyr).spatialReference
myextent = arcpy.Extent( ((x0,y0), (x0,y1), (x1, y1), (x1, y0)), sr )

# Method 1, using Select By Location
arcpy.env.extent = None
geom_ext = arcpy.Polygon(arcpy.Array(
    [getattr( myextent, _ ) for _ in 
    ('lowerLeft', 'upperLeft', 'upperRight', 'lowerRight') ] )
arcpy.SelectLayerByLocation(feat_lyr, 'INTERSECT', geom_ext )
geoms1 = [row[0] for row in arcpy.SearchCursor(feat_lyr, ['SHAPE@'])]

# Method 2, copy feature
arcpy.env.extent = myextent
arcpy.CopyFeatures_management(feat_lyr, feat_fc)
geoms2 = [row[0] for row in arcpy.SearchCursor(feat_fc, ['SHAPE@'])]

yosukekimura

This is what I wanted, seems like working the way I wanted. Thanks!

yosukekimura

Thanks, all my feature/layer are in_momory now.

bixb0012

Use arcpy.SelectLayerByLocation_management first to create a selection on your feature layer, then create the search cursor. If you are using ArcGIS 10.3 and already have an extent object created, you can pass the extent object's polygon property to arcpy.SelectLayerByLocation_management. If you are using ArcGIS 10.2.2 or earlier, you can manually create a polygon using the extent object's bounds.

DanPatterson_Retired

you can redirect your copy to in_memory if you want, if not, reuse the same filename during your process and delete when all done.

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.