Wie entferne ich den Ozean mit OpenStreetMap?

Post a reply

Smilies
:) :( :oops: :chelo: :roll: :wink: :muza: :sorry: :angel: :read: *x) :clever:
View more smilies

BBCode is ON
[img] is ON
[flash] is OFF
[url] is ON
Smilies are ON

Topic review
   

Expand view Topic review: Wie entferne ich den Ozean mit OpenStreetMap?

by Guest » 03 Jan 2025, 19:16

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)
Image

Top