Select features using geometry objects

[ScottOatley]

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

[ScottOatley]

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

[T__WayneWhitley]

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

[ScottOatley]

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

[T__WayneWhitley]

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.)

[ScottOatley]

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

Scott

[T__WayneWhitley]

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?

[ScottOatley]

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

[NeilAyres]

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

[T__WayneWhitley]

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

[ScottOatley]

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