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 :p.

[ChrisSnyder]

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

[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]

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

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

[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

[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]

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.

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

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

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