#——————————————————————————-
# 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 = ‘E:/Exeter/Research/AHRC_LongestLine/’
arcpy.env.workspace = gdb
fcPoly = ‘E:/Exeter/Research/AHRC_LongestLine/GB_area_10m_Thames_Forth.shp’
arcpy.RepairGeometry_management(fcPoly)
fcOut = ‘GB_LL_Thames_Forth_out.shp’
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