Interpolate missing z values in polyline

1917
12
Jump to solution
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?

Reply
0 Kudos
1 Solution

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

View solution in original post

12 Replies
Highlighted
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?

Reply
0 Kudos
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...

Reply
0 Kudos
Highlighted
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:

https://community.esri.com/people/xander_bakker/blog/2016/07/08/working-with-3d-and-m-aware-geometri... 

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

Reply
0 Kudos
Highlighted
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
Esri Esteemed Contributor

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

Reply
0 Kudos
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...

Reply
0 Kudos
Highlighted
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.

Reply
0 Kudos
Highlighted
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.

View solution in original post