topic Re: Select features using geometry objects in Python Questions

Select features using geometry objects

ScottOatley — Fri, 07 Jun 2013 19:44:21 GMT

Greetings,

I'm writing some arcpy to select polygons that contain the centroid of my current data frame extent. Here's my code:

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
myPt = arcpy.Point((df.extent.XMax+df.extent.XMin)/2, (df.extent.YMax+df.extent.YMin)/2)
myPtGeometry = arcpy.PointGeometry(myPt)

arcpy.SelectLayerByLocation_management ("DOQQImport", "COMPLETELY_CONTAINS", myPtGeometry)

The code runs without error, and the coordinates of myPt/myPtGeometry are correct. Yet, nothing gets selected.

Am I missing some coding step, or can you not use geometry this way?

Thanks,
Scott

Re: Select features using geometry objects

danielchaboya — Sat, 11 Dec 2021 11:47:31 GMT

Are you running it with your python window on (in ArcMap) to see if it's doing anything. Write some print statements for myPt and myPtGeometry variables to see if your geometry point has coordinates.

Maybe try also including ...

DOQQImport = arcpy.mapping.ListLayers(mxd, "DOQQImport", df)

then this...

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
DOQQImport = arcpy.mapping.ListLayers(mxd, "DOQQImport", df)
myPt = arcpy.Point((df.extent.XMax+df.extent.XMin)/2, (df.extent.YMax+df.extent.YMin)/2)
myPtGeometry = arcpy.PointGeometry(myPt)
arcpy.SelectLayerByLocation_management (DOQQImport, "COMPLETELY_CONTAINS", myPtGeometry)

I think you need to specifically reference a dataset. My thought is that because "DOQQImport" is not referenced the tool fails. The python window will give you an error if this is so.

Hope this helps.

Re: Select features using geometry objects

ScottOatley — Sat, 11 Dec 2021 11:47:33 GMT

I updated my code to include printing out x and y values of myPt and myPtGeometry with an AddMessage (I'm running the script from a toolbox). Both are correct, matching what I read on the screen with my cursor at the center of my data frame.

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
myPt = arcpy.Point((df.extent.XMax+df.extent.XMin)/2, (df.extent.YMax+df.extent.YMin)/2)
arcpy.AddMessage("{0}  {1}".format(myPt.X, myPt.Y))
myPtGeometry = arcpy.PointGeometry(myPt)
arcpy.AddMessage("{0}  {1}".format(myPtGeometry.firstPoint.X, myPtGeometry.firstPoint.Y))
arcpy.SelectLayerByLocation_management ("DOQQImport", "COMPLETELY_CONTAINS", myPtGeometry)

The layer I'm trying this with is in a group. I tried moving it outside the group, but still nothing gets selected.

Re: Select features using geometry objects

danielchaboya — Sat, 11 Dec 2021 11:47:36 GMT

I updated my code to include printing out x and y values of myPt and myPtGeometry with an AddMessage (I'm running the script from a toolbox). Both are correct, matching what I read on the screen with my cursor at the center of my data frame.
import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
myPt = arcpy.Point((df.extent.XMax+df.extent.XMin)/2, (df.extent.YMax+df.extent.YMin)/2)
arcpy.AddMessage("{0}  {1}".format(myPt.X, myPt.Y))
myPtGeometry = arcpy.PointGeometry(myPt)
arcpy.AddMessage("{0}  {1}".format(myPtGeometry.firstPoint.X, myPtGeometry.firstPoint.Y))
arcpy.SelectLayerByLocation_management ("DOQQImport", "COMPLETELY_CONTAINS", myPtGeometry)
The layer I'm trying this with is in a group. I tried moving it outside the group, but still nothing gets selected.

Did you run the tool within ArcMap and with the python window open? I still think the issue has to do with your "DOQQImport" dataset. It has to first be referenced before it can be used.

Re: Select features using geometry objects

T__WayneWhitley — Fri, 07 Jun 2013 22:01:00 GMT

You may have to use arcpy.RefreshActiveView ??

Re: Select features using geometry objects

T__WayneWhitley — Sun, 12 Dec 2021 16:13:02 GMT

Daniel is right, I wasn't paying attention although you will likely also need RefreshActiveView to make the selection visible in your map's data frame...
If this is your complete code, looks like you did a good job setting the objects for map, data frame (the 1st one, or only data frame in the map?), and the geometry object. But what is "DOQQImport"? I didn't test the following, but think if you get a reference to this (as Daniel said), this may work (is DOQQImport the exact layer name in the map's TOC?) -- also, I believe in this case the overlap type INTERSECT will give you virtually the same results as COMPLETELY_CONTAINS, so I'd just go with the default:

Highlighted in red are the critical changes to your code:

import arcpy

# need a 'handle' on all 3:  'current' mxd, 1st data frame, DOQQ layer
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
DOQQImport = arcpy.mapping.ListLayers(mxd, "DOQQImport", df)

# your geometry object to select by looks fine
myPt = arcpy.Point((df.extent.XMax+df.extent.XMin)/2, (df.extent.YMax+df.extent.YMin)/2)
arcpy.AddMessage("{0}  {1}".format(myPt.X, myPt.Y))
myPtGeometry = arcpy.PointGeometry(myPt)
arcpy.AddMessage("{0}  {1}".format(myPtGeometry.firstPoint.X, myPtGeometry.firstPoint.Y))

# execute using the fetched layer and defined geom:
arcpy.SelectLayerByLocation_management (DOQQImport, "INTERSECT", myPtGeometry)
arcpy.RefreshActiveView()

arcpy.AddMessage("\nIs this your intended polygon selection result?")
arcpy.AddMessage("(If there are overlapping polys at this location, the result should be a multiple selection.)")

del mxd

Hope that helps, of course test as necessary!

Enjoy,
Wayne

Re: Select features using geometry objects

T__WayneWhitley — Mon, 10 Jun 2013 15:17:01 GMT

I'm sorry, I made an error on the ListLayers line...on all functions returning a Python list, always have to remember to specify a member of the list even if only 1 member is returned, so this is correct and think Daniel already specified the line -- I'll write it again below, also showing the use of wildcards in case you get nothing returned because the layer name in the TOC is not exact, for example:

DOQQImport = arcpy.mapping.ListLayers(mxd, "*DOQQ*", df)[0]

Again, sorry if that was misleading - and also it is desirable sometimes to loop on the list in case multiple layers are returned that you want to handle.

Enjoy,
Wayne

Re: Select features using geometry objects

ScottOatley — Sat, 11 Dec 2021 11:47:39 GMT

Hello!

I tested making a reference as suggested with no change in results - nothing gets selected.

I don't believe the "DOQQImport" needs to be referenced. The layer is in the TOC, and I've used layers in other scripts without making such a reference. As a test I tried the following. I used an existing point layer ("somePtLayer") and manually selected one of its features that is covered by my DOQQImport polygon layer. I changed the code to use the point layer and the script successfully selected the correct polygons.

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
##myPt = arcpy.Point((df.extent.XMax+df.extent.XMin)/2, (df.extent.YMax+df.extent.YMin)/2)
##arcpy.AddMessage("{0}  {1}".format(myPt.X, myPt.Y))
##myPtGeometry = arcpy.PointGeometry(myPt)
##arcpy.AddMessage("{0}  {1}".format(myPtGeometry.firstPoint.X, myPtGeometry.firstPoint.Y))
##arcpy.SelectLayerByLocation_management ("DOQQImport", "COMPLETELY_CONTAINS", myPtGeometry)
arcpy.SelectLayerByLocation_management ("DOQQImport", "COMPLETELY_CONTAINS", "somePtLayer")

I'm thinking that while geometries can be used for geoprocessing, maybe they can't be used this way to make selections. Technically, selectByLocation isn't geoprocessing.

Thanks,
Scott

Re: Select features using geometry objects

NeilAyres — Tue, 11 Jun 2013 07:48:00 GMT

When creating the PointGeometry object, include the spatial reference of the coordinate that you are creating. Normally I do this from the sr of a another similar fc using the Describe object.
Cheers and good luck,
Neil

Re: Select features using geometry objects

T__WayneWhitley — Sat, 11 Dec 2021 11:47:42 GMT

Yes, a very important distinction, Neil...

Scott, glad you solved your problem - there are however some subtleties you missed and Neil hit it right on the head. About the spatial reference, if your data frame and feature class spatial references do not match, that of course could be a reason no selection was in the result. (Another obvious reason could be no feature is otherwise available for selection, i.e., no existing feature at the select location.)

In my testing my data frame and layers were in the same 'frame' of reference --- that is to say, they shared the same spatial reference. So, as a simulation of the case where the spatial references do not match, I set the data frame to something 'off the wall' just to prove a point...to a Greenland UTM zone (I am in Florida.)

As suspected, the select by location failed to select anything because the data frame centerpoint is on a different coord system...however, the documentation for Select Layer By Location allows you to set the output coord sys environment variable to overcome this:

arcpy.env.outputCoordinateSystem

Some minor modifications of your code allow the select by location to execute successfully to make the desired selection, with the critical line in bold text below. With the output coord sys environment variable set to the spatial reference of the data frame, the select layer by location tool essentially projects on-the-fly the feature layer I use (parcels) from ArcMap's TOC (of the current map doc) to match the data frame spatial ref, which is of course the same 'reference' of the point object.

import arcpy
mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd)[0]
 
parcels = arcpy.mapping.ListLayers(mxd, '*parcel*', df)[0]
sr_dataframe = df.spatialReference
sr_featureclass = arcpy.Describe(parcels).spatialReference

arcpy.env.outputCoordinateSystem = sr_dataframe

myPt = arcpy.Point((df.extent.XMax + df.extent.XMin)/2.0, (df.extent.YMax + df.extent.YMin)/2.0)
myPtGeometry = arcpy.PointGeometry(myPt)

arcpy.AddMessage("\nThe fetched parcel layer from the TOC, {0}".format(parcels.name))
arcpy.AddMessage("\nThe data frame spatial reference name:  {0}".format(sr_dataframe.name))
arcpy.AddMessage("\nThe feature class (to select by) spatial reference name:  {0}".format(sr_featureclass.name))
arcpy.AddMessage("\nThe data frame centerpoint (myPt) - X: {0}  Y: {1}".format(myPt.X, myPt.Y))
arcpy.AddMessage("\nThe geometry (myPtGeometry) - {0}  {1}".format(myPtGeometry.firstPoint.X, myPtGeometry.firstPoint.Y))
arcpy.AddMessage("\n...now proceeding with the selection, please wait...\n")
 
arcpy.SelectLayerByLocation_management(parcels, 'INTERSECT', myPtGeometry)
arcpy.RefreshActiveView()
 
del mxd

So in 'my GIS', these are the relevant 'AddMessage' execution msgs I got running as a script tool from ArcMap with the parcel layer loaded in the current doc - notice the 'read' spatial ref names (obviously different) and that the actual point geometry obj is in the data frame reference as expected:
**********************************************
The fetched parcel layer from the TOC, PARCEL_PUBLIC

The data frame spatial reference name: Greenland_1996_UTM_Zone_23N

The feature class (to select by) spatial reference name: NAD_1983_HARN_StatePlane_Florida_East_FIPS_0901_Feet

The data frame centerpoint (myPt) - X: -3358568.40446 Y: 3292414.02514

The geometry (myPtGeometry) - -3358568.40446 3292414.02514

...now proceeding with the selection, please wait...
**********************************************

Hope that helps, an important lesson I think in the spatial ref obj technique and geoprocessing -- great catch, Neil!

Enjoy,
Wayne

...the critical documentation:
Select Layer By Location (Data Management)
Desktop Â» Geoprocessing Â» Tool reference Â» Data Management toolbox
http://resources.arcgis.com/en/help/main/10.1/index.html#//001700000072000000

Output Coordinate System (Environment setting)
Desktop Â» Geoprocessing Â» Environment settings
http://resources.arcgis.com/en/help/main/10.1/index.html#//001w00000005000000

Re: Select features using geometry objects

ScottOatley — Sat, 11 Dec 2021 11:47:44 GMT

Greetings and thanks for the help so far.

I had high hopes regarding the spatial reference issue, but I'm still not getting this to work. I've got the code updated to match Wayne's example.

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

DOQQImport = arcpy.mapping.ListLayers(mxd, 'DOQQImport', df)[0]
sr_dataframe = df.spatialReference
sr_featureclass = arcpy.Describe(DOQQImport).spatialReference

arcpy.env.outputCoordinateSystem = sr_dataframe

myPt = arcpy.Point((df.extent.XMax + df.extent.XMin)/2.0, (df.extent.YMax + df.extent.YMin)/2.0)
myPtGeometry = arcpy.PointGeometry(myPt)

##fc = 'points'
##cursor = arcpy.da.InsertCursor(fc, ["SHAPE@XY"])
##xy =(myPt.X, myPt.Y)
##cursor.insertRow([xy])

arcpy.AddMessage("\nThe fetched polygon layer from the TOC, {0}".format(DOQQImport.name))
arcpy.AddMessage("\nThe data frame spatial reference name:  {0}".format(sr_dataframe.name))
arcpy.AddMessage("\nThe feature class (to select by) spatial reference name:  {0}".format(sr_featureclass.name))
arcpy.AddMessage("\nThe data frame centerpoint (myPt) - X: {0}  Y: {1}".format(myPt.X, myPt.Y))
arcpy.AddMessage("\nThe geometry (myPtGeometry) - {0}  {1}".format(myPtGeometry.firstPoint.X, myPtGeometry.firstPoint.Y))
arcpy.AddMessage("\n...now proceeding with the selection, please wait...\n")

arcpy.SelectLayerByLocation_management(DOQQImport, 'INTERSECT', myPtGeometry)
arcpy.RefreshActiveView()

del df, mxd

It successfully runs; here's my output:

[INDENT]Executing: gotoSiteLatLong
Start Time: Tue Jun 11 14:38:19 2013
Running script gotoSiteLatLong...

The fetched polygon layer from the TOC, DOQQImport

The data frame spatial reference name: NAD_1983_UTM_Zone_15N

The feature class (to select by) spatial reference name: GCS_WGS_1984

The data frame centerpoint (myPt) - X: 466834.903914 Y: 4975689.29185

The geometry (myPtGeometry) - 466834.903914 4975689.29185

...now proceeding with the selection, please wait...

Completed script gotoSiteLatLong...
Succeeded at Tue Jun 11 14:38:24 2013 (Elapsed Time: 5.00 seconds)[/INDENT]

However, I still don't get any selected polygons. I tried changing the data frame coordinate system to match the feature class with similar results. Using the commented out code,

##fc = 'points'
##cursor = arcpy.da.InsertCursor(fc, ["SHAPE@XY"])
##xy =(myPt.X, myPt.Y)
##cursor.insertRow([xy])

I was able to add a point at the center of my current extent, so my geometry is indeed creating a point in the correct location. And there are polygons in my feature class (DOQQImport) that cover the point.

I'm stumped as to what I'm missing.

Scott

Re: Select features using geometry objects

T__WayneWhitley — Tue, 11 Jun 2013 19:09:23 GMT

Hold on a sec, is not a DOQQ a raster dataset?
...or are you simply referring to a grid outline of it?

Re: Select features using geometry objects

ScottOatley — Tue, 11 Jun 2013 19:19:26 GMT

DOQQImport is a polygon layer whose polygons represent footprints of DOQQs.

Re: Select features using geometry objects

T__WayneWhitley — Sat, 11 Dec 2021 11:47:47 GMT

Whew, well that's what I thought but thought I'd better check.

I've got a min, do you think you can export a small extent and attach it here, and I'll take a quick look at it?

EDIT: It may also be a good thing to get the number of selected features - add these 2 lines:

count = arcpy.GetCount_management(DOQQImport).getOutput(0)
arcpy.AddMessage("The number in the selection:  {0}".format(count))

Also, it could take a bit for the screen to refresh in order to show you the selection, particularly if a large number of layers have to draw at a considerable extent, etc. Sorry, have to ask: is the layer turned on, or is it perhaps obfuscated somehow? At least the get count command will tell you if a selection was successfully made. (And this is alluded to in the web help on Select Layer By Location as well but I forgot about it.)

Re: Select features using geometry objects

ScottOatley — Tue, 11 Jun 2013 19:44:00 GMT

I've attached a snippet of the DOQQImport layer. They should show up around the Minneapolis area.

Scott

Re: Select features using geometry objects

T__WayneWhitley — Tue, 11 Jun 2013 20:01:07 GMT

The code worked fine...I zoomed to the layer and because of the overlapping poly nature of the layer, I got 4 selected polygons after running the code.

Did you add the 2 lines for GetCount and to report back to you the number of features selected? Then, if a feature selection is reported, are you 'zoomed in' too far to see the selected outline?

Re: Select features using geometry objects

ScottOatley — Wed, 12 Jun 2013 02:15:59 GMT

I've been running this as a script in a toolbox. I tried again here at home starting from scratch and get the same result, no selection.

Then I opened the python window from the geoprocessing menu and pasted the code in there as is and it successfully selected some polygons!

So why does it work in the python window, but not as a script?

Thanks,
Scott

Re: Select features using geometry objects

NeilAyres — Sat, 11 Dec 2021 11:47:50 GMT

Ahm.. Gentlemen,

in the code above I still do not see you doing this -

myPtGeometry = arcpy.PointGeometry(myPt, sr)

I believe that without this, the select and any other geometry method will not work.
Yes, the object has properties (X, Y etc etc) and it looks like a point, but its not quite what it seems.

N

Re: Select features using geometry objects

T__WayneWhitley — Sat, 11 Dec 2021 11:47:53 GMT

Noted Neil, and try that if you wish. However, you are not correct on this one, particularly since the same script executes correctly from the Python window, and this is documented in the webhelp for Select Layer By Location and proven in my posted code above in setting the output extent environment, honored in my test.

...think Scott has a different problem, although I don't know at this point. Could it have something to do with running 'in process' and use of 'CURRENT' - perhaps that isn't operating properly for you Scott, and the direct cause just escapes me. Just to confirm, you are running the script tool within an ArcMap document and when you check the properties of the tool it is checked 'in process'?

Wouldn't hurt to try Neil's approach, inserting the param for the output extent into the param he noted for the geom, although likely that won't work any differently. I tried it at 10.0 SP 5 and the output extent var works there -- whatever version you are at, could be yours is exhibited a different behavior. I duplicated your setup in my 10.0 and no problem.... confounding little problem, maybe someone else can contribute that has had a similar problem.

Just curious, if you place your DOQQ layer export from yesterday that you posted, place that in a new mxd and save it at full extent of the layer (with the data frame set to the GCS coord sys, not the projected one of the layer), then run the script tool, do you get a selection?

-Wayne

...this is the full code, adding Neil's suggestion - for you to test. I thought it should work either way, setting the outputCoordinateSystem env variable or setting the spatial ref param directly on the geom object itself. It doesn't -- with the data frame set to a different spatial reference than the layer, I observed that the setting of the outputCoordinateSystem environment variable was necessary and that the spatial reference parameter setting on the geometry object was optional.

Nevertheless, both are incorporated in the code below. Also, it has the 2 lines added for GetCount/AddMessage to report back in the dialog the number of features (if any) in the selection result:

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]

DOQQImport = arcpy.mapping.ListLayers(mxd, 'DOQQImport', df)[0]
sr_dataframe = df.spatialReference
sr_featureclass = arcpy.Describe(DOQQImport).spatialReference

# Setting outputCoordinateSystem env var was necessary,
# but see myPtGeometry below as well.
arcpy.env.outputCoordinateSystem = sr_dataframe

myPt = arcpy.Point((df.extent.XMax + df.extent.XMin)/2.0, (df.extent.YMax + df.extent.YMin)/2.0)

# Neil suggested this...in testing, setting this was optional.
# The spatial reference obj is used on the geometry itself:
myPtGeometry = arcpy.PointGeometry(myPt, sr_dataframe)

arcpy.AddMessage("\nThe fetched polygon layer from the TOC, {0}".format(DOQQImport.name))
arcpy.AddMessage("\nThe data frame spatial reference name:  {0}".format(sr_dataframe.name))
arcpy.AddMessage("\nThe feature class (to select by) spatial reference name:  {0}".format(sr_featureclass.name))
arcpy.AddMessage("\nThe data frame centerpoint (myPt) - X: {0}  Y: {1}".format(myPt.X, myPt.Y))
arcpy.AddMessage("\nThe geometry (myPtGeometry) - {0}  {1}".format(myPtGeometry.firstPoint.X, myPtGeometry.firstPoint.Y))
arcpy.AddMessage("\n...now proceeding with the selection, please wait...\n")

arcpy.SelectLayerByLocation_management(DOQQImport, 'INTERSECT', myPtGeometry)
arcpy.RefreshActiveView()

count = arcpy.GetCount_management(DOQQImport).getOutput(0)
arcpy.AddMessage("The number in the selection:  {0}".format(count))

del df, mxd

Re: Select features using geometry objects

ScottOatley — Wed, 12 Jun 2013 21:30:14 GMT

I'm currently using 10.1 SP 1.

Ahm.. Gentlemen,

in the code above I still do not see you doing this -
Code:

myPtGeometry = arcpy.PointGeometry(myPt, sr)

Even though I didn't post code with pointGeometry(myPoint, sr), I did try that with a variety of coordinate system setups and with and without the env.outputCoordinateSystem setting. No combination has worked when running the script from a toolbox.

Incidentally, my intention is to run this from an add-in toolbar which I tried today. It doesn't run there either. I also tried importing my toolbox to the python window and running it by calling interactively. That also didn't run. Only pasting the code directly to the python prompt in the python window runs correctly. In all cases I'm using the identical code.

Just to confirm, you are running the script tool within an ArcMap document and when you check the properties of the tool it is checked 'in process'?

Yes.

Did you add the 2 lines for GetCount and to report back to you the number of features selected? Then, if a feature selection is reported, are you 'zoomed in' too far to see the selected outline?

I tried this, and 0 is my result.

Maybe it's a 10.1 issue. I'll see if I've got 10.0 available on some machine at work.

Scott