nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Using centroid geometry to get values of intersecting feature

ChrisMathers

I know I can get the centroid for a row from the geometry object and I know in theory that I can use that to intersect with another FC and grab some of its values. I just dont know how

Find more posts tagged with

Accepted answers

All comments

ChrisMathers

Using the polygon.centroid gives you the true centroid, your red x, if it is inside the polygon. It if falls outside the polygon you are given the label point which is always inside by default. You can use the polygon.trueCentroid method to get the actual centroid regardless of its placement.

DonovanCameron

I've noticed when determining centroids for polygon features, they are calculated based on the bounding feature of the polygon (for some features, but not others...)

Either envelope, or convex hull.

This is problematic, because for a polygon that resembles the shape of the moon (in crescent), the centroid would be placed outside the polygons actual surface area.

[ATTACH=CONFIG]3845[/ATTACH]
(red crosshair is where the centroid is placed, green is where it would be preferred)

Of course, if using the centroid geometry of a point as reference, then this discrepency can be ignored.

&#109.zip

ChrisMathers

You can save arcpy.Geometry objects to a dictionary. If you create the object instead of just reading it out from the cursor. Using an arcpy.array (I havnt tried using a normal python array) to store points to use with arcpy.Polygon will make an object you can save to and use from a dictionary. Populating a dictionary with the geometry objects taht way works and doesnt take very long even for a large FC. I used it successfully to compare the centroids of one FC to the polygons of another. I just used the example from the Reading Geometries page to populate the arrays.

KimOllivier

I finally reproduced your errors, and some more using your data. There seems to be something about data in proper metric units that works ok but not US Feet 😉 So I thought that it must be a data corruption. I have not found a fault that would explain that, and now after fiddling around have ended up with some similar results.

I think the fundamental problem is similar to one I hit on another project where the shape object is only referenced in memory when it is extracted with a cursor. Later when we try to access it, there are no methods because it is not a proper Python object. Try doing a dir(shape). If it was a compliant Python object it would return all the methods listed on the Geometry object.

Even with interactive mode you cannot do any useful request on a shape retrieved from the dictionary, you get a crash, even when editing a script and invoking a method. It only works while in the cursor.

The only solution at 9.3 is to extract out all the vertices into an object that Python can use, like a string based GeoJSON dictionary! (At 10 this is practical because there is a function to convert a shape to a GeoJSON string.)

See Bruce Harold's example trying to comprehensively detect the changes between two feature datasets on the Resource Centre.

Although the script is in 10 format, the same applies at 9.3. My shortcut solution for a difference machine was to test the perimeter, length or centroid instead of every vertex because if even one changes it will affect these values.

Can you save gp.Geometry objects in a dictionary? It appears to be intended at 10, but I know there were issues at 9.3. So I don't think that there is much point in proceeding with more testing at 9.3.

ChrisSnyder

Kim,

Seems weird but I ran your edited code on my own FC (with some minor edits for the field name),
It worked twice (flawlessly), and bombed 8 times with one of those bad memory errors that says something like "Can't read memory at BG55XXX05... Application must close."

When it bombs it's always on the 1st attempt to read some aspect of the geometry object (where it prints the .extent property of the geom obj - I edited it to spit out the .xmin BTW).

Even manually I can load the thing into the dictionary, and I start typing the '.' (as in fcDictionary[idx][fieldNameDict[shapeFieldName]].*) denoting I want a property of the geom object, and POOF. But not always, which is really weird.

I will try some other stuff tommorrow (like running it on another machine).

Just curious, but what type of processor do you have? I have a somewhat older (4 or 5 years) AMD Opteron, and I know that it has led to some other quirky ArcGIS-related issues in the past related to raster processing, but nothing like this... i think I am the only goofball that uses an AMD. New machine will be Intel!

KimOllivier

The code example worked for me. I could not get it to fail on a complex polygon featureclass such as the pine forest landcover for a small country. I changed the test to a random shape and printed all vertices and handled possible donuts. I did not compare the vertices with a normal cursor, but they look plausible. There are a few more convenient touches at 10 for iterating over parts and of course the spatial operators but otherwise its the same for 9.3 or 10.
The featurecount is 47,713 and the script took 14 seconds.

THIS CODE LOADS A FC INTO A PYTHON DICTIONARY IN RAM 
# http://forums.arcgis.com/threads/17956-Using-centroid-geometry-to-get-values-of-intersecting-feature#post60607
# edits by Kimo
import arcgisscripting, sys
import random,datetime,pprint
start = datetime.datetime.now()
gp = arcgisscripting.create(9.3)
#THE INPUT FC
fc = "d:/workspace/datetest.gdb/pineforest"
#Process: Build a field list and a mini dictionary to look up the field order (aka index) given a field name
fieldList = gp.listfields(fc)
fieldIndex = 0
fieldNameDict = {}
for field in fieldList:
    fieldNameDict[field.name] = fieldIndex
    fieldIndex = fieldIndex + 1
shapeFieldName = gp.Describe(fc).shapeFieldName #  fix missing definition
pprint.pprint(fieldNameDict)
#Process: Load the fc into a dictionary
fcDictionary = {}
oidFieldName = gp.describe(fc).oidfieldname
searchRows = gp.searchcursor(fc)
searchRow = searchRows.next()
while searchRow:
    oidFieldValue = searchRow.getvalue(oidFieldName)
    fcDictionary[oidFieldValue] = []
    for field in fieldList:
        fcDictionary[oidFieldValue].append(searchRow.getvalue(field.name))
    searchRow = searchRows.next()
del searchRow,searchRows
print
print "Size of dictionary",len(fcDictionary)
for sample in range(10):
    print
    idx = random.randint(0,len(fcDictionary))
    #Print the values 
    print fcDictionary[idx]
    #Print a name field value 
    print fcDictionary[idx][fieldNameDict["LCDB1NAME"]]
    #Print the extent rectangle 
    print fcDictionary[idx][fieldNameDict[shapeFieldName]].extent
    #Print ALL the coordinate pairs including parts and interior rings
    shapeObj = fcDictionary[idx][fieldNameDict[shapeFieldName]]
    for p in range(shapeObj.partCount):
        part = shapeObj.getPart(p)
        pnt = part.next()
        while pnt:
            print p,pnt.X,pnt.Y
            pnt = part.next()
            if not pnt:
                pnt = part.next()
                if pnt:
                    interiorRing = True
                    print "Hole"
print "Completed",datetime.datetime.now() - start

[12139, <geoprocessing describe geometry object object at 0x16A388D8>, 66.0, u'Pine Forest - Closed Canopy', 66.0, u'Pine Forest - Closed Canopy', u'no', 64.70760971, 182.0, u'Rodney District', 16.0, u'Auckland Region', 11067.655349413637, 646803.66190707707]
Pine Forest - Closed Canopy
1811685.23719998 5594544.72269999 1811781.95419998 5594599.59829999 NaN NaN NaN NaN
0 1811685.2372 5594544.9228
0 1811688.5206 5594544.9241
0 1811733.3278 5594561.747
0 1811781.9542 5594599.5983
0 1811780.0872 5594596.3042
0 1811757.3881 5594569.9774
0 1811726.5661 5594549.7439
0 1811689.9306 5594544.7227
0 1811685.2372 5594544.9228
Completed 0:00:14

ChrisMathers

Its something about copying the whole shape field value and not an issue with the polygon class its self. If you the points of a polygon and add them to an array you can store sucessfully the polygon they create or just the array and make the polygon sucessfully on the fly. Making them on the fly however is rather slow. Making the polys from array and storing those takes a couple minutes for a large number of polys but not too long. Comparing a large number of geometries to another large group of geometries is rather slow however using this method.

ChrisSnyder

In my experiences, points, simple lines (a few vertices or so), and simple polygons (again - a few vertices or so), were fine. Try this:

1) Copy the geometry object from a large polygon into a dictionary (it will let you do it).
2) Spit out the vertices

Are they the same coordinates as the original FC (when you read from the FC directly)? In my adventures with this, they were not always the same, and seemed to be corrupted somehow. Also, sometimes running this code would also give me random memory errors and make PythonWin bomb. Maybe I'm doing something wrong. Can you get this to work? Like I said, point layers, and very small polygonns/lines don't seem to be a problem, but test it on a "real" FC:

#THIS CODE LOADS A FC INTO A PYTHON DICTIONARY IN RAM (AT LEAST IT SHOULD)
import arcgisscripting, sys
gp = arcgisscripting.create(9.3)
fc = r"D:\csny490\overlay_20100115\gis_layers\county.gdb\county" #THE INPUT FC

#Process: Build a field list and a mini dictionary to look up the field order (aka index) given a field name
fieldList = gp.listfields(fc)
fieldIndex = 0
fieldNameDict = {}
for field in fieldList:
    fieldNameDict[field.name] = fieldIndex
    fieldIndex = fieldIndex + 1

#Process: Load the fc into a dictionary
fcDictionary = {}
oidFieldName = gp.describe(fc).oidfieldname
searchRows = gp.searchcursor(fc)
searchRow = searchRows.next()
while searchRow:
    oidFieldValue = searchRow.getvalue(oidFieldName)
    fcDictionary[oidFieldValue] = []
    for field in fieldList:
        fcDictionary[oidFieldValue].append(searchRow.getvalue(field.name))
    searchRow = searchRows.next()
del searchRow
del searchRows

#Print the values of OBJECTID = 1
print fcDictionary[1]

#Print the "COUNTY_NM" field value of OBJECTID = 1
print fcDictionary[1][fieldNameDict["COUNTY_NM"]]

#Print the extent rectangle of OBJECTID = 1
print fcDictionary[1][fieldNameDict[shapeFieldName]].extent

#Print the coordinate pairs of OBJECTID = 1 (assuming it is a single part feature)
shapeObj = fcDictionary[1][fieldNameDict[shapeFieldName]]
pointArray = shapeObj.getpart(0)
pointObj = pointArray.next()
while pointObj:
    print str(pointObj.x) + ", " + str(pointObj.y)
    pointObj = pointArray.next()

KimOllivier

But is it corrupted?
This works for me:

# create a dictionary of centroids
sr = gp.Describe(fcAU).SpatialReference
cur = gp.SearchCursor(fcAU,"land = 1")
row = cur.next()
dAU = {}
pt = gp.CreateObject("Point")
while row:
    pt = row.shape.centroid
    au = row.AU_NO
    dAU[au] = pt
    row = cur.next()
del row,cur
print "OD records",len(dAU)

ChrisSnyder

Anyone know why the an ESRI geometry object isnt transferable (at least in an uncorrupted format) to a Python dictionary???

ChrisMathers

I can see what you mean about corrupting the geometry. I can write it all out to a dictionary, but trying to use the geometry object crashes python. I first tried comparing it in the normal manner and that didnt work so I tried making it a polygon object with arcpy.Polygon but that crashed as well. I think im going to try writing out the vertex array as a value in teh dictionary and creating the polygon objects on the fly from those. Sorry if this post has turned into a place to deposit my thoughts. I hope to actually post something that not only works but isnt wildly ineffiecient.

ChrisMathers

I used to do a method that used a spatial join but after upgrading to 10 it started to fail repeatedly without returning any errors. It was as if while loading the data into the parcels a sys.exit() signal was sent. The whole shell restarted instead of a normal crash regardless of what error catching I tried to include.

I think it may save time to loop once through the parcels for each 'i' instead of the other way around. Im testing that now. I may consider the dictionaries as an option if this doesnt work.

The reason I am using networked data is that it is on an SDE and I would rather not have to have the code reconcile the FGDB with the SDE at the start of each run though.

ChrisSnyder

Two ideas that don't really relate to your code structure directly:

1. Copy the FCs locally to a FGDB (e.g. your C:\ drive). Network/SDE = slow, local FGDB = fast.
2. If possible, store the "i" features/attribute values in a dictionary object. This would eliminate the need for the embedded cursor (which I have found to be very slow). In your code I see that for each parcel you open up a whole gaggle of cursors for the "i" FCs. This would be way faster if all this stuff was in a hash table of some sort.

Unfortunately, I have had mixed results loading the geometry objects (the shape field) into a Python dictionary in v9.3. I found that for simple geometry it seemed to always work, but more complicated shapes it seemed to get corrupt somehow... See: http://forums.arcgis.com/threads/9555-Modifying-Permanent-Sort-script-by-Chris-Snyder?p=30010&viewfull=1#post30010 and http://forums.arcgis.com/threads/9555-Modifying-Permanent-Sort-script-by-Chris-Snyder?p=30024&viewfull=1#post30024

Once something is in a dictionary (assumingly uncorrupted) it seems to be the fastest method of emulating embedded cursor behavior (hundreds/thousands of times faster). It has worked out extremely well for me using tabular data, but the shape stuff is different story...

Maybe this is fixed (or at least works better) in v10?

Also, maybe improvements could be made to the looping structure, but I don't know your data well enough... What exactly is your work flow? For each parcel see if it intersects any of the other "i" FCs, right? Would a looping SelectByLocation work the same?

ChrisMathers

Well the biggest problem I was having was that I assumed it was more complicated than it is. To get the geometry object for a row you only have to set a variable equal to the value of the shape field. Then you can use any of the geometry functions available to that geometry type. This is a simple example that iterates over my parcels and looks for a polygon in the FLU layer that its centroid is inside of.

parcel_desc=arcpy.Describe("parcels") #describe the parcels and get the shape field name
parcel_shapefield=parcel_desc.ShapeFieldName
 
flu_desc=arcpy.Describe("flu")#describe the flu and get the shape field name
flu_shapefield=flu_desc.ShapeFieldName 
 
parcel_cursor=arcpy.UpdateCursor("parcels")
for row in parcel_cursor: #iterate over the parcels
    parcel_row_shape=row.getValue(parcel_shapefield) #get the geometry object
    centroid=parcel_row_shape.centroid #get the centroid from the geometry object
    flu_cursor=arcpy.SearchCursor("flu") 
    for flu_row in flu_cursor: #iterate over the flu
        flu_row_shape=flu_row.getValue(flu_shapefield) #get the geometry object
        if centroid.within(flu_row_shape)==True: #compare the two geometries
            print "RENUM: %s  FLU: %s" % (parcel_row.A1RENUM, flu_row.FLU_CODE) #do some junk with the two rows
    del flu_cursor
del parcel_cursor

ChrisSnyder

Cool! Makes you wonder why there's hardly any ESRI help and/or code examples items with this...

P.S. I am wanting to adapt some very inefficient v9.3 code to use this stuff (arcpy.geometry). Can't wait! 😄

Need to get the "important" projects done so I can get into v10 finally 😛

ChrisMathers

I have figured it out! Will post code later today once I get it nailed down.

ChrisMathers

After some playing around I know how to make a point-geometry object. All you do is get the centroid from shape field.

parceldesc=arcpy.Describe('parcels')
parcelShapeField=parceldesc.ShapeFieldName
parcel_cursor=arcpy.UpdateCursor('parcels')
for parcel_row in parcel_cursor:
    parcel_row_shape=parcel_row.getValue(parcelShapeField) #
    centroid=parcel_row_shape.centroid

and to compare that to another geometry I just need to do

centroid.crosses(other geometry object)

I dotn know how to get a polygon geometry from a row in a way taht is useful though. I can get a convex hull and an extent envelope but neither are helpful.

ChrisMathers

Thanks Chris. I think this is what I want to do. Granted this is extremely metaish code lol.

For update_row in update_cursor:
  geo_update_row=create a geometry object for this row's label point
  for search_row in search_cursor:
    geo_search_row=create a geometry object for this row
    if geo_update_row.crosses(geo_search_row) == True:
      update_row.something = search_row.something
      update_cursor.updaterow(update_row)
    del geo_search_row
    del geo_update_row

Thats the general idea. Now I just have to learn geometry objects.

ChrisSnyder

I can't find any documentation in the help about that stuff (arcpy.geometry)...

You've probaly seen it, but here's a good example script though:
http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=7DCD71C0-1372-4ECD-2860-A023289A7EC5

Funny that in the code they acknowledge some bug issues and have a work around...

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.