使用 Polars 从结构类型和 list[list[list[f64]]] 列创建 MultiPolygon 对象

Question

我已经下载了 NYC Taxi Zones dataset（从 SODA Api 下载并保存为 json 文件 - 不是 GeoJson 或 Shapefile）。数据集相当小，因此我使用包含的全部信息。为了方便发帖，我展示了数据集的前 2 行：

• 原始（结构类型值为

  the_geom

  ）。
在极坐标中使用 unnest() 命令解压结构类型后的
• 数据集。 --使用

  write_ndjson()

  命令更新

原始数据集：

应用

unnest()

并选择一些列和前 2 行后，位于数据集下方。

您可以使用以下命令使用极坐标导入数据

import polars as pl
poc = pl.read_json("./data.json"))

我对多边形感兴趣。我实际上正在尝试使用 multipolygon 和 wkt（众所周知的文本表示） -

shape_area

模块使用的方法来重新计算

shapely

。

到目前为止我所做的是使用列

coordinates

并将其转换为 MultiPolygon() 对象 - 由 Shapely 模块读取。

def flatten_list_of_lists(data):
    return [subitem3 for subitem1 in data for subitem2 in subitem1 for subitem3 in subitem2]

该函数将

list[list[list[list[f64]]]]

对象作为输入，并转换为

list[list[f64]]

对象。

flattened_lists = [flatten_list_of_lists(row) for row in poc["coordinates"].to_list()]

print(flattened_lists)

[[-74.18445299999996, 40.694995999999904], [-74.18448899999999, 40.69509499999987], [-74.18449799999996, 40.695184999999 87]、[-74.18438099999997、40.69587799999989]、[-74.18428199999994、40.6962109999999]、[-74.18402099999997、40.69707499999989]...

然后我使用下面的函数来应用字符串连接，基本上：

• 将

  list[list[f64]]

  对象转换为字符串。
• 在字符串前面添加关键字 MultiPolygon。
• 分别用 '('., ')' 替换 '[' 和 ']'。

def polygon_to_wkt(polygon):
    # Convert each coordinate pair to a string formatted as "lon lat"
    coordinates_str = ", ".join(f"{lon} {lat}" for lon, lat in polygon)
    # Format into a WKT Multipolygon string (each polygon as a single polygon multipolygon)
    return f"""MULTIPOLYGON ((({coordinates_str})))"""

formatted_wkt = [polygon_to_wkt(polygon) for polygon in flattened_lists]
poc = poc.with_columns(pl.Series("WKT", formatted_wkt))

最后，我使用方法

wkt.loads("MultiPolygon ((()))").area

来计算Multipolygon对象的形状面积

def convert_to_shapely_area(wkt_str):
    try:
        return wkt.loads(wkt_str).area
    except Exception as e:
        print("Error converting to WKT:", e)
        return None
poc = poc.with_columns(
    pl.col("WKT").map_elements(convert_to_shapely_area).alias("shapely_geometry")
)
print(poc.head())

即使对于第一个形状，WKT 正确返回对象的面积，而对于第二个 MultiPolygon，该方法返回以下错误：

IllegalArgumentException：LinearRing 的点未形成闭合线串

我在两行之间注意到的是，纽瓦克机场的多边形是 list[list[f64]]] 坐标的连续对象。然而，牙买加湾有多个子列表 [list[list[f64]]] 元素（检查列坐标以验证这一点）。下面的截图也验证了这一说法。

因此，在应用 WKT 之前，有没有办法将牙买加湾的多边形统一为一个几何对象？

P.S：GitHub 上的许多解决方案都使用形状文件，但我想定期使用 SODA API 自动重新下载 NYC 区域。

从 SODA API 下载原始 .json 文件（省略

logger_object

让我们假装它是

print()

）

import requests
params = {
        "$limit": geospatial_batch_size, #i.e. 50_000
        "$$app_token": config.get("api-settings", "app_token")
    }
    response = requests.get(api_url, params=params)
    if response.status_code == 200:
        data = response.json()
        if not data:
            logger_object.info("No data found, please check the API connection.")
            sys.exit()
        with open("{0}/nyc_zone_districts_data.json".format(geospatial_data_storage), "w") as f:
            json.dump(data, f, indent=4)
    else:
        logger_object.error("API request failed.")
        logger_object.error("Error: {0}".format(response.status_code))
        logger_object.error(response.text)

Answer 1

我最终找到了解决这个问题的方法。正如已经描述的，我正在展平坐标列表（又名点）。地理空间数据中的点是 (x,y) 坐标。因此，多多边形是多个点的组合。

def flatten_list_of_lists(data) -> list:
    return [subitem2 for subitem1 in data for subitem2 in subitem1]

，上述函数很重要，因为用作输入参数的对象

data

的类型为

[list[list[list[f64]]]]

，并且

MultiPolygon

每个点都有特定的基数级别。

然后我使用

shapely

将其展平列表转换为 MultiPolygon 对象

from shapely.geometry import Polygon, MultiPolygon, Point
def transform_polygons_to_multipolygons(flat_list:list) -> list:
    return [ MultiPolygon( [Polygon(coord) for coord in polygon]).wkt for polygon in  flat_list]

您会注意到，创建每个

MultiPolygon

对象后，我将其保存为

WKT

格式（因此作为字符串）。

最后，我计算多边形的面积

from shapely import wkt
def compute_geo_areas(multipolygons:MultiPolygon) -> float:
    return wkt.loads(multipolygons).area

最终代码

flattened_lists:list = [flatten_list_of_lists(row) for row in df_geo["coordinates"].to_list()]
multipolygons:list = transform_polygons_to_multipolygons(flattened_lists)

df_geo = df_geo.with_columns(pl.Series("multipolygons", multipolygons)) \
.with_columns(polygon_area=pl.col("multipolygons").map_elements(compute_geo_areas)) \
.drop("coordinates")

，其中

df_geo

是从 SO 问题上附加的 JSON 文件加载的 Polars DataFrame。