Select to view content in your preferred language

Using centroid geometry to get values of intersecting feature

4064
19
11-23-2010 07:08 AM
ChrisMathers
Deactivated User
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.
Tags (2)
0 Kudos
19 Replies
ChrisSnyder
Honored Contributor
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-A0232...

Funny that in the code they acknowledge some bug issues and have a work around...
0 Kudos
ChrisMathers
Deactivated User
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.
0 Kudos
ChrisMathers
Deactivated User
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.
0 Kudos
ChrisMathers
Deactivated User
I have figured it out! Will post code later today once I get it nailed down.
0 Kudos
ChrisSnyder
Honored Contributor
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 😛
0 Kudos
ChrisMathers
Deactivated User
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
0 Kudos
ChrisSnyder
Honored Contributor
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&viewfu... and http://forums.arcgis.com/threads/9555-Modifying-Permanent-Sort-script-by-Chris-Snyder?p=30024&viewfu...

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?
0 Kudos
ChrisMathers
Deactivated User
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.
0 Kudos
ChrisMathers
Deactivated User
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.
0 Kudos