How to interpolate line Z from sparse attribute values

Question

I have a series of end-to-end connected 3d lines which I'd like to add elevation values for all vertices on the lines which are currently null or 0. How to accomplish this?

In the following image there are 3 lines, with 1 selected and vertices marked (green squares). Of these 3 lines there 2 vertices which have a Z value. Each Z-bearing vertex is on a different line. How to calculate and assign the intermediate vertex Z values between A and B? Again referring to the image, the resultant point A.15 would be 661, A.14 is 662, ...

Highlighted line with vertex coordinates showing only 1 populated Z

(For simplicity's sake this answer doesn't need to take X,Y distance between vertices into account. It's okay to simply divide A.Z & B.Z difference by number of points. Bonus kudos if that's included anyway ;-)

Both ArcGIS and QGIS/OGR solutions welcome; ditto ready to run tools and/or python.

• matt wilkie: interpolate Z between polyline's vertices (IZ.InterpolateZsBetween ?) looks relevant, if can figure out how to use from python, and span lines.
• Michael Stimson: I've done this in C# using IZ, it's not that difficult but there's a bit of a learning curve ahead if you don't know ArcObjects or either VB.net or C#. It can be done in python I think... if you've only got a few good Z values you can create a terrain from the points and use Interpolate Shape resources.arcgis.com/en/help/main/10.1/index.html#/… to get the values onto the line (3d Analyst required). Also a good post gis.stackexchange.com/questions/119249/… that could be modified
• matt wilkie: Thanks for the pointers @MichaelMiles-Stimson, but generating a surface to interpolate from isn't viable. The known-Z points are too sparse. I know how to get straight to ArcObjects from python, it's the syntax for using the objects themselves that stymies me. I guess it's time I cracked that nut though. ;-)
• Michael Stimson: There's a few good posts that utilize ArcObjects in python, see gis.stackexchange.com/questions/80/… for the basics and gis.stackexchange.com/questions/142734/… for an example (unrelated) and a 'how to' read the ArcObjects API documentation.

Answer 1

Very long but working:

# Tool assumes:
# a) 1st layer in Table of Conternt is TARGET polylineZ feature
# b) records in table sorted from upstream down 
# c) selected features maintain direction, i.e. upstream line END point equal downstream line START point
# d) missing z value for vertex is 0
import arcpy, traceback, os, sys
from arcpy import env
env.outputZFlag = "Enabled"
env.outputMFlag = "Enabled"
try:
 def showPyMessage():
 arcpy.AddMessage(str(time.ctime()) + " - " + message)
 mxd = arcpy.mapping.MapDocument("CURRENT")
 destLR = arcpy.mapping.ListLayers(mxd)[0]
 dDest=arcpy.Describe(destLR)
 destProj=dDest.spatialReference 
# copy all lines to manageable list and merge
 g=arcpy.Geometry()
 geometryList=arcpy.CopyFeatures_management(destLR,g)
 nLines=len(geometryList)
 allPoints=geometryList[0].getPart(0)
 for i in range (1,nLines):
 part=geometryList[i].getPart(0)
 allPoints.extend(part)
 longLine=arcpy.Polyline(allPoints)
# calc chainages
 dictFeatures,m = {},0
 for p in allPoints:
 dist=longLine.measureOnLine(p)
 dictFeatures[m]=(dist,p.Z)
 m+=1
# find pairs
 first,nPoints,pairs=-1,len(allPoints),[]
 for i in range(nPoints):
 (d,z)=dictFeatures[i]
 if abs (z-0.0)>0.001: first=i; break
 if first==-1 or first==(nPoints-1):
 arcpy.AddMessage("Not enough z info")
 sys.exit(0)
 pairs.append((first,d,z))
 while True:
 second =-1
 for i in range(first+1,nPoints):
 (d,z)=dictFeatures[i]
 if abs (z-0.0)>0.001: second=i; break
 if second==-1:break
 first=second
 pairs.append((second,d,z))
# interpolate
 n=len(pairs)
 for j in range(1,n):
 first,d1,z1=pairs[j-1]
 second,d2,z2=pairs[j]
 dz=(z2-z1)/(d2-d1)
 for i in range(first+1,second):
 d,z=dictFeatures[i]
 z=z1 + dz*(d-d1)
 dictFeatures[i]=(d,z)
# update 
 with arcpy.da.UpdateCursor(destLR,"Shape@") as cursor:
 aKey,m=0,0
 pz=arcpy.Point()
 newL=arcpy.Array()
 for row in cursor:
 shp=geometryList[m];m+=1
 part=shp.getPart(0)
 n=len(part)
 for i in range(n):
 p=part[i]
 d,z=dictFeatures[aKey]
 pz.X,pz.Y,pz.Z=p.X,p.Y,z
 newL.add(pz)
 aKey+=1
 newLine=arcpy.Polyline(newL,destProj,True)
 newL.removeAll() 
 arcpy.AddMessage(newLine.length)
 cursor.updateRow((newLine,))
except:
 message = "\n*** PYTHON ERRORS *** "; showPyMessage()
 message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
 message = "Python Error Info: " + str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage() 

• matt wilkie: Thank you very much. This provides a solid starting point, and is much better than the rabbit hole I was in, trying to adapt @Tomek's gis.stackexchange.com/a/18655/108 to my circumstance.
• FelixIP: you still have to figure out how to extrapolate points at the very start and end lines