0

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
  • Try to add spatial reference in arcpy.Polyline(). Also consider itertools to iterate through points easily. – FelixIP Aug 13 at 5:26
  • Also within has second optional parameter, might play a role with points on a boundary. – FelixIP Aug 13 at 5:41
  • The SpatialReference of the two shapes must be identical, or all comparison operators will always return False. Use the .spatialReference of the Describe on the source feature class in the Polyline constructor. – Vince Aug 13 at 10:13
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.

  • 1
    This works because the spatial_reference parameter of the Polyline constructor is not ignored (it's only optional if you want all comparison operators to fail and for the coordinates to be truncated to assumed millimeters). – Vince Aug 13 at 10:20
  • See @Vince comments on what is wrong in your workflow, it is failure to define spatial reference. – FelixIP Aug 13 at 22:19
  • This answer sort of includes the solution, but is lacking in explanation. It is the way the lines are defined (with/without spatial reference), NOT the way the vertices are constructed that was causing the problem. I will accept it if it gets updated with a suitable explanation. Otherwise, I've posted a more precise answer, including explanation. Thanks for your help. – Son of a Beach Aug 13 at 23:27
  • Post yours. Don't forget to mention that your follow up on my first comment is sort of misleading, to say the least. I see it is deleted already. Nicely done mate. – FelixIP Aug 13 at 23:52
0

As explained in comments, and 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, privacy policy and cookie policy

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