如何在OSMnx中移除萨尔瓦多市边界的人工海洋三角区域？

我完全理解你遇到的困扰——OSM的行政边界偶尔会和海域边界出现奇怪的重叠，导致减去水体后留下这种突兀的三角区域。与其用手动缓冲这种治标不治本的方法，咱们可以从优化标签覆盖、完善裁剪逻辑这两个方向入手，找到更通用的解决方案：

一、优化水体标签，覆盖更多海洋区域 #

你当前的水体标签只包含了bay和maritime=yes，这大概率是漏减了部分海域的原因。OSM里标记海洋/海域的标签不止这些，咱们可以扩展标签集合，确保所有相关的海洋区域都被纳入裁剪范围：

修改你的water_tags为：

water_tags = {
    "natural": ["bay", "sea", "ocean"],
    "maritime": "yes",
    "water": ["sea", "ocean"],
    "place": ["sea"]
}

这样能覆盖OSM中绝大多数海洋相关的地理特征，减少漏减的概率。

二、结合海岸线进行二次裁剪 #

如果单纯优化水体标签还解决不了问题，咱们可以引入海岸线来进一步修正边界。你之前尝试过海岸线但效果不好，可能是没有用对方法——咱们可以先获取海岸线，再用海岸线围成的陆地范围和城市边界做交集，确保最终的边界只保留陆地区域：

在你的函数中添加海岸线处理的逻辑：

# 新增：获取海岸线并构建陆地裁剪范围
try:
    coastline_tags = {"natural": "coastline"}
    coastline_gdf = ox.features_from_polygon(bounding_box, tags=coastline_tags)
    if not coastline_gdf.empty:
        # 将海岸线转为线集合，然后生成陆地多边形
        coastline_lines = coastline_gdf.geometry.union_all()
        # 用shapely的polygonize生成由海岸线围成的多边形
        land_polygons = list(shapely.ops.polygonize(coastline_lines))
        # 找到和城市边界重叠的陆地多边形并合并
        land_union = shapely.ops.unary_union([poly for poly in land_polygons if poly.intersects(city_polygon)])
        city_polygon = city_polygon.intersection(land_union)
except Exception as e:
    print(f"Error processing coastline: {e}")

这个逻辑的核心是：用海岸线把陆地范围圈出来，再让城市边界只保留和陆地重叠的部分，从根源上排除海洋区域。

三、修正边界的拓扑问题 #

有时候这种三角区域的出现，也可能是因为原始行政边界存在拓扑错误（比如自相交）。你已经用了buffer(0)来修正，这很好，还可以加上shapely.make_valid()来进一步确保几何图形的有效性：

city_polygon = shapely.make_valid(city_polygon)
city_polygon = city_polygon.buffer(0)

整合后的完整代码 #

把这些优化点整合到你的函数中，修改后的代码如下：

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", "sea", "ocean"],
            "maritime": "yes",
            "water": ["sea", "ocean"],
            "place": ["sea"]
        }
        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}")

    # 第二步：用海岸线二次裁剪
    try:
        coastline_tags = {"natural": "coastline"}
        coastline_gdf = ox.features_from_polygon(bounding_box, tags=coastline_tags)
        if not coastline_gdf.empty:
            coastline_lines = coastline_gdf.geometry.union_all()
            land_polygons = list(shapely.ops.polygonize(coastline_lines))
            land_union = shapely.ops.unary_union([poly for poly in land_polygons if poly.intersects(city_polygon)])
            city_polygon = city_polygon.intersection(land_union)
    except Exception as e:
        print(f"Error processing coastline: {e}")

    # 修正拓扑错误
    city_polygon = shapely.make_valid(city_polygon)
    city_polygon = city_polygon.buffer(0)

    # 保留最大的多边形（去掉岛屿）
    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)

额外提示 #

如果以上方法还是解决不了，可能是OSM中萨尔瓦多的行政边界本身存在数据问题。这时候你可以去OpenStreetMap官网手动查看该区域的边界数据，确认是否有错误的节点/线段，必要时可以提交修正（当然这是更进阶的操作了）。

备注：内容来源于stack exchange，提问作者enriicoo