Polars处理NYC Taxi Zones数据构建MultiPolygon遇面积计算错误求助

问题背景 #

从SODA API获取NYC Taxi Zones的JSON数据集（非GeoJSON/Shapefile），用Polars的read_json读取后通过unnest展开the_geom字段。尝试转换coordinates列为Shapely可识别的MultiPolygon对象时，纽瓦克机场的多面数据能正常计算面积，但牙买加湾的数据触发「IllegalArgumentException: Points of LinearRing do not form a closed linestring」错误——原因是牙买加湾的坐标包含多个子列表，且存在未闭合的线性环，之前的展平操作破坏了几何结构的层级关系。

解决方案 #

核心思路是保留坐标的层级结构，先修复未闭合的LinearRing，再正确构造MultiPolygon对象，同时用Shapely的有效性校验处理潜在的几何问题。

步骤1：下载并读取数据

import polars as pl
import requests
from shapely.geometry import MultiPolygon, Polygon
from shapely.validation import make_valid

# 从SODA API自动下载数据
soda_url = "https://data.cityofnewyork.us/resource/i8iw-xf4u.json"
response = requests.get(soda_url)
raw_data = response.json()

# Polars读取并展开the_geom结构体
df = pl.read_json(raw_data).unnest("the_geom")

步骤2：编写坐标处理函数

这个函数会处理坐标层级，修复未闭合的环，并构造有效的MultiPolygon：

def process_coordinates(coords):
    polygons = []
    # 判断是否为MultiPolygon结构（外层包含多个Polygon）
    is_multi_poly = len(coords) > 0 and isinstance(coords[0][0][0], list)
    
    # 遍历每个Polygon的坐标
    target_coords = coords if is_multi_poly else [coords]
    for poly_coords in target_coords:
        # 处理外环
        exterior_ring = poly_coords[0]
        # 检查并闭合外环（首尾坐标不一致则补全）
        if exterior_ring[0] != exterior_ring[-1]:
            exterior_ring.append(exterior_ring[0])
        
        # 处理内环（如果有）
        interior_rings = []
        for ring in poly_coords[1:]:
            if ring[0] != ring[-1]:
                ring.append(ring[0])
            interior_rings.append(ring)
        
        # 构造单个Polygon并加入列表
        polygons.append(Polygon(exterior_ring, interior_rings))
    
    # 构造MultiPolygon并确保几何有效（处理自交等问题）
    multi_polygon = MultiPolygon(polygons)
    valid_geometry = make_valid(multi_polygon)
    return valid_geometry

步骤3：在Polars中应用函数计算面积

# 添加WKT格式几何列和计算面积列
df = df.with_columns(
    # 生成WKT字符串（可选）
    pl.col("coordinates")
      .map_elements(lambda x: process_coordinates(x).wkt, return_dtype=str)
      .alias("valid_wkt"),
    # 计算修复后的面积
    pl.col("coordinates")
      .map_elements(lambda x: process_coordinates(x).area, return_dtype=pl.Float64)
      .alias("calculated_shape_area")
)

# 查看牙买加湾的数据（示例）
print(df.filter(pl.col("zone") == "Jamaica Bay").select("zone", "shape_area", "calculated_shape_area"))

关键说明 #

1. 层级保留：不再直接展平坐标，而是根据MultiPolygon的层级结构（[Polygon1, Polygon2...]，每个Polygon包含[外环, 内环1, 内环2...]）逐个处理，避免破坏几何结构。
2. 环闭合处理：每个线性环（外环/内环）都会检查首尾坐标，不一致则自动补充尾坐标闭合环，解决报错根源。
3. 有效性校验：使用make_valid函数处理几何对象可能存在的自交、重叠等问题，确保后续计算正常。

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