Python如何按米级绝对数值缩放GeoJSON多边形

我需要编写Python代码，实现输入包含多边形要素的GeoJSON数据后，输出以米为单位、按固定绝对数值向外扩展的多边形结果。
举个实际场景：如果输入是代表房屋的正方形多边形，我需要得到每一条边都恰好向外扩大20米的新正方形。这类逻辑对规则正方形实现简单，但处理复杂不规则多边形时难度较高。

我当前参考已有方案实现的代码如下，目前采用的是相对比例缩放逻辑（默认放大10%），无法满足固定米数的绝对缩放需求：

import json
from shapely import affinity
from shapely.geometry import shape, Point, mapping

def handler(event, context):
    with open('polygon.geojson', encoding='utf-8') as f:
        js = json.load(f)
    polygon = shape(js['features'][0]['geometry'])
    polygon_nearby = affinity.scale(polygon, xfact=1.1, yfact=1.1)
    print(json.dumps({"type": "FeatureCollection", "features": [{"type": "Feature", 'properties': {}, 'geometry': mapping(polygon)}]}))
    print(json.dumps({"type": "FeatureCollection", "features": [{"type": "Feature", 'properties': {}, 'geometry': mapping(polygon_nearby)}]}))

if __name__ == '__main__':
    handler(None, None)

我需要把上述逻辑修改为支持按米为单位设置绝对扩展距离的版本。

你当前用的shapely.affinity.scale方法仅对坐标值做比例乘算，完全不考虑坐标对应的实际地理距离，既无法保证固定米数的扩展量，在不规则多边形上也会出现各边外扩距离不一致的问题，不适用这个场景。
你需要的“多边形各边向外扩展固定米数”本质是地理要素的正缓冲区操作，按以下步骤实现即可得到精确结果：

• 原始GeoJSON一般为WGS84（EPSG:4326）经纬度坐标，单位是度，不能直接做米制距离计算，需要先转换为单位为米的投影坐标系
• 对投影后的米制多边形调用缓冲区生成方法，传入需要扩展的米数作为参数，该方法会自动处理所有复杂多边形的边角，保证新多边形边界与原多边形边界的垂直距离恰好为设定值
• 将生成的缓冲区多边形重新转回WGS84经纬度坐标，即可导出为符合要求的GeoJSON

首先安装需要的依赖：

pip install shapely pyproj

可直接使用的修正后代码如下，默认设置为向外扩展20米，可自行修改expand_distance参数调整扩展距离：

import json
from shapely.geometry import shape, mapping
from shapely.ops import transform
from pyproj import Transformer

def expand_polygon_by_meter(input_geojson_path, output_geojson_path, expand_distance=20):
    # 读取原始GeoJSON
    with open(input_geojson_path, encoding='utf-8') as f:
        js = json.load(f)
    
    polygon = shape(js['features'][0]['geometry'])
    # 计算多边形中心经纬度，自动匹配对应UTM投影带（米制单位，精度满足绝大多数场景）
    center_lon = polygon.centroid.x
    center_lat = polygon.centroid.y
    utm_zone = int((center_lon + 180) // 6) + 1
    utm_epsg = 32600 + utm_zone if center_lat >=0 else 32700 + utm_zone

    # 定义坐标转换器：WGS84转UTM，UTM转回WGS84
    to_meter_transformer = Transformer.from_crs("EPSG:4326", f"EPSG:{utm_epsg}", always_xy=True).transform
    to_wgs84_transformer = Transformer.from_crs(f"EPSG:{utm_epsg}", "EPSG:4326", always_xy=True).transform

    # 转米制坐标、生成外扩缓冲区、转回经纬度
    polygon_meter = transform(to_meter_transformer, polygon)
    # join_style=2 表示边角用斜接方式，处理矩形等直角要素时会保持直角，不会生成圆角
    expanded_meter = polygon_meter.buffer(expand_distance, join_style=2)
    expanded_wgs84 = transform(to_wgs84_transformer, expanded_meter)

    # 导出结果
    result = {
        "type": "FeatureCollection",
        "features": [
            {"type": "Feature", "properties": {}, "geometry": mapping(polygon)},
            {"type": "Feature", "properties": {}, "geometry": mapping(expanded_wgs84)}
        ]
    }
    with open(output_geojson_path, 'w', encoding='utf-8') as f:
        json.dump(result, f, ensure_ascii=False)
    return result

if __name__ == '__main__':
    # 参数依次为：输入文件路径、输出文件路径、外扩米数
    expand_polygon_by_meter('polygon.geojson', 'expanded_polygon.geojson', expand_distance=20)

注意事项 #

• 上述代码自动根据多边形中心匹配UTM投影带，在南北纬84度之间的区域，距离误差小于1米，满足房屋、地块这类要素的精度要求
• 如果处理的多边形跨多个UTM带，可替换为当地适用的本地米制投影坐标系，精度会更高
• 如果需要外扩后的多边形边角为圆角，去掉buffer方法中的join_style=2参数即可

内容的提问来源于stack exchange，提问作者horin