Skip to main content
2 of 2
added 5347 characters in body
Bera
  • 77.8k
  • 14
  • 78
  • 188

Create a point every meter along the line from each end to the start, and buffer each point starting at one meter and finish at the middle points buffer width.

Then union all the buffers.

import numpy as np
from shapely import wkb
from shapely.ops import unary_union
import geopandas as gpd
import matplotlib.pyplot as plt

riverwkb = b'\x01\x02\x00\x00\x00H\x00\x00\x00\x10\xdf\xd0\xdf\xd99"A\xe8@%\x97!U\\Ac\x8aM\xba\xd89"A\xc4iN\xbe$U\\A\x18\xddZP\xdd9"Au\xbd\x16\x9c\'U\\A/\xd7\xf8\xa1\xe79"A\x91\xa6.U*U\\A\xec\'\x90\xa8\xef9"A\xd8d\xa7W-U\\A\x0b \xb8j\xfd9"A^\xe3\x0e\xec/U\\A2\x16j\x9d\x0e:"A""e\x122U\\AY\x0c\x1c\xd0\x1f:"AR\xf6\n\x144U\\A\xd3\xadJ\xdd/:"A\xac\x9f\x11_6U\\ALOy\xea?:"A\x07I\x18\xaa8U\\AdI\x17<J:"A\x8d\xc7\x7f>;U\\A\x19\x9c$\xd2N:"A>\x1bH\x1c>U\\A\x11\x9e\x9aaK:"A\xefn\x10\xfa@U\\A\xaf\xf6\t\xa6E:"A\xa0\xc2\xd8\xd7CU\\A\xfa\xa3\xfc\x0fA:"A\xe7\x80Q\xdaFU\\A\xf2\xa5r\x9f=:"A\x98\xd4\x19\xb8IU\\A\xf2\xa5r\x9f=:"AI(\xe2\x95LU\\A\x9f\xfa\xf5\xc4>:"A\xfa{\xaasOU\\AdI\x17<J:"A\x80\xfa\x11\x08RU\\A\x8a?\xc9n[:"AE9h.TU\\A\xb93\x05\x12p:"A\x1dcL\x9dUU\\A\xf0%\xcb%\x88:"A\nx\xbeTVU\\Ay\xc3\r\x14\x9f:"A\xf6\x8c0\x0cWU\\A\x02aP\x02\xb6:"AM7\xf2\x9eWU\\A8S\x16\x16\xce:"A\xa3\xe1\xb31XU\\AoE\xdc)\xe6:"A\xcf\xb6\x14{XU\\A\xa57\xa2=\xfe:"AQ67WYU\\A\x89~\xebv\x17;"A\xd3\xb5Y3ZU\\A\x12\x1c.e.;"AT5|\x0f[U\\A9\x12\xe0\x97?;"A\x19t\xd25]U\\AA\x10j\x08C;"A`2K8`U\\AA\x10j\x08C;"A=[t_cU\\A|\xc1H\x917;"AXD\x8c\x18fU\\AU\xcb\x96^&;"A\x87\x182\x1ahU\\A\'\xd7Z\xbb\x11;"A\xf6\xac\xc6\xadiU\\A\xf0\xe4\x94\xa7\xf9:"A\xf9\xab\x0bfkU\\A\xc2\xf0X\x04\xe5:"A\xd2\xd5\xef\xd4lU\\A\xe6\xa7\x99;\xcf:"A\x15\x95#\x1fnU\\A]\nWM\xb8:"A\x83)\xb8\xb2oU\\A.\x16\x1b\xaa\xa3:"A\xc6\xe8\xeb\xfcpU\\A\xff!\xdf\x06\x8f:"A4}\x80\x90rU\\A\xd0-\xa3cz:"A8|\xc5HtU\\A\xfc\xe2m\x0bh:"A\xd2\xe5\xba%vU\\A\xe5\xe8\xcf\xb9]:"A\xed\xce\xd2\xdexU\\A8\x94L\x94\\:"A\x9e"\x9b\xbc{U\\A8\x94L\x94\\:"APvc\x9a~U\\AO\x8e\xea\xe5f:"A\xd5\xf4\xca.\x81U\\A\xc11\x8f\x82s:"A[s2\xc3\x83U\\A\x95|\xc4\xda\x85:"A\x8aG\xd8\xc4\x85U\\A\xc4p\x00~\x9a:"A\xcd\x06\x0c\x0f\x87U\\AM\x0eCl\xb1:"AO\x86.\xeb\x87U\\A\xd6\xab\x85Z\xc8:"A\x10\xc6?Y\x88U\\A\r\x9eKn\xe0:"Agp\x01\xec\x88U\\A\x96;\x8e\\\xf7:"A\x92Eb5\x89U\\A\x1f\xd9\xd0J\x0e;"A\xe9\xef#\xc8\x89U\\A\xa8v\x139%;"A~Z\xd4\xec\x89U\\A1\x14V\'<;"A\xaa/56\x8aU\\A\xbb\xb1\x98\x15S;"A?\x9a\xe5Z\x8aU\\A\xf1\xa3^)k;"A\xd5\x04\x96\x7f\x8aU\\A\xd5\xea\xa7b\x84;"A\x96D\xa7\xed\x8aU\\A\x0b\xddmv\x9c;"AW\x84\xb8[\x8bU\\A\x94z\xb0d\xb3;"A\x82Y\x19\xa5\x8bU\\A\x1e\x18\xf3R\xca;"Aon\x8b\\\x8cU\\A\xa7\xb55A\xe1;"A[\x83\xfd\x13\x8dU\\A8Q\x02\xa0\xfb;"A\xdd\x02 \xf0\x8dU\\AnC\xc8\xb3\x13<"A\xc9\x17\x92\xa7\x8eU\\A\xf7\xe0\n\xa2*<"AK\x97\xb4\x83\x8fU\\A.\xd3\xd0\xb5B<"A\xf8\xeb7\xa9\x90U\\A\xb7p\x13\xa4Y<"AzkZ\x85\x91U\\A@\x0eV\x92p<"A\x92U-\x86\x92U\\Ao\x02\x925\x85<"Aj\x7f\x11\xf5\x93U\\A\xe8\xa3\xc0B\x95<"A/\xbeg\x1b\x96U\\A'
points = [b'\x01\x04\x00\x00\x00\x01\x00\x00\x00\x01\x01\x00\x00\x00\xf0\xe4\x94\xa7\xf9:"A\xf9\xab\x0bfkU\\A',
          b'\x01\x04\x00\x00\x00\x02\x00\x00\x00\x01\x01\x00\x00\x00\x9f\xfa\xf5\xc4>:"A\xfa{\xaasOU\\A\x01\x01\x00\x00\x008Q\x02\xa0\xfb;"A\xdd\x02 \xf0\x8dU\\A']

river = wkb.loads(riverwkb)
riverdf = gpd.GeoDataFrame(geometry=[river], crs="epsg:3006")

pointgeoms = [wkb.loads(p) for p in points]
pointdf = gpd.GeoDataFrame(data={"type":["start", "end"]}, geometry=pointgeoms, crs="epsg:3006")
#     type                                           geometry
# 0  start                MULTIPOINT (597372.827 7427501.594)
# 1    end  MULTIPOINT (597279.385 7427389.807, 597501.813...

#Multipoints to singlepoints
pointdf["points"] = pointdf.apply(lambda x: [pnt for pnt in x.geometry.geoms], axis=1)
pointdf = pointdf.explode(column="points", ignore_index=True).set_geometry("points").drop(columns="geometry").rename_geometry("geometry")


centerbuffer = pointdf.loc[pointdf["type"]=="start"].buffer(70)

#Calculate distance along the river for each of the three points
pointdf["distance_along_line"] = pointdf.apply(lambda x: river.line_locate_point(x.geometry), axis=1)
#     type                        geometry  distance_along_line
# 0  start  POINT (597372.827 7427501.594)           423.724419
# 1    end  POINT (597279.385 7427389.807)           203.194798
# 2    end  POINT (597501.813 7427639.752)           770.110404

d1 = abs(pointdf.distance_along_line.iloc[0]-pointdf.distance_along_line.iloc[1]) #Number of points between top point and middle
d2 = abs(pointdf.distance_along_line.iloc[0]-pointdf.distance_along_line.iloc[2]) #Number of points between top point and bottom

#Create a point every 1 m from top to center point
p1 = [river.line_interpolate_point(x) for x in np.arange(203, 423+1, 1)]
#Buffer them starting at 1 m up to 70 m
buffers1 = [pnt.buffer(dist) for pnt, dist in zip(p1, np.arange(start=1, stop=70, step=(70-1)/len(p1)))]

p2 = [river.line_interpolate_point(x) for x in np.arange(423, 770+1, 1)]
buffers2 = [pnt.buffer(dist) for pnt, dist in zip(p2, np.arange(start=70, stop=1, step=-(70-1)/len(p2)))]

uu = unary_union(buffers1+buffers2)
buffdf = gpd.GeoDataFrame(geometry=[uu], crs="epsg:3006")

#Plot
fig, ax = plt.subplots(figsize=(15, 15))

#buffdf.plot(ax=ax, color="green", alpha=0.2)
buffdf.boundary.plot(ax=ax, color="purple", linewidth=3, linestyle=(0, (5, 10)))
pointdf.plot(ax=ax, markersize=200, color="black")
centerbuffer.boundary.plot(ax=ax, color="black", linewidth=4)
riverdf.plot(ax=ax, color="blue")

plt.show()

enter image description here

Bera
  • 77.8k
  • 14
  • 78
  • 188