1

I'm using Shapely to interpolate points every 500m along a linestring. For each interpolated point I'd like to find the nearest original points (immediately before and after) in the linestring.

(I'd like to do this so I can go back to work out the bearing of the linestring at each interpolated point.)

The only way I can think of doing this is as below. Is there a better way?

import matplotlib.pyplot as plt
from shapely import wkt
from shapely.geometry import Point, LineString
from shapely.ops import snap

river_wkt = "LINESTRING Z (242385.8 793780.99 0, 242422.39 793747.8 0, 242475.35 793754.62 0, 242533.52 793818.67 0, 242553.39 793853.04 0, 242649.87 793863.5700000001 0, 242686.89 793899.5700000001 0, 242746.79 793862.9300000001 0, 242772.02 793915.87 0, 242794.64 793905.02 0, 242808.22 793929.6900000001 0, 242911.26 793966.6 0, 242955.88 794023.61 0, 242982.75 794027.24 0, 243030.93 794075.77 0, 243072.94 794093.25 0, 243086.54 794087.1 0, 243186.3 794110.14 0, 243200.66 794160.13 0, 243225.4 794177.72 0, 243297.58 794186.88 0, 243350.03 794153.73 0, 243408.16 794220.9 0, 243476.18 794223.24 0, 243481.16 794251.97 0, 243479.31 794268.84 0, 243539.32 794383.92 0, 243579.45 794422.34 0, 243634.53 794441.3300000001 0, 243672.78 794486.98 0, 243678.1 794522.7000000001 0, 243738.22 794541.54 0, 243789.61 794503.21 0, 243818.01 794518.98 0, 243824.57 794552.21 0, 243811.55 794600.37 0, 243831.56 794632.54 0, 243851.77 794636.88 0, 243875.8 794618.6800000001 0, 243914.97 794552.4300000001 0, 243959.58 794567.67 0, 244007.48 794606.98 0, 244082.58 794613.9500000001 0, 244151.77 794684.5600000001 0, 244264.71 794739.8 0, 244272.96 794767.9400000001 0, 244243.45 794807.62 0, 244246.26 794828.3100000001 0, 244259.04 794841.2000000001 0, 244304.59 794840.37 0, 244354.74 794861.65 0, 244400.94 794838.17 0, 244471.7 794828.98 0, 244567.55 794883.3200000001 0, 244581.59 794929.38 0, 244632.91 794967.4300000001 0, 244702.02 794994.35 0, 244725.11 795058.16 0, 244751.86 795066.3100000001 0, 244817.33 795058.52 0, 244846.51 795095.09 0, 244949.85 795107.22 0, 245068.35 795178.92 0)"
river_line = shapely.wkt.loads(river_wkt)

# Interpolate
river_interp_points = []
for i in range(0, int(river_line.length), 500):
    river_interp_points.append(river_line.interpolate(i))
river_interp_line = LineString(river_interp_points)

# Find point on river closest to interpolated point #3
closest_point = snap(Point(river_interp_line.coords[3]), river_line, 50)
river_points = [(x,y) for x, y in zip(river_line.coords.xy[0], river_line.coords.xy[1])]
closest_index = river_points.index((closest_point.x, closest_point.y))

# Find points before and after point closest to interpolated point
point_before = river_points[closest_index-1]
point_after = river_points[closest_index+1]

fig, ax = plt.subplots()
ax.plot(river_line.coords.xy[0], river_line.coords.xy[1])
ax.plot(river_interp_line.coords.xy[0], river_interp_line.coords.xy[1], "o")
ax.plot(closest_point.coords.xy[0], closest_point.coords.xy[1], "*")
ax.plot(point_before[0], point_before[1], "x")
ax.plot(point_after[0], point_after[1], "x")

enter image description here

1

Not a complete solution, but maybe construct a KDTree from the points in the LineString and then use query to find nearest neighbors to points within the interpolated line.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.