4

I have a polygon and the polygon's centroid. I want to select the most distant vertex from centroid (picture: point A). I used arcpy.FeatureVerticesToPoints_management to create a feature class containing points generated from polygon, but I don't know what to do next.

Could you help me?

Also I want to find intersection point (point B) of the outline of polygon and line created by centroid and the most distant point (line OA).

Example

1
  • As many of the answers have pointed out, you really have two questions here: (1) the most distant vertex, (2) finding intersection point B. I suggest that you break off #2 into a new Question, to give it the focus and attention it deserves (it is a very interesting problem in its own right!) :)
    – Erica
    Commented Feb 4, 2015 at 22:13

2 Answers 2

7

Here's a pared down version of @crmackey's answer. The polygon layer is called 'POLY1', and should be the only thing you need to change to get an output point file of farthest vertices - it creates centroids on-the-fly:

>>> points = []
>>> with arcpy.da.SearchCursor("POLY1",['SHAPE@']) as cursor:
...     for row in cursor:
...         centroid = row[0].centroid
...         dist = 0
...         for part in row[0]:
...             for pnt in part:
...                 cent_vert_dist = arcpy.PointGeometry(pnt).distanceTo(centroid)
...                 if cent_vert_dist > dist:
...                     dist = cent_vert_dist   
...                     far_point = arcpy.PointGeometry(pnt)
...         points.append(far_point)
...         
>>> arcpy.CopyFeatures_management(points,'in_memory\points')

Back-tracking to the intersection point opposite the farthest vertex is possible, but will require some additional trigonometry that I'm not prepared to get into, atm.

0
2

This worked for me. The script will create an output point feature class that returns the point that is farthest from the centroid for all polygons:

import arcpy
import os
import sys
import traceback
import math
from datetime import datetime as d
arcpy.env.overwriteOutput = True

def Message(msg):
    print str(msg)
    arcpy.AddMessage(msg)

def findDistance(a,b):
    x = abs(a[0] - b[0])
    y = abs(a[1] - b[1])
    return math.sqrt((x*x) + (y*y))

def iter_geom(g):
    for i in xrange(g.partCount):
        yield i, g.getPart(i)

def polyPoints(in_polys, out_points):
    """function to find the minimum or maximimum distance between points

    Required:
        in_polys -- source points for near analysis
        out_points -- output nearest or farthest point from centroid
    """

    startTime = d.now()
    # grab xy coords
    point_dict, cent_dict = {}, {}
    with arcpy.da.SearchCursor(in_polys, ['OID@','SHAPE@']) as rows:
        for row in rows:
            points = []
            for i, part in iter_geom(row[1]):
                for pt in iter(lambda: part.next(), None):
                    points.append((pt.X, pt.Y))
            point_dict[row[0]] = points
            cent_dict[row[0]] = (row[1].centroid.X, row[1].centroid.Y)

    # grab attributes
    desc = arcpy.Describe(in_polys)
    fields = [f.name for f in desc.fields if not f.required]
    with arcpy.da.SearchCursor(in_polys, ['OID@'] + fields) as rows:
        att_dict = {r[0]: r[1:] for r in rows}



    # create dictionary to find nearest or farthest vertex from centroid
    newpt_dict = {}
    for oid, points in point_dict.iteritems():
        c = cent_dict[oid]
        dist_dict = {coords: findDistance(c,coords) for coords in points}
        query_pt = max(dist_dict.values())
        for k,v in dist_dict.iteritems():
            if v == query_pt:
                newpt_dict[oid] = (k, v)

    # create output table
    path, name = os.path.split(out_points)
    if arcpy.Exists(out_points):
        arcpy.Delete_management(out_points)
    sr = desc.spatialReference
    arcpy.CreateFeatureclass_management(path, name, 'POINT', in_polys, 'DISABLED', 'DISABLED', sr)
    arcpy.management.AddField(out_points, 'POLY_ID', 'LONG')
    arcpy.management.AddField(out_points,'POINT_DIST','DOUBLE')
    add_fields = ['POLY_ID', 'POINT_DIST']
    ifields = ['SHAPE@'] + fields + add_fields
    with arcpy.da.InsertCursor(out_points, ifields) as irows:
        for oid, pt_dist in newpt_dict.iteritems():
            all_atts = (pt_dist[0],) + att_dict[oid] + (oid, pt_dist[1])
            irows.insertRow(all_atts)
    Message('Created: %s' %os.path.basename(out_points))
    Message('(Elapsed time: %s)' %(str(d.now() - startTime)[:-3]))
    return out_points

if __name__ == '__main__':

##    # stand alone
    polys = r'C:\TEMP\Utility\Utilities.gdb\Parcels_small'
    out = polys + '_farPoints'
    polyPoints(polys, out)

Not the answer you're looking for? Browse other questions tagged or ask your own question.