# Interpolate missing z values in polyline

3775
12
08-25-2016 05:27 AM
New Contributor III

I have a polylinezm imported from data put together from various sources. One issue I have is that some parts of the line are missing the Z value for the vertices:

To the left, a profile graph showing the Z values along the line, to the right sketch properties showing some vertices with a dummy -9999 Z value.

I've been looking for a way to interpolate the values for the missing vertices, but have not managed to find a way. Anyone got a tip on how to do this?

Tags (3)
1 Solution

Accepted Solutions
Esri Esteemed Contributor

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

# get polyline (first feature)
print "get polyline (first feature)"
polyline = arcpy.da.SearchCursor(fc, ('SHAPE@')).next()[0]

# first loop to get the values
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

# second loop to correct values
print "second loop to correct values"
print "points without valid z:", len(lst_nones)
for pnt_id in lst_nones:
try:
# extract previous, current and next point info
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]
##        print " - pnt_info_prev:", pnt_info_prev
##        print " - pnt_info_curr:", pnt_info_curr
##        print " - pnt_info_next:", pnt_info_next

# determine relative m position
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)
##        print "   - m:", m_prev, m_curr, m_next, m_pos

# get z value for current point
z_prev = pnt_info_prev[3]
z_next = pnt_info_next[3]
z_curr = m_pos * (z_next - z_prev) + z_prev
##        print "   - z:", z_prev, z_next, z_curr

# update lst
pnt_info_curr[3] = z_curr
##        print " - pnt_info_curr (2):", pnt_info_curr
lst_xyz[pnt_id] = pnt_info_curr

except:
# in case of error the invalid z will remain
pass

# create the list of corrected points
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]))

# generate polyline zm
print "generate polyline zm"
polyline_z = arcpy.Polyline(arcpy.Array(lst_pnts), sr, True, True)

# write to output shp
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.

12 Replies
MVP Legendary Contributor

your profile graph suggests a simple average of the two surrounding points, however, your table has removed this sequencing from easy solutions.  Are the points with the missing Z values actually grouped in one big lump?

New Contributor III

The missing z are several small lumps, og varying size. In some cases it's just one vertex, in other cases it's more...

Esri Esteemed Contributor

It might be easier to go back to the source (the raster used to interpolate the values of the polyline), correct the raster by interpolating the missing values and interpolate shape using the new surface.

The alternative is coding it in python. To get you started, in this post you will find some basic information:

New Contributor III

Well, this is more or less the source. It's not based on interpolation from a raster, but vector data from photogrammetric mapping where the troublesome parts are collected by digitizing from other sources (and hence no z) ...

I have some minor experience with handling z/m in arcpy before, but wanted to see if there are any easier way to do this that trying to write up a tool myself. Will have to look more into it, I guess...

Esri Esteemed Contributor

To bad, that this is the source. If you are planning to use Python, you could mimic an M value as 2D distance along the line for each point (first run). In the next you you would take the known Z values, like in your example:

ID 19292 has height 846.870
ID 19308 has height 838.360

And perform a linear "interpolation" to deduce the intermediate Z values based on the relative M position between the known Z value M positions.

Esri Esteemed Contributor

Could you post a single line. If I have some time later today, I can have a look.

New Contributor III

What I have is joined to one line, abt. 32k points, if that's manageable I can zip a shape-file of it...

Esri Esteemed Contributor

Sure no problem, just ZIP the shapefile and attach to the thread (you can use the link "use advanced editor" (upper right corner of your reply) and the attach the ZIP file using the "Attach" link in the lower left corner of you replay in the advanced editor.

Esri Esteemed Contributor

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

# get polyline (first feature)
print "get polyline (first feature)"
polyline = arcpy.da.SearchCursor(fc, ('SHAPE@')).next()[0]

# first loop to get the values
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

# second loop to correct values
print "second loop to correct values"
print "points without valid z:", len(lst_nones)
for pnt_id in lst_nones:
try:
# extract previous, current and next point info
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]
##        print " - pnt_info_prev:", pnt_info_prev
##        print " - pnt_info_curr:", pnt_info_curr
##        print " - pnt_info_next:", pnt_info_next

# determine relative m position
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)
##        print "   - m:", m_prev, m_curr, m_next, m_pos

# get z value for current point
z_prev = pnt_info_prev[3]
z_next = pnt_info_next[3]
z_curr = m_pos * (z_next - z_prev) + z_prev
##        print "   - z:", z_prev, z_next, z_curr

# update lst
pnt_info_curr[3] = z_curr
##        print " - pnt_info_curr (2):", pnt_info_curr
lst_xyz[pnt_id] = pnt_info_curr

except:
# in case of error the invalid z will remain
pass

# create the list of corrected points
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]))

# generate polyline zm
print "generate polyline zm"
polyline_z = arcpy.Polyline(arcpy.Array(lst_pnts), sr, True, True)

# write to output shp
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.