Monthly Archives: September 2014

Arcpy: Polygons to Centroids (within Polygons, and with all Attributes)

Here’s a fairly simple python script that creates centroid points (constrained to fall within polygons) and attaches all the attributes (plus a link field, “ORIG_ID”) from the polygons to the centroids. Notably, it uses cursors to read raw geometry objects – no geoprocessing tools (except CreateFeatureclass). Also, while Advanced licensees may use the Feature to Point tool for this procedure, the method presented here is completely free, using no extensions or non-basic licensing – I don’t know of another free way to achieve this in ArcGIS. Wouldn’t it be easier to run a spatial join on the points to attach the attributes? Yes, but spatial join is slow. In my test of 1000 polygons, the cursor method completed in about 2 seconds, while adding a spatial join call extended the run time to more than 10 seconds. This could mean big savings on large datasets.

  1. # import libraries
  2. import arcpy, os
  3. # set input/output parameters
  4. polyFC = arcpy.GetParameterAsText(0)
  5. outCentroids = arcpy.GetParameterAsText(1)
  6. # set overwrite environment
  7. arcpy.env.overwriteOutput = True
  8. # if the output file does not exist, create it. Add “ORIG_ID” field.
  9. if not arcpy.Exists(outCentroids):
  10.     arcpy.CreateFeatureclass_management(os.path.dirname(outCentroids),
  11.                                     os.path.basename(outCentroids),
  12.                                     “POINT”,
  13.                                     polyFC,
  14.                                     “”,
  15.                                     “”,
  16.                                     polyFC)
  17.     arcpy.AddField_management(outCentroids,‘ORIG_ID’‘LONG’)
  18. # create an InsertCursor containing all the fields in ourCentroids, which are all the fields in polyFC plus “ORIG_ID”
  19. iCursor = arcpy.da.InsertCursor(outCentroids, [‘*’])
  20. # read all features in polyFC
  21. with arcpy.da.SearchCursor(polyFC, [“SHAPE@”,‘*’,‘OID@’]) as sCursor:
  22.     for row in sCursor:
  23.          # create an array to hold a new, modified row
  24.         rowArray = []
  25.          # read all the fields except “SHAPE@”
  26.         for fieldnum in range(1,len(row)):
  27.             # if this is the 2nd field [i.e. SHAPE], replace it with the centroid coordinates
  28.             if fieldnum == 2:
  29.                 rowArray.append(row[0].centroid)
  30.             else:
  31.                 rowArray.append(row[fieldnum])
  32.         # write the new row to the cursor
  33.         iCursor.insertRow(rowArray)
  34. del iCursor

Looping through the fields in this way (lines 29 – 36) feels cumbersome – if anyone knows a better way to do this, please let me know!

————————————————-

EDIT (January 19, 2016): Following some discussion on LinkedIn’s GIS and Geography group regarding this post, here is an example of why you must use the SHAPE@ token (and Polygon.centroid) rather than either the SHAPE@XY or SHAPE@TRUECENTROID tokens to guarantee that the point falls inside the given polygon:

… points = []
… with arcpy.da.SearchCursor(“MY_FEATURE_LAYER”,[‘SHAPE@’,’Shape@XY’,’SHAPE@TRUECENTROID’]) as cursor:
…   for row in cursor:
…     points.append(arcpy.PointGeometry(row[0].centroid)) # SHAPE@
…     points.append(arcpy.PointGeometry(arcpy.Point(row[1][0],row[1][1]))) # SHAPE@XY
…     points.append(arcpy.PointGeometry(arcpy.Point(row[2][0],row[2][1]))) # SHAPE@TRUECENTROID
… arcpy.CopyFeatures_management(points, r’in_memory\mypoints’)

centroid_tokens

 

Using ArcObjects in Python

Recently, I found the need for more fine-grained control over a map document than the Arcpy library alone could afford, so I turned to the comtypes library (downloaded here) and associated Snippets file (ArcGIS 10.1 version here), which allows access to VBA-based ArcObjects. The Mark Cederholm presentation I followed is from the 2010 ESRI Dev Summit, so it is possible that this method is somewhat outdated, but it does exactly what I need.

Here is a simple example showing how to add some guides to the layout rulers in a map document, which is not possible solely using Arcpy:

# Import the modules
from comtypes.client import GetModule, CreateObject
from Snippets import GetStandaloneModules, InitStandalone

# First time through, need to import the “StandaloneModules”. Can comment out later.
# GetStandaloneModules()
InitStandalone()

# Get the Carto module
esriCarto = GetModule(r”C:\Program Files (x86)\ArcGIS\Desktop10.1\com\esriCarto.olb”)

# Create a map document object
mxdObject = CreateObject(esriCarto.MapDocument, interface=esriCarto.IMapDocument)

# Create new mxd file
mxdObject.New(r”map.mxd”)

# Create some new layout guides
for i in range(1,21):
    mxdObject.PageLayout.VerticalSnapGuides.AddGuide(i)
    mxdObject.PageLayout.HorizontalSnapGuides.AddGuide(i)

# Change current view to layout
mxdObject.SetActiveView(mxdObject.PageLayout)

# Save the mxd
mxdObject.Save()

# Delete the map document object
del mxdObject

…and the result:

layout