# Interpolate missing z values in polyline

1917
12
08-25-2016 05:27 AM
Highlighted
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
Highlighted
by 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()    # 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            m_prev = pnt_info_prev            m_next = pnt_info_next            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            z_next = pnt_info_next            z_curr = m_pos * (z_next - z_prev) + z_prev    ##        print "   - z:", z_prev, z_next, z_curr            # update lst            pnt_info_curr = 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, pnt_info, pnt_info, pnt_info))    # 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
Highlighted
by MVP Esteemed 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?

Highlighted
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...

Highlighted
by 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:

Highlighted
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...

Highlighted
by 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.

Highlighted
by Esri Esteemed Contributor

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

Highlighted
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...

Highlighted
by 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.

Highlighted
by 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()    # 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            m_prev = pnt_info_prev            m_next = pnt_info_next            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            z_next = pnt_info_next            z_curr = m_pos * (z_next - z_prev) + z_prev    ##        print "   - z:", z_prev, z_next, z_curr            # update lst            pnt_info_curr = 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, pnt_info, pnt_info, pnt_info))    # 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. 