topic Using Geometry Objects in Python Questions

Using Geometry Objects

JoeBorgione — Tue, 06 Mar 2018 23:19:51 GMT

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

Re: Using Geometry Objects

DanPatterson_Retired — Sat, 11 Dec 2021 08:20:12 GMT

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‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Using Geometry Objects

JoshuaBixby — Wed, 07 Mar 2018 00:20:59 GMT

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.

Re: Using Geometry Objects

XanderBakker — Sat, 11 Dec 2021 08:20:14 GMT

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]
    AddFieldByTemplate(fc_pnt, fld_res, fld_template)

    # 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


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

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

Resulting labels match:

Re: Using Geometry Objects

JoeBorgione — Wed, 07 Mar 2018 01:46:21 GMT

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!

Joshua Bixby‌

Dan Patterson‌

Xander Bakker‌

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

Re: Using Geometry Objects

DanPatterson_Retired — Wed, 07 Mar 2018 02:14:54 GMT

numpy is still faster

Re: Using Geometry Objects

JoeBorgione — Sat, 11 Dec 2021 08:20:17 GMT

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]
#
inpolys = 'quads'  
#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])
                    arcpy.SelectLayerByAttribute_management("view","ADD_TO_SELECTION",select)
                    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.

bixb0012‌

Dan_Patterson‌

xander_bakker‌