by
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?
1 Solution
Accepted Solutions
by
XanderBakker
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.
by
by
by
by
XanderBakker
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.
