footprints.py

import json
from cjio.cityjson import CityJSON
from cjio import cityjson, subset
from shapely import geometry
from shapely.ops import triangulate
from shapely.geometry import mapping

def create_footprints(inpath, outpath):
    cmf = open(inpath, 'r')
    cm = cityjson.reader(file=cmf, ignore_duplicate_keys=True)
    
    geojson = {
       "type": "FeatureCollection",
       "features": []
       }
    
    cm.decompress()

    cj_vertices = cm.j["vertices"]
    
    for co in cm.j["CityObjects"]:
        for geom in cm.j["CityObjects"][co]["geometry"]:
            geom_2d = []
            geom_to_shapely(cj_vertices, geom["boundaries"], geom_2d)
            geom_2d.append(geom_2d[0])
            
            p = geometry.Polygon(geom_2d).convex_hull
            buf = p.buffer(0)
            
            mapped = mapping(buf)
            feature = {"type": "Feature", "geometry": mapped, "properties": cm.j["CityObjects"][co]["attributes"]}
            geojson["features"].append(feature)
            
    fout = open(outpath, "w")
    # separators given as argument to remove whitespace
    json.dump(geojson, fout)
        
        
def geom_to_shapely(cj_vertices, boundaries, geom_2d):
    for b in boundaries:
        if any(isinstance(i, list) for i in b):
            geom_to_shapely(cj_vertices, b, geom_2d)
        elif b is not None:
            for v_i in b:
                v = cj_vertices[v_i]
                geom_2d.append([v[0], v[1]])
                
                
inpath = "C:/Users/jordi/Documents/3dgeoinfo/data/subset_utrecht.json"
outpath = "C:/Users/jordi/Documents/3dgeoinfo/data/subset_utrecht_geojson.json"
create_footprints(inpath, outpath)