5

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

I would be very grateful for any help.

  • 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 Feb 4 '15 at 22:13
6

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.

| improve this answer | |
  • This is nice. I always forget about the distanceTo() method. In my version, I thought it would be good to bring all the attributes over as well because I tested on Parcels. Since they have coincidental boundaries, it made it easier to determine which point corresponded to which polygon. – crmackey Feb 5 '15 at 14:23
1

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)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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