Select features using geometry objects

9823
21
06-07-2013 12:44 PM
ScottOatley
New Contributor II
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
Tags (2)
0 Kudos
21 Replies
danielchaboya
New Contributor III
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.
0 Kudos
ScottOatley
New Contributor II
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.
0 Kudos
danielchaboya
New Contributor III
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.
0 Kudos
T__WayneWhitley
Frequent Contributor
You may have to use arcpy.RefreshActiveView ??
0 Kudos
T__WayneWhitley
Frequent Contributor
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
0 Kudos
T__WayneWhitley
Frequent Contributor
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
0 Kudos
ScottOatley
New Contributor II
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
0 Kudos
NeilAyres
MVP Alum
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
0 Kudos
T__WayneWhitley
Frequent Contributor
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
0 Kudos