Блог

Могу ли я просто вынуть части и сделать их своими собственными функциями или это потребует чего-то более сложного?

Я пытаюсь разделить одну из этих карт на более мелкие части, чтобы проиндексировать их: https://github.com/simonepri/geo-maps

На верхнем уровне они представляют собой геометрическую коллекцию, но на более низких уровнях есть мультиполигоны.

Я думал о том, чтобы просто перебирать части MultiPolygon, но я недостаточно осведомлен о них, чтобы знать, действительно ли это

Я пробовал geojsplit , но, похоже, это не работает для GeometryCollections

Извлечение дочерних геометрий из GeometryCollection в файле GeoJSON просто, но MultiPolygon немного сложнее.

Геометрия MultiPolygon, определенная в RFC7946, представляет собой массив массивов координат полигонов, поэтому необходимо выполнить итерацию по этому массиву массивов, чтобы получить каждый из полигонов. Обратите внимание, что каждый многоугольник может иметь несколько частей, первая из которых является внешним кольцом, а остальные — внутренними кольцами или отверстиями.

 {
       "type": "MultiPolygon",
       "coordinates": [
           [
               [
                   [180.0, 40.0], [180.0, 50.0], [170.0, 50.0],
                   [170.0, 40.0], [180.0, 40.0]
               ]
           ],
           [
               [
                   [-170.0, 40.0], [-170.0, 50.0], [-180.0, 50.0],
                   [-180.0, 40.0], [-170.0, 40.0]
               ]
           ]
       ]
   }
  

Следующее будет извлекать каждый дочерний полигон из MultiPolygon в отдельный объект и аналогичным образом извлекать каждую дочернюю геометрию из GeometryCollection в отдельный объект. Свойства родительского объекта будут унаследованы и помещены в извлеченные объекты. Свойство «parent» добавляется к флагу, если функция была получена из MultiPolygon или GeometryCollection.

 import json
import sys

# Split apart geojson geometries into simple geometries

features = []

def add_feature(props, geom):
    feature = {
        "type": "Feature",
        "geometry": geom
    }
    if props:
        feature["properties"] = props
    features.append(feature)

def process_polygon(props, coords, parent):
    print(">> export polygon")
    if parent:
        if props is None:
            props = dict()
        else:
            props = props.copy()
        # mark feature where it came from if extracted from MultiPolygon or GeometryCollection
        props["parent"] = parent
    add_feature(props, {
            "type": "Polygon",
            "coordinates": coords
        })

def check_geometry(f, geom_type, geom, props, parent=None):
    if geom_type == "Polygon":
        coords = geom["coordinates"]
        print(f"Polygon > parts={len(coords)}")
        process_polygon(props, coords, parent)
    elif geom_type == "MultiPolygon":
        coords = geom["coordinates"]
        print(f"MultiPolygon > polygons={len(coords)}")
        for poly in coords:
            process_polygon(props, poly, "MultiPolygon")
    elif geom_type == 'GeometryCollection':
        print("GeometryCollection")
        """
        "type": "GeometryCollection",
        "geometries": [{
             "type": "Point",
             "coordinates": [101.0, 1.0]
         }, {
            "type": "Polygon",
            "coordinates": [                    
                        [[ 100, 0 ], [ 100, 2 ], [ 103, 2 ], [ 103, 0 ], [ 100, 0 ]]
            ]
            }]
        """
        for g in geom["geometries"]:
            check_geometry(f, g.get("type"), g, props, 'GeometryCollection')
    else:
        # other geometry type as-is: e.g. point, line, etc.
        print(">> add other type:", geom_type)
        if parent:
            if props is None:
                props = dict()
            else:
                props = props.copy()
            props["parent"] = parent
        add_feature(props, geom)

def check_feature(f):
    geom = f.get("geometry")
    if geom is None:
        print("--nFeature: type=unknown")
        return
    geom_type = geom.get("type")
    props = f.get("properties")
    check_geometry(f, geom_type, geom, props)

def output_geojson():
    if len(features) == 0:
        print("no features to extract")
        return

    print("nNumber of extracted features:", len(features))
    print("Output: out.geojson")
    geo = {
        "type": "FeatureCollection",
    }
    for key in ["name", "crs", "properties"]:
        if key in data:
            geo[key] = data[key]
    geo["features"] = features
    with open("out.geojson", "w") as fp:
        json.dump(geo, fp, indent=2)

##############################################################################
# main #
##############################################################################

if len(sys.argv) >= 2:
    input_file = sys.argv[1]
else:
    print("usage: geojson-split.py <file>")
    sys.exit(0)

print("Input:", input_file)

with open(input_file, encoding="utf-8") as file:
    data = json.load(file)
data_type = data.get("type")
if data_type == "FeatureCollection":
    for f in data['features']:
        obj_type = f.get("type")
        if obj_type == "Feature":
            check_feature(f)
        else:
            print("WARN: non-feature:", obj_type)
elif data_type == "Feature":
    check_feature(data)
elif data_type in ['GeometryCollection', 'Point', 'LineString', 'Polygon',
                   'MultiPoint', 'MultiLineString', 'MultiPolygon']:
    print("geometry type", data_type)
    check_geometry(data, data_type, data, None)
else:
    print("WARN: unknown/other type:", data_type)

output_geojson()
  

Если вы хотите создать отдельный файл GeoJSON для каждой геометрии, тогда add_feature() функция может создать новый файл, а не добавлять объект в список.

Если вы запустите скрипт на earth-seas-10km.geo.json, загруженный из geomaps, то результат будет следующим:

Вывод:

 Input: earth-seas-10km.geo.json
geometry type GeometryCollection
GeometryCollection
MultiPolygon > polygons=21
>> export polygon
>> export polygon
...
>> export polygon

Number of extracted features: 21
Output: out.geojson