topic Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor in Python Questions

Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Sat, 11 Dec 2021 17:47:17 GMT

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

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

DanPatterson_Retired — Mon, 09 Mar 2015 04:31:27 GMT

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.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

JoshuaBixby — Mon, 09 Mar 2015 14:36:59 GMT

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.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Mon, 09 Mar 2015 15:39:17 GMT

Thanks, all my feature/layer are in_momory now.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Mon, 09 Mar 2015 16:58:03 GMT

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

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Sat, 11 Dec 2021 17:47:20 GMT

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@'])]

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

JoshuaBixby — Mon, 09 Mar 2015 23:28:19 GMT

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.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Mon, 09 Mar 2015 23:51:03 GMT

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.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

JoshuaBixby — Tue, 10 Mar 2015 01:14:42 GMT

If I understand you correctly, it is one polygon on or near the periphery that is not being handled the way you think. A few questions.

Is the missing polygon a single-part or multi-part polygon?
You appear to be working at continental map scale, i.e., 1:60,000,000. What spatial reference are you using?
Which part specifically isn't working? Does the arcpy.SelectLayerByLocation_management tool miss the polygon or does that step work and it is the search cursor that is missing it?
Regarding cursors, you should really consider using the newer ArcPy Data Access (arcpy.da) cursors instead of the older/original ArcPy (arcpy) cursors.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Tue, 10 Mar 2015 03:30:05 GMT

Thank you so much again for the help!

1. It is a single part. Most of features are, and this particular one definitely is.

2. I am using WGS84, and everything is made to WGS84. I figured how to set reference to arcpy.Polygon, so the extent polygon knows the spatial reference now.

3. I actually tested what I am doing in my ArcMap. Specifically, I run the script but I also exported the extent polygon, the feature class, feature layers, while I set the extent to be global first. The image I had was showing those exported feature (shape file). I then manually tried to select by geometry in two way. One is from the ArcMap's tool (Selection ->Select By Location ...). I also run "Select Layer By Location" tool from the tool box. Neither of them did not select this blue feature, even though blue is included in the pink (extent). Extent's top is 56.565833 dd, and the feature's top is 56.557500 dd , so off by 0.0833 degree (~ 0.5 minutes).

I also have a new observation. I made the problematic feature into a feature class lone feature. Then this select by location feature works, it picks up the feature. But if I apply the tool to the original feature which has ~1600 features spreading across the continent, that particular feature is not selected. The features are wide spread across the area, and each polygon is tiny in comparison to the extent. So I am wondering ArcGIS somehow sets criteria to lower for the spatial accuracy and then this polygon got dropped somehow...?

I made three shape files, and posted here

https://utexas.box.com/s/a6q45925jrpxzknxdp357fremtjhprhz

It has three features, ext1, feat1 and feat1_sel. ext1 is the extent polygon, feat1 the feature class of ~1600 features. feat1_sel is the one having trouble.

4. sorry for the type for cursor. I always use arcpy.da.cursor. They should have just overwritten old cursor.

yosuke

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

NeilAyres — Tue, 10 Mar 2015 06:52:50 GMT

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

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

JoshuaBixby — Tue, 10 Mar 2015 18:14:09 GMT

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.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Tue, 10 Mar 2015 19:55:55 GMT

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.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Tue, 10 Mar 2015 20:16:20 GMT

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

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

JoshuaBixby — Sat, 11 Dec 2021 17:47:23 GMT

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.

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

yosukekimura — Wed, 11 Mar 2015 20:34:23 GMT

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!

Re: Set spatial extent to arcpy.da.SerachCursor and UpdateCursor

JoshuaBixby — Mon, 23 Mar 2015 18:07:13 GMT

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.