分割多面体

Question

我可以把零件拿出来作为自己的功能吗？还是会涉及更复杂的事情？

我正在尝试将其中一张地图拆分为更小的部分以对其进行索引：https://github.com/simonepri/geo-maps

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

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

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

Answer 1

从 GeoJSON 文件中的 GeometryCollection 中提取子几何很简单，但 MultiPolygon 有点棘手。

RFC7946 中定义的 MultiPolygon 几何是多边形坐标数组的数组，因此必须遍历该数组数组以获取每个多边形。请注意，每个多边形可能有多个部分，第一个部分是外环，其他部分是内环或孔。

{
       "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() 函数可以创建一个新文件，而不是将特征附加到列表中。

如果在geomaps下载的earth-seas-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