The Ohsome website allows generating statitics but expects the input as a list of coordinates rather than geojson. The first step is to download a goejson shape of the project from the Task Server – look for the button highlighted below (Download AOI):

The geojson file may be a little too big/complex to use on the ohsome website so next I used a python script to minify it a little:

import json
from shapely.geometry import shape, mapping

INPUT = "input.geojson"
OUTPUT = "simplified.geojson"
TOLERANCE = 0.001  # increase for more simplification

with open(INPUT, "r", encoding="utf-8") as f:
    data = json.load(f)

def simplify_geom(geom_obj, tol):
    g = shape(geom_obj)
    g_simplified = g.simplify(tol, preserve_topology=True)
    return mapping(g_simplified)

if data["type"] == "FeatureCollection":
    for feat in data["features"]:
        feat["geometry"] = simplify_geom(feat["geometry"], TOLERANCE)
elif data["type"] in ("Feature",):
    data["geometry"] = simplify_geom(data["geometry"], TOLERANCE)
else:
    # bare geometry
    data = simplify_geom(data, TOLERANCE)

with open(OUTPUT, "w", encoding="utf-8") as f:
    json.dump(data, f)

Then I used another python script to convert the simplified.geojson to the bbox coords needed by Ohsome:

import json
import sys

def iter_coords_multi_polygon(coords):
    # coords structure: [ [ [ [lon, lat], ... ] ] , ... ]
    for polygon in coords:
        for ring in polygon:
            for lon, lat in ring:
                yield lon, lat

def main(path):
    with open(path, "r", encoding="utf-8") as f:
        data = json.load(f)

    # If the file is just a geometry object like in your example:
    if data.get("type") == "MultiPolygon":
        coords = data["coordinates"]
        coord_iter = iter_coords_multi_polygon(coords)
    # If it's a Feature with geometry:
    elif data.get("type") == "Feature" and data.get("geometry", {}).get("type") == "MultiPolygon":
        coords = data["geometry"]["coordinates"]
        coord_iter = iter_coords_multi_polygon(coords)
    else:
        raise ValueError("Expected a MultiPolygon geometry")

    parts = []
    for lon, lat in coord_iter:
        parts.append(f"{lon},{lat}")

    # Single line output
    print(",".join(parts))

if __name__ == "__main__":
    # Usage: python script.py input.geojson
    if len(sys.argv) != 2:
        print(f"Usage: {sys.argv[0]} input.geojson", file=sys.stderr)
        sys.exit(1)
    main(sys.argv[1])

The resulting output can be pasted into the Ohsome website, first choose Polygon:

Then look for the ‘Bounding Polygon’ form field and paste the coords as shown below:


0 Comments

Leave a Reply

Avatar placeholder