我已经下载了NYC Taxi Zones dataset(从SODA Api下载并保存为json文件-不是GeoJson或Shapetime).该数据集相当小,因此我使用了包含的全部信息.为了文章的方便,我展示了数据集的前2行:

  • original( struct 类型值为the_geom).
  • 在polars中使用unnest()命令解包 struct 类型后dataset.--更新了write_ndjson()命令

The original dataset: enter image description here

应用unnest()并 Select 一些列和前2行后,在数据集下方.

enter image description here

您可以使用以下命令使用Polars导入数据

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

我对多边形感兴趣.实际上,我正在try 使用shapely模块使用的多多边形和wkt(众所周知的文本表示)方法重新计算shape_area.

到目前为止,我所做的是使用列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.69499599999904], [-74.184488999999,40.695094999999987], [-74.18449799999996,40.6951849999987], [-74.184380999999997,40.695877999999989], [-74.1842819999994,40.696210999999999], [-74.1840209999997,40.697074999999989].

然后我使用下面的函数应用字符串连接,基本上是:

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

enter image description here

最后,我使用方法wkt.loads("MultiPolygon ((()))").area来计算多多边形对象的形状区域

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

enter image description here

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

IllegalArgument例外:LinearRing的点不形成封闭的行串

What I have noticed between the two rows, is that the multipolygon of Newark Airport is a continues object of list[list[f64]]] coordinates. Whereas, the Jamaika Bay has multiple sublists [list[list[f64]]] elements (check column coordinates to verify this). Also the screenshot below verify this statement. enter image description here

因此,在应用WKT之前,有没有办法将牙买加湾的多多边形统一到一个地理对象中?

PS: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)

推荐答案

我最终找到了这个问题的解决方案.正如已经描述的那样,我正在拉平坐标列表(又名点).地理空间数据中的点是(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问题上附加的SON文件加载的Polars DataFrame.

Python相关问答推荐

从流程获取定期更新

有没有办法清除气流中的僵尸

过载功能是否包含Support Int而不是Support Int?

遵循轮廓中对象方向的计算线

Chatgpt API不断返回错误:404未能从API获取响应

Locust请求中的Python和参数

Python会扔掉未使用的表情吗?

Pandas 第二小值有条件

如何比较numPy数组中的两个图像以获取它们不同的像素

将jit与numpy linSpace函数一起使用时出错

从收件箱中的列中删除html格式

删除所有列值,但判断是否存在任何二元组

使可滚动框架在tkinter环境中看起来自然

使用setuptools pyproject.toml和自定义目录树构建PyPi包

如何在Raspberry Pi上检测USB并使用Python访问它?

如何在Polars中从列表中的所有 struct 中 Select 字段?

如何指定列数据类型

try 检索blob名称列表时出现错误填充错误""

如何排除prefecture_related中查询集为空的实例?

跳过嵌套JSON中的级别并转换为Pandas Rame