拆分多重多边形

问题描述

我能不能只是把零件拿出来作为它们自己的特征，或者这是否会涉及到更复杂的东西？

我正在尝试将其中一个地图拆分成更小的部分，以便为它们编制索引：https://github.com/simonepri/geo-maps

在顶层，它们是GeometryCollection，但在较低级别，它们是多重多边形。

我想只循环遍历多重多边形的各个部分，但我对它们的了解不够，无法知道这是否有效

我已尝试geojsplit，但它似乎不适用于GeometryCollection

---

解决方案

从GeoJSON文件中的GeometryCollection中提取子几何图形很简单，但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]
               ]
           ]
       ]
   }

下面将把多重多边形中的每个子多边形提取到单独的要素中，并同样将GeometryCollection中的每个子几何图形提取到单独的要素中。父特征的属性将被继承并放置在提取的特征上。如果要素来自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("--
Feature: 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("
Number 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()函数可以创建新文件，而不是将要素追加到列表。

如果您在从geomaps下载的Earth-Sea-10km.Geo.json上运行该脚本，则输出如下：

输出：

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