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

[NeilAyres]

There may be something wrong with the geometry of this one feature. Run Check Geometry or Repair Geometry on your polygons.

[yosukekimura]

Thank you for advice.

I ran CheckGeometry on both extent polygon and features polygon, neither returned anything.  I ran RepearGeometry tool on both, but the symptom was the same.

[JoshuaBixby]

Thanks for uploading data, it helps.

If I use your data as-is, and choose Intersect_3D instead of Intersect, the polygon in question gets selected for a total of 1358.  If I project your data from WGS84 to a North American continental projection like Lamberts Conformal Conic, I get the same result.

Since using the 3D intersect option and projecting your data explicitly gets the same, expected result; I can only surmise something is going on with default values or tolerances when using geographic coordinates and a 2D intersect.  Maybe the projection-on-the-fly is having an issue, I am not sure.

[yosukekimura]

Joshua,

Thank you for taking time testing this out.  Interesting that INTERSECT_3D behave differently, but I don't know how to interpret this behavior.  I never used 3D feature in ArcGIS, and this particular feature class was created by me starting from arcpy.Polygon() based on input coordinate values I got.

The fact that making feature class of only this single problematic polygon let the ArcGIS to behave as expected make me speculate that somewhere ArcGIS is setting analysis resolution of coordinates and the extent of original feature class is so large compared distance of the polygon to extent.  But the SelectLayerByLocation doenst observe much environmental settings, maybe too dangerous for user to arbitrary set some threshold for vector data.  Hopefully in version 10.3, passing Extent object will "fix" the problem I see.

For the time being I will stick with CopyFeatures method, because it is working and robustness is more important since this is part of bigger script.    Part of input for my script is your employer, by the say, US Forest Service

[JoshuaBixby · Accepted solution]

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("xxxx_feat", "INTERSECT","xxxx_ext1", "", "NEW_SELECTION")
<Result 'xxxx_feat'>
>>> arcpy.GetCount_management("xxxx_feat")
<Result '1357'>
>>> arcpy.SelectLayerByAttribute_management("xxxx_feat","CLEAR_SELECTION")
<Result 'xxxx_feat'>
>>>
>>> i = 0
>>> ext = next(arcpy.da.SearchCursor("xxxx_ext1","SHAPE@"))[0]
>>> with arcpy.da.SearchCursor("xxxx_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.

[yosukekimura]

This sounds a good work around.  It doesn't resolve the issue of Select Layer By Location, but what I wanted to do from the first place was to cycle through feature within extent.  I really don't have to select anything.  I will try this in my code see if performance is good.   Thank you so much for all the help!

[JoshuaBixby]

Turns out, existing bug:

Esri Support:

BUG-000083370 : Select by Location does not select all features that meet the criteria. Workround - Import the extent of the Source layer into the Target layer

To workaround this issue, you have a couple of options. The first is to create a new feature above the square that is not being selected out side of the source layer. This will force the extent of the feature class out and allow the square feature to be selected. Once the extent has been extended you can go ahead and delete the other polygon. The other option is to temporarily load the large square (xxxx_ext1) into the xxxx_feature feature class to pull out the extent. This will guarantee that the extent is correct in every direction.