Nothing p*sses me off more than to clip data in a geographic coordinate system and then project to UTM or MTM since the clipped polygon doesn't account for the latitudinal curvature. The following helper script and the main script will allow you to create a tool which will do this. Densify.py and its helper script takes an input and creates an output with the points that form the line segments/polygon boundaries with an increased point density. At some stage I will upload a full toolset to the ESRI website. So on to the helper script
'''
DensifyHelper.py
Author
Dan Patterson
Dept of Geography and Environmental Studies
Carleton University, Ottawa, Canada
Dan_Patterson@carleton.ca
Date March 2008
Modified June 2013
Purpose
helper scripts for Densify.py which densifies points forming
line segments. currently requires point input for testing
'''
#-------------------------------------------------------------
import os, sys, math
#constants
degToRad = math.pi/180.0
radToDeg = 180.0/math.pi
def p1p2Angle(pnt1, pnt2):
#Requires two distinctly different points, returns the angle
dY=pnt2[1] - pnt1[1]
dX=pnt2[0] - pnt1[0]
if dY == 0.0 and dX == 0.0:
return -9999.99
theAngle = math.atan2(dY,dX)*radToDeg
return theAngle
def polyLength(pnts):
#returns the total length of a poly* feature given its points
dist = 0.0
for i in range(1,len(pnts)):
pnt1 = pnts[i-1]; pnt2 = pnts[i]
dist += pntDist(pnt1, pnt2)
return dist
def pntCenter (p1,p2):
#returns the centre point
centPnt=[(p1[0]+p2[0])/2,(p1[1]+p2[1])/2]
return centPnt
def pntDist(pnt1, pnt2):
dx = pnt2[0] - pnt1[0]
dy = pnt2[1] - pnt1[1]
dist = math.hypot(dx,dy)
return dist
def pntXY(pnt):
#Gets X,Y coordinates given a point object
# Returns a string
XY= str(pnt.X) + ", " + str(pnt.Y)
return XY
def transRotatePnts(pnts, pntCent, angle):
#requires a point pair, point center and an angle in degrees
# [[0,0],[1,1] ] point pair
# 45 rotation angle
X0=pntCent[0]; Y0=pntCent[1]
N=len(pnts)
#translate and rotate shapes
angle = angle*(-1.0) #reverse the rotation
cosXY = math.cos(degToRad * angle)
sinXY = math.sin(degToRad * angle)
rotPnts = []
for j in range(0,N):
X1 = pnts[j][0] - X0; Y1 = pnts[j][1] - Y0
X = (X1 * cosXY) - (Y1 * sinXY)
Y = (X1 * sinXY) + (Y1 * cosXY)
pnt = [X,Y]
rotPnts.append(pnt)
return rotPnts
def polyToPnts(shapeList):
#requires the shape type, shapeList is OID and a shape
#returns:
# list of lists of all points
# ([OIDval, X, Y]), None is used for interior rings and parts
# simply decide which ones you need upon retrieval
lineList = []
OIDval = shapeList[0] #object ID
aShape = shapeList[1] #the shape
i = 0
aPart = 0
while aPart < aShape.partCount:
anArray = aShape.getPart(aPart)
pnt = anArray.next()
j = 0
while pnt:
X = pnt.X; Y = pnt.Y
aList = [OIDval, X, Y]
lineList.append(aList)
pnt = anArray.next()
if not pnt:
lineList.append(None)
pnt = anArray.next()
aPart = aPart + 1
return lineList
def densifyPoly(pnts, stepBy):
#requires a point list ([X, Y] and a step for densification
densified = []
remainder = 0.0
for i in range(1,len(pnts)):
pnt1 = pnts[i-1]; pnt2 = pnts[i]
if (pnt2 == None):
densified.append(pnt1)
densified.append(None)
continue
elif (pnt1 == None):
continue
else:
#str(pnt1) != str(pnt2): #added last bit
densified.append(pnt1)
dist = pntDist(pnt1, pnt2)
(remainder, steps) = math.modf(dist/stepBy)
if steps > 0.0: #segment can be split
centPnt = pntCenter(pnt1,pnt2)
angle = p1p2Angle(pnt1, pnt2)
pairPnts = transRotatePnts([pnt1,pnt2],centPnt, angle)
outPnts = newPnts(pairPnts, stepBy)
angle = angle * (-1)
rotPnts = transRotatePnts(outPnts, [0,0], angle)
for pnt in rotPnts:
XY = [pnt[0] + centPnt[0], pnt[1] + centPnt[1]]
densified.append(XY)
#if (pnts[len(pnts)-1]) != None:
# densified.append(pnts[len(pnts)-1])
densified.append(pnts[len(pnts)-1])
return densified
def newPnts(pairPnts, stepBy):
#pairPnts a pair of points
#stepBy increment used to densify
#x1 = float(pairPnts[0][0]); y1 = float(pairPnts[0][1])
#x2 = float(pairPnts[1][0]); y2 = float(pairPnts[1][1])
x1 = pairPnts[0][0]; y1 = pairPnts[0][1]
x2 = pairPnts[1][0]; y2 = pairPnts[1][1]
dx = x2 - x1; dy = y2 - y1
pntsList = []
if dx == 0.0:
dist = abs(dy); slope = math.radians(90.0)
elif dy == 0.0:
dist = abs(dx); slope = 0.0
else:
dist = math.hypot(dx, dy);
slope = math.atan2(dy, dx) #atan2, returns angles in correct quad
inter = dy - slope*dx
i = stepBy
while i <= dist:
if dx == 0.0:
xNew = x1
if dy > 0:
yNew = y1 + i
else:
yNew = y1 - i
else:
xNew = x1 + math.cos(slope)*i #adjacent
yNew = y1 + math.sin(slope)*i
pntsList.append([xNew, yNew])
i = i + stepBy
return pntsList
def createPointFile(outFC, outType, SR, pnts, fieldsToAdd, arcpy):
#Requires
# outFC output feature class name "c:/temp/myfile.shp"
# outType output type "Polyline"
# SR Spatial reference from input feature class
# pnts a list of points (x,y,ID)
# eg [ [10,10,1], [12,15,2], [7,3,3] ]
# gp the Geoprocessor Object
#Create the output filename and feature class, use the Spatial Reference
# of the input feature class
#import os, sys #required if used in other modules
fullName = os.path.split(outFC)
outFolder = fullName[0].replace("\\","/")
outFName = fullName[1].replace(" ", "_")
outFullName = outFolder + "/" + outFName
try:
arcpy.CreateFeatureclass_management(outFolder, outFName, outType, "#", "Disabled", "Disabled", SR)
arcpy.AddMessage("\n" + "Creating output " + str(outType) + " file " + \
"\n" + "\t" + str(outFC)+ "\n")
except:
arcpy.AddMessage("Failed to create feature class" + arcpy.GetMessages())
sys.exit()
try:
arcpy.AddField_management(outFullName, "OldID", "LONG", "9")
except:
arcpy.AddMessage("Failed to add fields: " + arcpy.GetMessages())
#sys.exit()
#
#Create the insert cursor
aPnt = arcpy.CreateObject('Point')
try:
cur = arcpy.InsertCursor(outFC)
except:
arcpy.AddMessage("Failed to create cursor")
sys.exit()
anID=0
aCount = 0
for i in range(0,len(pnts)):
aList = pnts[i][0] #list of points
anID = pnts[i][1] #original ID value
for aPair in aList:
if aPair != None:
try:
aPnt.X = aPair[0]
aPnt.Y = aPair[1]
aPnt.ID = aCount
feat = cur.newRow()
feat.shape = aPnt
feat.setValue("ID", aCount)
feat.setValue("OldID", anID)
cur.insertRow(feat)
aCount = aCount + 1
except:
arcpy.AddMessage("Cannot create feature" + arcpy.GetMessages())
del cur
#------------------------------------------------------------------------------
def polyFromPoints (outFC, outType, SR, outPntsList, theFields, arcpy):
'''
Requires
outFC output feature class name "c:/temp/myfile.shp"
outType output type "Polygon" or "Polyline"
SR Spatial reference from input feature class
outPntsList a list of lists of points
eg [ [10,10], [12,15], [7,3] ]
theFields the fields to add and their format
gp the Geoprocessor Object
Create the output filename and feature class, use the Spatial Reference
of the input feature class
'''
#import os, sys #required if used in other modules
fullName = os.path.split(outFC)
outFolder = fullName[0]
outFName = fullName[1].replace(" ", "_")
outFullName = (outFolder + "/" + outFName).replace("\\","/")
#Create a polygon shapefile
try:
arcpy.CreateFeatureclass_management(outFolder, outFName, outType, "#", "Disabled", "Disabled", SR)
arcpy.AddMessage("\n" + "Creating output " + str(outType) + " file " + \
"\n" + "\t" + str(outFC))
except:
arcpy.AddMessage("Failed to create feature class")
sys.exit()
#
#Add any required fields
#
if len(theFields) != 0:
for i in range(0, len(theFields)):
aField = theFields[i]
try:
arcpy.AddField_management(outFullName, aField[0], aField[1], aField[2], aField[3])
except:
arcpy.AddMessage("Cannot add field " + str(aField[0]) + "\n" + arcpy.GetMessages())
try:
cur = arcpy.InsertCursor(outFC)
except:
arcpy.AddMessage("Failed to create cursor")
sys.exit()
anID=0
aPnt = arcpy.CreateObject("Point")
partArray = arcpy.CreateObject("Array")
polyArray = arcpy.CreateObject("Array")
for i in range(0,len(outPntsList)):
outShape = outPntsList[i][0] #list of points
anID = outPntsList[i][1] #original ID value
for aPair in outShape:
if aPair != None: #check for rings and parts
aPnt.X = aPair[0]
aPnt.Y = aPair[1]
partArray.add(aPnt)
else:
polyArray.add(partArray)
partArray.removeAll()
feat = cur.newRow()
feat.shape = polyArray
feat.ID = anID
if len(theFields) != 0:
for j in range(0, len(theFields)):
aField = theFields[j]
feat.setValue(aField[0],outPntsList[i][j+1])
cur.insertRow(feat)
polyArray.removeAll()
anID = anID+1
del polyArray
del aPnt
No comments:
Post a Comment