Longest Line Within a Polygon

12074
5
04-16-2014 01:45 PM
RafaelRuas_Martins
New Contributor III
Hi,
I'd like to model the longest possible line within a (irregular ou regular) polygon feature.
The line should be a 100% inside polygon!!!
Someone knows a tool for that task?
I am using ArcGIS for Desktop 10.2.1.
Tks.
Rafael
0 Kudos
5 Replies
WilliamCraft
MVP Regular Contributor
It sounds like you are wanting to select the longest line feature which is completely within a specific polygon.  Are you looking for a model, a python script, or the instructions for how to do this manually within ArcMap?  In general, I think you would need several tools to address this requirement.  Here are the high-level steps:

** Assumes you have a line feature class and a polygon feature class **


  1. Select a single polygon feature, either manually or based on some attribute criteria. 

  2. Select by Location: select lines from the line feature class which are completely within the polygon feature class (the previously-selected polygon from the step above will honor the selection set of 1 feature)

  3. Find the longest feature - this can be handled multiple ways:


  • Select by Attribute: select lines from the currently-selected features of the line feature class using MAX(SHAPE.LENGTH) or MAX(SHAPE_LENGTH) depending on your type of data source

  • Open the attribute table of the line feature class, choose the option at the bottom to view only the selected features, and sort the SHAPE.LENGTH field (or SHAPE_LENGTH field) in descending order

  • Create a new field of DOUBLE (38,8) type and use the Calculate Geometry tool (right click the field from the attribute table) or the Calculate Field geoprocessing tool to calculate the true length of the line features based on the units from the coordinate system of the feature class. 


There are likely more ways to handle step 3 than the bulleted options I specified above.  Let me know more about your requirements and I can be more specific about the approach.
0 Kudos
KimOllivier
Regular Contributor II
Just generate all the lines by connecting vertex pairs of the polygon and find the longest.
Try to avoid iterating around each point combination, think of a process that handles all lines at once.

You can leave out the first and last pair of each point iteration because they will be sides, and you only have to do half because the direction does not matter. You can use python collection permutations to get a unique list of pairs. Then turn them into a list of polylines.

Because a polygon may not be regular, say a banana shape, then you will need a filter to exclude potential lines crossing the boundary for my simple first stab, but there is a test [ a_poly.overlaps(a_line) ] for that in shape objects. Even better you could test them all in one step, using SelectLayerByLocation(..."WITHIN_CLEMENTINI"...) to select remaining candidates to find the max(length). All the tools will take a geometry list instead of a featureclass so it is fast, in memory and efficient enough to iterate over a large number of polygons.
0 Kudos
markdenil
Regular Contributor II
I used to do this by building a near table on the vertices and
selecting the pair with the maximum distance; noting that an
isocolese triangle would have two such pairs,
so attention paid to getting the correct pair
(and not the two points across the base)
was worthwhile.

There is a new tool called
MinimumBoundingGeometry In the Management toolbox under Features.

The CONVEX_HULL option adds (among others) the fields
MBG_Length�??The longest distance between any two vertices of the convex hull; these vertices are called antipodal pairs or antipodal points. (It may be found between more than one pair of vertices, but the first found will be used.)
MBG_APodX1�??The x coordinate of the first point of the antipodal pairs.
MBG_APodY1�??The y coordinate of the first point of the antipodal pairs.
MBG_APodX2�??The x coordinate of the second point of the antipodal pairs.
MBG_APodY2�??The y coordinate of the second point of the antipodal pairs.

It does NOT, however, test for interiority...

It does sound like constructing a set of geometry object lines;
testing them WITHIN_CLEMENTINI,
and picking the longest may be what you have to do.
0 Kudos
KimOllivier
Regular Contributor II
I couldn't resist hacking up a script to see what happened.
There were of course a few issues.
The MakeFeatureLayer tool would not take a geometry list, so I resorted to using in_memory featureclasses
The WITHIN_CLEMENTINI query worked fine interactively in ArcMap but chokes on a polygon with more than a dozen vertices in a script. So a new workaround is needed for that. I test each line as it is generated, and then I could also test its length cumulatively so would only have to test for contains if it is a replacement candidate and I only need to keep the longest current line.

This reduced the time per polygon to about 2 seconds. Fairly slow, but it will handle very complex polygons now.

Second Hack:

#-------------------------------------------------------------------------------
# Name:        LongestLine.py
# Purpose:     find the longest inside line across a polygon
#              Limitations: only single part polygons, no donuts
# Author:      kimo
#
# Created:     29/04/2014
#              30/04/2014 changed crossing test to each features
#               exclude triangles (use trigonometry for this degenerate case)
#               longest line only tested between vertices may need to densify first
#
# Copyright:   (c) kimo 2014
# Licence:     Creative Commons 3.0 New Zealand
#-------------------------------------------------------------------------------
import arcpy
import sys
import itertools
import datetime
gdb = sys.argv[1]
fcPoly = 'forest'
fcOut = 'fan'
arcpy.env.workspace = gdb
arcpy.env.overwriteOutput = True
begintime = datetime.datetime.now()
desc = arcpy.Describe(fcPoly)
sr = desc.spatialReference
arcpy.env.outputCoordinateSystem = sr
if arcpy.Exists(fcOut):
    arcpy.management.Delete(fcOut)
arcpy.management.CreateFeatureclass(gdb,fcOut,"POLYLINE",spatial_reference=sr)
arcpy.management.AddField(fcOut,"ORIGID","LONG")
icur = arcpy.da.InsertCursor(fcOut,['ORIGID','SHAPE@'])
with arcpy.da.SearchCursor(fcPoly,["OBJECTID","SHAPE@"]) as cur:
    for row in cur:
        gPoly = row[1]
        id = row[0]
        print "Begin id {} parts {} points {} type {}".format(id,gPoly.partCount,gPoly.pointCount,gPoly.type)
        if gPoly.partCount == 1 and gPoly.pointCount > 4: # exclude triangles
            starttime = datetime.datetime.now()
            lstV = []
            for part in gPoly:
                for v in part:
                    lstV.append(v)
            oCount = 0
            pCount = 0
            rCount = 0
            maxlength = -1
            # use itertools to get unique combinations of pairs
            # leave in lines coincident with boundary segments
            # because they will always be less than max length if > triangle
            for vPair in itertools.combinations(lstV[:-1],2): # exclude last duplicated poly vertex
                aLine = arcpy.Polyline(arcpy.Array(vPair),sr)
                # dont even bother doing spatial test if already shorter for speed
                if aLine.length > maxlength :
                    if gPoly.contains(aLine): # not lines crossing edges
                        mLine = aLine
                        maxlength = aLine.length
                        pCount+=1
                    else :
                        oCount+=1
                else:
                    rCount+=1
            print "Max swapped {} Overlap rejects {} Already shorter {}".format(pCount,oCount,rCount)
            rec = (id,mLine)
            icur.insertRow(rec)
            msg = "{} {} time {}".format(id,maxlength,datetime.datetime.now() - starttime)
            print msg
            arcpy.AddMessage(msg)

del icur

print "done",datetime.datetime.now() - begintime



I added the original Objectid to be able to relate back to each polygon, geometries do not allow me to keep the original ID so I keep it as a python variable.
0 Kudos
KimOllivier
Regular Contributor II
Since writing my solution I have found an ArcView 3.x solution that shows I have not covered all the cases. You could have a line running through a pair of vertices on a polygon with rentrant sides that extent out each way to the polygon boundary between vertices. I assumed that the longest line would always be coincident with vertices at the ends.

http://www.jennessent.com/arcview/longest_lines.htm

Well done to Jeff Jenness for handling this case in Avenue.
0 Kudos