Splitting A Multipolygon

Can I just take out the parts and take them out to be their own features or would this involve something more complicated? I'm trying to split one of these maps to smaller pieces t

Solution 1:

Extracting the child geometries from GeometryCollection in a GeoJSON file is straightforward but the MultiPolygon is a little tricky.

MultiPolygon geometry as defined in RFC7946 is an array of Polygon coordinate arrays, so must iterate over this array of arrays to get each of the polygons. Note that each polygon may have multiple parts with the first being the exterior ring and the others as interior rings or holes.

{"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]]]]}

The following will extract each child polygon from a MultiPolygon into a separate feature and likewise extract each child geometry from a GeometryCollection into a separate feature. Properties of the parent feature will be inherited and put on the extracted features. A "parent" property is added to flag if the feature came from a MultiPolygon or a GeometryCollection.

import json
import sys

# Split apart geojson geometries into simple geometries

features = []

defadd_feature(props, geom):
    feature = {
        "type": "Feature",
        "geometry": geom
    }
    if props:
        feature["properties"] = props
    features.append(feature)

defprocess_polygon(props, coords, parent):
    print(">> export polygon")
    if parent:
        if props isNone:
            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
        })

defcheck_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 isNone:
                props = dict()
            else:
                props = props.copy()
            props["parent"] = parent
        add_feature(props, geom)

defcheck_feature(f):
    geom = f.get("geometry")
    if geom isNone:
        print("--\nFeature: type=unknown")
        return
    geom_type = geom.get("type")
    props = f.get("properties")
    check_geometry(f, geom_type, geom, props)

defoutput_geojson():
    iflen(features) == 0:
        print("no features to extract")
        returnprint("\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
    withopen("out.geojson", "w") as fp:
        json.dump(geo, fp, indent=2)

############################################################################### main ###############################################################################iflen(sys.argv) >= 2:
    input_file = sys.argv[1]
else:
    print("usage: geojson-split.py <file>")
    sys.exit(0)

print("Input:", input_file)

withopen(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()

If you wanted to create a separate GeoJSON file for each geometry then the add_feature() function could create a new file rather than append feature to the list.

If you run the script on earth-seas-10km.geo.json downloaded from geomaps then output is as follows:

Output:

Input: earth-seas-10km.geo.json
geometry type GeometryCollection
GeometryCollection
MultiPolygon > polygons=21>> export polygon
>> export polygon
...
>> export polygon

Number of extracted features:21Output: out.geojson