OK, so here's the result (in green as 3D view):
The code I used is:
def main():
import arcpy
fc = 'D:\Xander\GeoNet\InterpolatePolylineZ\shp\BRGB_Ny_20160825.shp'
fc_out = 'D:\Xander\GeoNet\InterpolatePolylineZ\shp\BRGB_Ny_20160825_corr.shp'
sr = arcpy.Describe(fc).spatialReference
print "get polyline (first feature)"
polyline = arcpy.da.SearchCursor(fc, ('SHAPE@')).next()[0]
print "first loop to get the values"
cnt = 0
lst_xyz = []
lst_nones = []
lst_oks = []
pnt_prev = None
for part in polyline:
for pnt in part:
x = pnt.X
y = pnt.Y
z = pnt.Z
if z < -9998:
z = None
lst_nones.append(cnt)
else:
lst_oks.append(cnt)
if pnt_prev == None:
dist = 0
m = dist
else:
m += getDist(pnt, pnt_prev)
lst_xyz.append([cnt, x, y, z, m])
pnt_prev = pnt
cnt += 1
if cnt % 10000 == 0:
print " - processing point:", cnt
print "second loop to correct values"
print "points without valid z:", len(lst_nones)
for pnt_id in lst_nones:
try:
pnt_info_curr = lst_xyz[pnt_id]
pnt_info_prev = lst_xyz[pnt_id-1]
pnt_id_next = min([p for p in lst_oks if p > pnt_id])
pnt_info_next = lst_xyz[pnt_id_next]
m_curr = pnt_info_curr[4]
m_prev = pnt_info_prev[4]
m_next = pnt_info_next[4]
m_pos = (m_curr - m_prev) / (m_next - m_prev)
z_prev = pnt_info_prev[3]
z_next = pnt_info_next[3]
z_curr = m_pos * (z_next - z_prev) + z_prev
pnt_info_curr[3] = z_curr
lst_xyz[pnt_id] = pnt_info_curr
except:
pass
print "create the list of corrected points"
lst_pnts = []
for pnt_info in lst_xyz:
lst_pnts.append(arcpy.Point(pnt_info[1], pnt_info[2], pnt_info[3], pnt_info[4]))
print "generate polyline zm"
polyline_z = arcpy.Polyline(arcpy.Array(lst_pnts), sr, True, True)
print "write to output shp"
arcpy.CopyFeatures_management([polyline_z], fc_out)
def getDist(pnt1, pnt2):
import math
return math.hypot(pnt1.X-pnt2.X, pnt1.Y-pnt2.Y)
if __name__ == '__main__':
main()
You should implement an update cursor to update the original geometry.
I have also attached the output shapefile (no attributes) for you te verify the result.