Since I already started with some Python code I decided to give it a go and I came up with some crappy code that actually results in a single 3D line. So this only processes a single line and the points that are within a distance of 1m tolerance to that line. It works for a test set and also the data Erin attached to this thread.It works like this:
- it converts the vertices of the line to a list (X,Y,Z,M,"line"), Z = None and M is distance from start
- the points are also converted to a list (X,Y,Z,M,"point"), Z has value and M is not set yet
- the lists are merged together (M for points is calculated)
- the first and last Z values are determined and set to the from and to points
- the intermediate vertices are provided with a Z value based on the known values and M locations of points
- the list is converted to a 3D polyline featureclass
import arcpy,math
fcPnt = r'C:\Project\_Forums\Erin\Pipetest\points.shp'
fcLine = r'C:\Project\_Forums\Erin\Pipetest\pipe.shp'
fcLine3D = r'C:\Project\_Forums\Erin\test.gdb\pipe3D' # output
tolerance = 1
fldZpnt = "abs_height" # "SHAPE@Z"
# fetch polyline
with arcpy.da.SearchCursor(fcLine, ("SHAPE@")) as cursor:
for row in cursor:
polyline = row[0]
del row
# get sr from input lines
sr = arcpy.Describe(fcLine).spatialReference
# put polyline crds in list
lstPnts1 = []
arr = arcpy.Array()
dist = 0
pnt = arcpy.Point(0,0)
for arr in polyline:
for i in range(0,arr.count):
prevPnt = arcpy.Point(pnt.X,pnt.Y)
pnt = arr.next()
# Add distance from startpoint (M)
if i == 0:
dist = 0
else:
dist += math.hypot(pnt.X-prevPnt.X,pnt.Y-prevPnt.Y)
lstPnt = []
lstPnt.extend((pnt.X,pnt.Y,pnt.Z,dist,"line"))
lstPnts1.append(lstPnt)
##print "lstPnts1: {0}\n".format(lstPnts1)
lstPnts2 = []
cnt = 0
with arcpy.da.SearchCursor(fcPnt, ("SHAPE@X","SHAPE@Y",fldZpnt)) as cursor:
for row in cursor:
cnt += 1
lstPnt = []
lstPnt.extend((row[0],row[1],row[2],-1,"point"))
lstPnts2.append(lstPnt)
del row
##print "lstPnts2: {0}\n".format(lstPnts2)
# merge lists together
lstPnts3 = []
# loop through segments of polyline
for i in range(1,len(lstPnts1)):
X1,Y1,Z1,M1,geom1 = lstPnts1[i-1]
X2,Y2,Z2,M2,geom2 = lstPnts1
pnt1 = arcpy.Point(X1,Y1)
pnt2 = arcpy.Point(X2,Y2)
arrLine = arcpy.Array()
arrLine.add(pnt1)
arrLine.add(pnt2)
line = arcpy.Polyline(arrLine)
if i==1:
lstPnts3.append(lstPnts1[i-1])
for j in range(0,len(lstPnts2)):
Xn,Yn,Zn,Mn,geomn = lstPnts2
pntn = arcpy.Point(Xn,Yn)
dist = line.distanceTo(pntn)
if dist < tolerance:
# determine M value, assume point is "on" line
Mn = M1 + math.hypot(X1-Xn,Y1-Yn)
lstPnt= []
lstPnt.extend((Xn,Yn,Zn,Mn,geomn))
lstPnts3.append(lstPnt)
lstPnts3.append(lstPnts1)
##print ""
##print "X\tY\tZ\tM\tGeometry"
##for p in lstPnts3:
## print "{0}\t{1}\t{2}\t{3}\t{4}".format(p[0],p[1],p[2],p[3],p[4]).replace('.',',')
##print ""
# no extrapolation, just propagate last know Z
firstZ = None
lastZ = None
cnt=-1
for p in lstPnts3:
cnt+=1
if p[2] != None:
lastZ = p[2]
if firstZ == None:
firstZ = p[2]
##print ""
##print "first Z: {0}".format(firstZ)
##print "last Z: {0}".format(lastZ)
##print ""
# update firstZ and lastZ for from and to point
lstPnts3[0][2] = firstZ
lstPnts3[len(lstPnts3)-1][2] = lastZ
##print ""
##print "X\tY\tZ\tM\tGeometry (lastZ en firstZ)"
##for p in lstPnts3:
## print "{0}\t{1}\t{2}\t{3}\t{4}".format(p[0],p[1],p[2],p[3],p[4]).replace('.',',')
##print ""
# now fill up the blanks
print "Fill up the blanks"
cnt=-1
for p in lstPnts3:
cnt+=1
if p[2] == None:
currM = p[3]
# get known Z,M before current point
for i in range(cnt-1,0,-1):
if lstPnts3[2] != None:
beforeM = lstPnts3[3]
beforeZ = lstPnts3[2]
break
# get known Z,M after current point
for i in range(cnt+1,len(lstPnts3),1):
if lstPnts3[2] != None:
afterM = lstPnts3[3]
afterZ = lstPnts3[2]
break
if afterM - beforeM > 0:
currZ = beforeZ + (afterZ - beforeZ) * (currM - beforeM) / (afterM - beforeM)
else:
# in case of 2 points on same location...
currZ = beforeZ
lstPnts3[cnt][2] = currZ
##print ""
##print "X\tY\tZ\tM\tGeometry"
##for p in lstPnts3:
## print "{0}\t{1}\t{2}\t{3}\t{4}".format(p[0],p[1],p[2],p[3],p[4]).replace('.',',')
##print ""
# time to write the 3D line to a featureclass
features =[]
arr = arcpy.Array()
for p in lstPnts3:
if p[4] == "line":
pnt = arcpy.Point(p[0],p[1],p[2],p[3])
pntg = arcpy.PointGeometry(pnt,sr,True,True)
arr.add(pnt)
polyline3D = arcpy.Polyline(arr,sr,True,True)
features.append(polyline3D)
arcpy.CopyFeatures_management(features, fcLine3D)
print "ready..."Maybe this is useful to someone...Kind regards,Xander
