分割多边形

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

如果你运行 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

分割多边形

Splitting a MultiPolygon

python

polygon

geojson

python-3.x

elasticsearch