Wie entferne ich den Ozean mit OpenStreetMap? ⇐ Python

[Guest]

Ich arbeite mit den Daten der OpenStreetMap-Bibliothek (osmnx, in Python), um die Grenze von Salvador, Brasilien, zu extrahieren. Der Code entfernt erfolgreich Gewässer und erhält das Festland und die Inseln (ich lasse sie der Einfachheit halber weg), aber es gibt ein Problem mit einer künstlichen dreieckigen Grenze, die bis in den Ozean reicht.
Dieses Dreieck entsteht durch einen Teil der Verwaltungsgrenze, der aus irgendeinem Grund mit der natürlichen Küstenlinie verschmilzt und nicht vom Meer subtrahiert wird. Wie würden Sie versuchen, diese künstliche dreieckige Grenze in OpenStreetMap zu entfernen? Gibt es einen anderen Filter, den ich ausprobieren sollte? Ich möchte eine allgemeinere Lösung finden, als hierfür manuell „Debuff, Buff“ durchzuführen. Ich habe versucht, die Küste zu nutzen, aber die Ergebnisse waren nicht großartig.

Code: Select all

import osmnx as ox
import shapely

def get_city_boundary(place_name):
city_gdf = ox.geocode_to_gdf(place_name)
city_polygon = city_gdf.loc[0, "geometry"]

expand = 0
minx, miny, maxx, maxy = city_polygon.bounds
bounding_box = shapely.geometry.box(minx - expand, miny - expand, maxx + expand, maxy + expand)

try:
water_tags = {"natural": ["bay"], "maritime": "yes"}
water_gdf = ox.features_from_polygon(bounding_box, tags=water_tags)
if not water_gdf.empty:
water_union = water_gdf.geometry.union_all()
city_polygon = city_polygon.difference(water_union)
except Exception as e:
print(f"Error processing water areas: {e}")

city_polygon = city_polygon.buffer(0)
### change around here a bit if you want the islands of the city to appear
if isinstance(city_polygon, shapely.geometry.MultiPolygon):
polygons = list(city_polygon.geoms)
city_polygon = max(polygons, key=lambda p: p.area)

return city_polygon

place_name = "Salvador, Bahia, Brazil"
processed_polygon = get_city_boundary(place_name)

processed_gdf = ox.geocode_to_gdf(place_name).copy()
processed_gdf.loc[0, "geometry"] = processed_polygon

if processed_gdf.crs is None:
processed_gdf.set_crs(epsg=4326, inplace=True)

output_file = "salvador_processed_boundary.shp"
processed_gdf.to_file(output_file)