Splitting A Multipolygon
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
Post a Comment for "Splitting A Multipolygon"