分割多面体

问题描述

我可以只取出零件并将它们取出作为自己的特征吗，或者这会涉及更复杂的事情吗？

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

顶层是GeometryCollection，而底层是Multipolygons。

我当时只是想遍历Multipolygon的各个部分，但是我对它们的了解不足以知道这是否有效

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

解决方法

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

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

{
       "type": "MultiPolygon","coordinates": [
           [
               [
                   [180.0,40.0],[180.0,50.0],[170.0,40.0]
               ]
           ],[
               [
                   [-170.0,[-170.0,[-180.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,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]
         },"coordinates": [                    
                        [[ 100,0 ],[ 100,2 ],[ 103,0 ]]
            ]
            }]
        """
        for g in geom["geometries"]:
            check_geometry(f,g.get("type"),g,'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,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(f"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