topic Re: Interpolate missing z values in polyline in Data Management Questions

Interpolate missing z values in polyline

HåvardMoe1 — Thu, 25 Aug 2016 12:27:14 GMT

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?

Re: Interpolate missing z values in polyline

DanPatterson_Retired — Thu, 25 Aug 2016 12:43:59 GMT

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?

Re: Interpolate missing z values in polyline

XanderBakker — Thu, 25 Aug 2016 12:51:44 GMT

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-geometries-in-arcpy

Re: Interpolate missing z values in polyline

HåvardMoe1 — Thu, 25 Aug 2016 12:55:17 GMT

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

Re: Interpolate missing z values in polyline

HåvardMoe1 — Thu, 25 Aug 2016 13:04:40 GMT

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

Re: Interpolate missing z values in polyline

XanderBakker — Thu, 25 Aug 2016 13:12:16 GMT

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.

Re: Interpolate missing z values in polyline

XanderBakker — Thu, 25 Aug 2016 13:14:31 GMT

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

Re: Interpolate missing z values in polyline

HåvardMoe1 — Thu, 25 Aug 2016 13:37:58 GMT

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

Re: Interpolate missing z values in polyline

XanderBakker — Thu, 25 Aug 2016 13:54:08 GMT

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.

Re: Interpolate missing z values in polyline

XanderBakker — Sat, 11 Dec 2021 06:04:08 GMT

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.

Re: Interpolate missing z values in polyline

XanderBakker — Thu, 25 Aug 2016 16:35:34 GMT

I forgot to mention that since the featureclass is m-enabled, I stored the 2D distance as M too.

Before: After:

Re: Interpolate missing z values in polyline

HåvardMoe1 — Fri, 26 Aug 2016 07:30:00 GMT

Oh, that looks very nice! Reading through the code I see some nice solutions to a couple of things I was struggling to wrap my head around when I was thinking through this yesterday.

Re: Interpolate missing z values in polyline

XanderBakker — Fri, 26 Aug 2016 11:49:04 GMT

I'm glad you like it. If you have any problems implementing the Update Cursor, just let me know.