Select to view content in your preferred language

# Using Geometry Objects

1050
6
03-06-2018 03:19 PM
MVP Emeritus

I would like to use geometry objects (I think I do anyway) to get better performance in a workflow that checks if a point is within a given polygon.  The online help indicates:

In many geoprocessing workflows, you may need to run a specific operation using coordinate and geometry information but don't necessarily want to go through the process of creating a new (temporary) feature class, populating the feature class with cursors, using the feature class, then deleting the temporary feature class. Geometry objects can be used instead for both input and output to make geoprocessing easier.

Currently I do exactly what the online help says I can avoid, but I'm having quite a time trying to understand what it is I need to do. Can anyone give me some clues as to how I might go about:

1. Creating point geometries from a point feature class

2. Create polygon geometries from a polygon feature class

3. See which polygon object a given point geometry is within

That should just about do it....
1 Solution

Accepted Solutions
MVP Emeritus
``````in_fc = r"C:\.....\r_sorted"
flds = ['SHAPE@']
cur = arcpy.da.SearchCursor(in_fc, flds)
pnts = [i for i in cur]

type(pnts[0])
tuple

type(pnts[0][0])
arcpy.arcobjects.geometries.PointGeometry

pnts[:3]
[(<PointGeometry object at 0x2375784a828[0x23754cbb350]>,),
(<PointGeometry object at 0x2375784a7b8[0x23754cbb3a0]>,),
(<PointGeometry object at 0x2375784a8d0[0x23754cbb3f0]>,)]

other_fc = r"C:\....\r_extent"
cur2 = arcpy.da.SearchCursor(other_fc, 'SHAPE@')

cur2
<da.SearchCursor object at 0x000002375703CE68>

poly = [i for i in cur2]

poly
[(<Polygon object at 0x23757744940[0x23754cbb238]>,)]

polygeom = poly[0][0]    # --- slice the tuple
for i in pnts[:10]:
pnt = i[0].centroid  # ---- strange but true... you need the centroid of a point
print("point {}, {} is in  {}".format(pnt.X, pnt.Y, pnt.within(polygeom)))

point 303770.0, 5029863.0 is in  True
point 304571.0, 5029840.0 is in  True
point 303841.0, 5029786.0 is in  True
point 303859.0, 5029760.0 is in  True
point 304014.0, 5029729.0 is in  True
point 304504.0, 5029782.0 is in  True
point 303776.0, 5029681.0 is in  True
point 303864.0, 5029630.0 is in  True
point 304978.0, 5029822.0 is in  True
point 304460.0, 5029697.0 is in  True‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````
6 Replies
MVP Emeritus
``````in_fc = r"C:\.....\r_sorted"
flds = ['SHAPE@']
cur = arcpy.da.SearchCursor(in_fc, flds)
pnts = [i for i in cur]

type(pnts[0])
tuple

type(pnts[0][0])
arcpy.arcobjects.geometries.PointGeometry

pnts[:3]
[(<PointGeometry object at 0x2375784a828[0x23754cbb350]>,),
(<PointGeometry object at 0x2375784a7b8[0x23754cbb3a0]>,),
(<PointGeometry object at 0x2375784a8d0[0x23754cbb3f0]>,)]

other_fc = r"C:\....\r_extent"
cur2 = arcpy.da.SearchCursor(other_fc, 'SHAPE@')

cur2
<da.SearchCursor object at 0x000002375703CE68>

poly = [i for i in cur2]

poly
[(<Polygon object at 0x23757744940[0x23754cbb238]>,)]

polygeom = poly[0][0]    # --- slice the tuple
for i in pnts[:10]:
pnt = i[0].centroid  # ---- strange but true... you need the centroid of a point
print("point {}, {} is in  {}".format(pnt.X, pnt.Y, pnt.within(polygeom)))

point 303770.0, 5029863.0 is in  True
point 304571.0, 5029840.0 is in  True
point 303841.0, 5029786.0 is in  True
point 303859.0, 5029760.0 is in  True
point 304014.0, 5029729.0 is in  True
point 304504.0, 5029782.0 is in  True
point 303776.0, 5029681.0 is in  True
point 303864.0, 5029630.0 is in  True
point 304978.0, 5029822.0 is in  True
point 304460.0, 5029697.0 is in  True‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````
MVP Esteemed Contributor

Can you elaborate a bit more on your workflow?  For example, do you have a feature class of points and want to see which polygons from another feature class they intersect?

Seeing your data is already in feature classes, there is a good chance geoprocessing tools will outperform converting all the geometries to geometry objects and processing them that way.

Esri Esteemed Contributor

Here an example of what you can do. See code below:

• just define the point featureclass and the output field to store the result from the "spatial join" (line 5 and 6)
• define the polygon featureclass and field you want to take over to the point featureclass (lines 7 and 😎

``````def main():
import arcpy

# settings
fc_pnt = r'C:\GeoNet\SpatialJoinPy\test.gdb\points'
fld_res = 'result'
fc_pol = r'C:\GeoNet\SpatialJoinPy\test.gdb\polygons'
fld_pol = 'SCaNombre'

# create a dct for the polygons OID:(SHAPE, FieldOfInterest)
flds = ('OID@', 'SHAPE@', fld_pol)
dct_pol = {r[0]: (r[1], r[2]) for r in arcpy.da.SearchCursor(fc_pol, flds)}

# Add output result field if it does not exist
fld_template = arcpy.ListFields(fc_pol, fld_pol)[0]

# update cursor for update the points
flds = ('SHAPE@', fld_res)
with arcpy.da.UpdateCursor(fc_pnt, flds) as curs:
for row in curs:
pntg = row[0]
result = GetValueFromPolygonThatContainsPoint(pntg, dct_pol)
if not result is None:
row[1] = result
curs.updateRow(row)

def GetValueFromPolygonThatContainsPoint(pntg, dct):
result = None
for oid, data in dct.items():
polygon = data[0]
if polygon.contains(pntg): # or if pntg.within(polygon):
result = data[1]
break
return result

if len(arcpy.ListFields(fc, fld_name)) == 0:
fld_template.precision, fld_template.scale,
fld_template.length)

if __name__ == '__main__':
main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````

Resulting labels match:

MVP Emeritus

Joshua- Yes; I have a point feature class and Polygon featureclass.  That quoted paragraph from the help could have said "Don't do what Joe does....  He creates feature layers and performs spatial selections as he loops through the points"  It's painfully slow and as it loops the screen flickers and drives (me and) the users nuts.

Ultimately what happes is I compare the values of a couple of point attributes to the value of an attribute in the [selected] polygon; if the values are equal that's a good thing, if not, it's a bad thing.  So there is a spatial component and then an attribute component.

I just found myself going in circles through the help pages starting with Geometry, then Point Geometry, etc etc. I'll check out what Dan and Xander provide herein when I get back to to office in the morning.

Thanks to all!

Edited moments later...  Just took a look at a couple geoprocessing tools, like intersect.  That could prove frutiful....

That should just about do it....
MVP Emeritus

numpy is still faster

MVP Emeritus

I came up with a workable solution using geometry objects: Joshua, I did try a spatial join to a feature class in memory and that worked too.  But really wanted to avoid creating even an in_memory feature class.

Using Dan's approach as a guide, I was able to get not only the spatial component but the attribute to attribute comparison.  Here goes:

``````flds = ['OID@', 'SHAPE@','HouseDir','StreetDir']
# housedir & streetdir are populated with N,S,E,W
inpoints = 'points'
# 'points' is a feature layer of points already added to mxd
#
cur = arcpy.da.SearchCursor(inpoints, flds)
pnts = [i for i in cur]
#
#quads is a feature layer of a polygon feature class
# of the four quadrants of the Salt Lake valley

flds = flds = ['OID@', 'SHAPE@','DirValues']
#DirValues describes the quadrant as ('N','E'), ('S','E') etc

cur2 = arcpy.da.SearchCursor(inpolys,flds)
polys = [i for i in cur2]

# nest a couple of for loops;
# what quad is a point within
# and do the HouseDir and StreetDir
# match the DirValues

for i in pnts:
pnt_centroid = i[1].centroid
for p in polys:
poly_geometry = p[1]
if pnt_centroid.within(poly_geometry):
hd = i[1]
sd = i[2]
polydirs = p[2]
if hd is None:
pass
elif hd == ' ':
pass
elif sd is None:
pass
elif sd == ' ':
pass
else:
if hd and sd in polydirs:
pass
else:
# here is where the real work is done
# 'view' is table view of the the 'points' feature layer
#  select based on the OID of the point/veiw record
#  and write it to an Errors table with with a
#  description if what's wrong in the 'expression' variable

select = 'OBJECTID = {}'.format(i[0])
if arcpy.Describe("view").FIDSet =='':
pass
elif int(arcpy.Describe("view").FIDSet[0]) >= 1:
arcpy.Append_management("view",temptable,"NO_TEST")
# temptable is a table in_memory that acts as a temporary
# holding bin
expression = "HouseDir or StreetDir in Wrong Quad"
arcpy.CalculateField_management(temptable,"ErrorType","'{}'".format(expression),"PYTHON")
outtable = r"C:\Replicas\Replicas.gdb\Errors"
arcpy.Append_management(temptable,outtable,"TEST")
arcpy.DeleteRows_management("tableinmemory")
arcpy.SelectLayerByAttribute_management("view","CLEAR_SELECTION")‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍``````

Take away from today's lesson: Geometry objects can be a bugger to figure out, but once that's accomplished, they are a powerful tool to use in your analyses.

That should just about do it....