3

I'm currently working on a project involving creating an isochrone map of historical water travel. At the moment, I have a geodataframe containing a starting point on the river, a multipoint of the upstream and downstream endpoints that can be reached within the given timeframe, and a buffer around the starting point representing the distance that can be walked from the starting point in the given time. I also have a linestring of the river. I'm trying to create a polygon representing the total area that can be covered through a combination of boats and walking, so I basically need to create a buffer around the river that is widest at the starting point (the width of my current circular buffer) and tapers towards the ends. In the picture below, I have the river in blue and the points shown in black, and I need to generate a shape similar to the purple dotted outline. Drawing of a river with start and endpoints, as well as desired buffer

From previous research, I know I could create a buffer around the river, but only of uniform thickness. I could also generate a convex hull, but that would not conform well to the shape of the river. To make matters worse, I need to do this for about 400 separate locations, so I need a solution that can be coded. I'm fairly new to this, so apologies for anything I've overlooked!

Vince
  • 20,017
  • 15
  • 45
  • 64
bseal
  • 31
  • 2
  • See possible ideas here: https://gis.stackexchange.com/questions/326681/getting-geometries-of-hurricanes-cone-of-uncertainty-using-shapely, https://gis.stackexchange.com/questions/340968/varying-size-buffer-along-line-with-postgis – Mike T May 18 '23 at 10:51
  • JTS has this feature, but it has not been implemented in GEOS/Shapely/GeoPandas yet. – Mike T May 18 '23 at 10:54

1 Answers1

0

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.UU\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;"A2K8U\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
  • 72,339
  • 13
  • 72
  • 161