1

I'm trying to develop a function which will find the longest straight line length that can be drawn within any of the polygons in a feature class. (To keep things simple, I'm only including lines that can be drawn between any two vertices of a polygon.)

My logic so far, is (for each polygon):

  • create Polyline geometries between all possible pairs of vertices of the polygon
  • for each Polyline created, test if it is fully within the polygon (eg, does not cross the boundary of a concave part of the polygon)
  • if it is fully within the polgyon and it's length is the largest found so far, update the maximum length variable.

Unfortunately, my arcpy line.within(poly) always returns false, even though the lines are created from pairs of points of that polygon.

I've commented out the test of line length to ignore that for now, but the "WITHIN" message (below) is never printed and the function always returns 0 (zero).

What I'm doing wrong in this code, below?

I can confirm that the function does create over 72,000 polylines for my test polygon, and that they all have realistic looking (greater than zero) lengths. (Confirmed by using arcpy.AddMessage() for these values.)

def PolyLength(polyFC):
    maxLength = 0
    with arcpy.da.SearchCursor(polyFC, ["SHAPE@"]) as cur:
        for row in cur:
            poly = row[0]
            featureMaxLength = 0
            points = []
            lines = []
            for part in poly:
                points.extend(part)
                arcpy.AddMessage("Points count:  {}".format(len(points)))
            for a in points:
                for b in points:
                    if not a.equals(b):
                        lines.append(arcpy.Polyline(arcpy.Array([a,b])))
            arcpy.AddMessage("Lines count:  {}".format(len(lines)))
            for line in lines:
                #if line.length > featureMaxLength and line.within(poly):
                if line.within(poly):
                    featureMaxLength = line.length
                    arcpy.AddMessage("WITHIN and greater length:  {}".format(line.length))
                else:
                    arcpy.AddMessage("NOT WITHIN or NOT greater length:  {}".format(line.length))
            if featureMaxLength > maxLength:
                maxLength = featureMaxLength
    return maxLength
0

2 Answers 2

2

Method "within" works fine, because this script

import arcpy
import itertools as itt
lines="LINES"
pgons="PGONS"
d=arcpy.Describe(pgons)
SR=d.spatialReference    
g=arcpy.Geometry()
curT=arcpy.da.InsertCursor(lines,"Shape@")
with arcpy.da.SearchCursor(pgons,"Shape@") as cursor:
    for row in cursor:
        pgon=row[0]
        gList=arcpy.FeatureVerticesToPoints_management(pgon, g, "ALL")
        lMax=-1
        for pair in itt.combinations(gList,2):
            arr=arcpy.Array()
            for pG in pair:arr.add(pG.firstPoint)
            line=arcpy.Polyline(arr,SR)
            lCur=line.length
            if lCur<lMax:continue
            if not line.within(pgon):continue
            lMax=lCur
            longest=line
        arcpy.AddMessage(longest.length)
        curT.insertRow((longest,))

produces expected results:

enter image description here

It is the way you are constructing the list of vertices.

0
1

As demonstrated in the other answer (but not fully explained there), the solution was to make sure that each line is created with the same spatial reference as the polygons.

(I had actually tried this earlier, and got an error stating that spatial_reference is a read-only attribute. Turns out I had specified it incorrectly, and of course that error message had led me down the garden path!)

So the specific answer to this question for my specific code issue is...

Replace this line:

lines.append(arcpy.Polyline(arcpy.Array([a,b])))

with this:

sr = arcpy.Describe(polyFC).spatialReference
...
...
lines.append(arcpy.Polyline(arcpy.Array([a,b]), sr))

My complete updated code is:

def lineLen(line):
    return line.length

def PolyLength(polyFC):
    maxLength = 0
    sr = arcpy.Describe(polyFC).spatialReference
    with arcpy.da.SearchCursor(polyFC, ["SHAPE@"]) as cur:
        for row in cur:
            poly = row[0]
            featureMaxLength = 0
            points = []
            lines = []
            for part in poly:
                points.extend(part)
                arcpy.AddMessage("Points count:  {}".format(len(points)))
            for a in points:
                for b in points:
                    if a and b and not a.equals(b):
                        lines.append(arcpy.Polyline(arcpy.Array([a,b]), sr))
            arcpy.AddMessage("Lines count:  {}".format(len(lines)))
            lines.sort(key=lineLen, reverse=True)
            for line in lines:
                if line.within(poly):
                    featureMaxLength = line.length
                    arcpy.AddMessage("WITHIN:  {}".format(line.length))
                    break
            if featureMaxLength > maxLength:
                maxLength = featureMaxLength
    return maxLength

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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