分裂一个MultiPolygon

Question

我能不能把这些零件拿出来，让它们成为自己的特征，还是会涉及到更复杂的事情？

我试图将这些地图中的一张分割成更小的部分来索引：https://github.com/simonepri/geo-maps

在顶层，它们是GeometryCollection，但在较低的级别，有MultiPolygons。

我正在考虑在MultiPolygon的各个部分循环，但是我对它们还不够了解，不知道这是否有效。

我试过geojsplit，但它似乎不适用于GeometryCollections

Answer 1

从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中提取成一个单独的特性。父特征的属性将被继承并放在提取的特征上。如果特性来自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()函数可以创建一个新文件，而不是在列表中添加特性。

如果在从风水下载的地球海洋-10 km.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