Создайте объекты MultiPolygon из типа структуры и столбцов list[list[list[f64]]] с помощью Polars

Я загрузил набор данных «Зоны такси Нью-Йорка» (загружен из SODA Api и сохранен как файл json, а не GeoJson или Shapefile). Набор данных довольно мал, поэтому я использую всю включенную информацию. Для удобства поста представляю первые 2 строки набора данных:

• исходный (со значением типа структуры the_geom).
• набор данных после распаковки типа структуры с помощью команды unnest() в полярах. --обновлено командой write_ndjson()

Исходный набор данных:

Ниже набора данных после применения unnest() и выбора некоторых столбцов и первых двух строк.

Вы можете импортировать данные, используя поляры, с помощью следующей команды

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

Меня интересуют мультиполигоны. На самом деле я пытаюсь пересчитать shape_area, используя мультиполигон и wkt (представление хорошо известного текста) - метод, используемый модулем 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.69518499999987], [-7], 4.18438099999997, 40.69587799999989], [-74.18428199999994, 40.6962109999999], [-74.18402099999997, 40.69707499999989]...

Затем я использую функцию ниже, которая применяет конкатенацию строк и в основном:

• Преобразует объект list[list[f64]] в String.
• Добавляет ключевое слово 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]]] (проверьте координаты столбца, чтобы убедиться в этом). Также скриншот ниже подтверждает это утверждение.

Таким образом, есть ли способ объединить мультиполигоны залива Ямайка в один объект GEOmetric перед применением WKT?

P.S. Многие решения на GitHub используют файл формы, но я хотел бы регулярно автоматически повторно загружать зоны Нью-Йорка с помощью SODA API.

Чтобы загрузить необработанный файл .json из SODA API (опустите 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)