Arcpy centroid within polygon

[DarrenWiens2]

I am attempting to make a script tool that  takes polygon input, and outputs centroid point (within the polygon). This knowledge base article seems to be incorrect, in that the SHAPE@XY token returns the center of gravity for the polygon, not the centroid within the polygon. Likewise, the centroid property in the Polygon class, used in the code below, returns the center of gravity rather than centroid within polygon. Any ideas?

# import libraries
import arcpy, os

# set input/output parameters
polyFC = arcpy.GetParameterAsText(0)
outCentroids = arcpy.GetParameterAsText(1)

# set overwrite environment
arcpy.env.overwriteOutput = True

# if the output file does not exist, create it. Add "ORIG_ID" field.
if not arcpy.Exists(outCentroids):
    arcpy.CreateFeatureclass_management(os.path.dirname(outCentroids),
                                    os.path.basename(outCentroids),
                                    "POINT",
                                    polyFC,
                                    "",
                                    "",
                                    polyFC)
    arcpy.AddField_management(outCentroids,'ORIG_ID', 'LONG')

# collect all of the polygon centroids into an array
pointArray = []
for row in arcpy.da.SearchCursor(polyFC, ["SHAPE@",'*','OID@']):
    rowArray = []
    rowArray.append(row[0].centroid)
    for field in range(1,len(row)):
        rowArray.append(row[field])
    pointArray.append(rowArray)
del row

# write the centroids to the output file
cursor = arcpy.da.InsertCursor(outCentroids, ["SHAPE@",'*'])
for point in pointArray:
    arcpy.AddMessage(point)
    cursor.insertRow(point)
del cursor

Example output below:

[OwenEarley]

Have you tried using the polygon labelPoint? - this is supposed to be within the polygon:

The point at which the label is located. The labelPoint is always located within or on a feature.

The documentation suggests that the labelPoint should be returned by centroid if the true centroid is not within the feature but this does not seem to be the case in your example.

[DanPatterson_Retired]

labelpoint is what you are after ... me-thinks

[DarrenWiens2]

Thanks for the replies, but still "no."

It is my understanding that "centroid" is supposed to return the labelPoint in the event that the first returned point does not fall inside the polygon.

Here is the updated, still not working, code:

# import libraries
import arcpy, os


# set input/output parameters
polyFC = arcpy.GetParameterAsText(0)
outCentroids = arcpy.GetParameterAsText(1)


# set overwrite environment
arcpy.env.overwriteOutput = True


# if the output file does not exist, create it. Add "ORIG_ID" field.
if not arcpy.Exists(outCentroids):
    arcpy.CreateFeatureclass_management(os.path.dirname(outCentroids),
                                    os.path.basename(outCentroids),
                                    "POINT",
                                    polyFC,
                                    "",
                                    "",
                                    polyFC)
    arcpy.AddField_management(outCentroids,'ORIG_ID', 'LONG')


# collect all of the polygon centroids into an array
pointArray = []
for row in arcpy.da.SearchCursor(polyFC, ["SHAPE@",'*','OID@']):
    rowArray = []
    rowArray.append(row[0].labelPoint)
    for field in range(1,len(row)):
        rowArray.append(row[field])
    pointArray.append(rowArray)
del row


# write the centroids to the output file
cursor = arcpy.da.InsertCursor(outCentroids, ["SHAPE@",'*'])
for point in pointArray:
    arcpy.AddMessage(point)
    cursor.insertRow(point)
del cursor

[DanPatterson_Retired]

labelpoint  The point at which the label is located. The labelPoint is always located within or on a feature.

centroid     The true centroid if it is within or on the feature; otherwise, the label point is returned. Returns a point object.

truecentroid   The center of gravity for a feature.

From the Geometry class, so I guess the most appropriate would be to specify truecentroid or labelpoint if centroid isn't up to snuff

[DarrenWiens2]

Thanks, Dan, although I'm reading the exact same help files.

...and I get the same result for centroid, trueCentroid, and labelPoint.

I'll admit I am not fully comfortable accessing geomtries with tokens - could this be the source of the issue? I've produced the same result (centroid outside of polygons) with both SHAPE@ and SHAPE@XY.

[DanPatterson_Retired]

Darren...I figured as much...this was more for the benefit of those that aren't aware of the power of the help files and how one can type in a keyword such as truecentroid and be given a list of appropriate links.  Similarily, on GeoNet, you can use the same truecentroid (2) and come up with a list of links some of which show that there have been errors in polygon geometry in the past  some have found other issues and if you search hard enough...you can find yourself again proving that you can never really get lost using the help files or GeoNet

[OwenEarley · Accepted solution]

After running a similar test I am getting the expected results - they match what the documentation says should be happening:

cursor = arcpy.da.InsertCursor(outCentroids, ["SHAPE@",'ORIG_ID'])
for row in arcpy.da.SearchCursor(polyFC, ["SHAPE@",'OID@']): 
    cursor.insertRow((row[0].centroid, row[1]))
del row

Is the spatial reference defined on your input data?

[DanPatterson_Retired]

If you check my thread about geometry errors, you will find that, that was the conclusion, is SR is not set, single precision math is used otherwise double, and I am certain that it has ramifications throughout the geom... chain

[DarrenWiens2]

Thank you! It goes to show that taking a step back can lend some clarity. In the end, I had misdiagnosed my problem, but your example definitely helped.

The problem was in how I was re-writing the geometry. The major difference in my working script (below) and yours is that mine takes along all (an unknown number) of fields and writes them to the new geometry. I had been assigning the row[0].centroid coordinates to the SHAPE@XY token, but it was later overwritten by the shape value coming in with all fields ('*'). At least, I think that's what was going on...

# import libraries
import arcpy, os

# set input/output parameters
polyFC = arcpy.GetParameterAsText(0)
outCentroids = arcpy.GetParameterAsText(1)

# set overwrite environment
arcpy.env.overwriteOutput = True

# if the output file does not exist, create it. Add "ORIG_ID" field.
if not arcpy.Exists(outCentroids):
    arcpy.CreateFeatureclass_management(os.path.dirname(outCentroids),
                                    os.path.basename(outCentroids),
                                    "POINT",
                                    polyFC,
                                    "",
                                    "",
                                    polyFC)
    arcpy.AddField_management(outCentroids,'ORIG_ID', 'LONG')

# create an InsertCursor containing all the fields in ourCentroids, which are all the fields in polyFC plus "ORIG_ID"
cursor = arcpy.da.InsertCursor(outCentroids, ['*'])

# read all features in polyFC
for row in arcpy.da.SearchCursor(polyFC, ["SHAPE@",'*','OID@']):
     # create an array to hold a new, modified row
    rowArray = []
     # read all the fields except "SHAPE@"
    for fieldnum in range(1,len(row)):
        # if this is the second field [i.e. SHAPE], replace it with the centroid coordinates
        if fieldnum == 2:
            rowArray.append(row[0].centroid)
        else:
            rowArray.append(row[fieldnum])
    # write the new row to the cursor
    cursor.insertRow(rowArray)
del row