#!/usr/bin/env python ############################################################################### # # Purpose: CityJSON OGR driver # Author: Even Rouault # ############################################################################### # Copyright (c) 2019, Even Rouault # SPDX-License-Identifier: MIT ############################################################################### # Metadata parsed by GDAL C++ code at driver pre-loading, starting with '# gdal: ' # gdal: DRIVER_NAME = "CityJSON" # gdal: DRIVER_SUPPORTED_API_VERSION = [1] # gdal: DRIVER_DCAP_VECTOR = "YES" # gdal: DRIVER_DCAP_VIRTUALIO = "YES" # gdal: DRIVER_DMD_LONGNAME = "CityJSON" # gdal: DRIVER_DMD_EXTENSIONS = "json" import json import os from gdal_python_driver import BaseDataset, BaseDriver, BaseLayer class Layer(BaseLayer): def __init__(self, filename, content): self.content = content self.name = os.path.splitext(os.path.basename(filename))[0] self.fields = [ {"name": "id"}, {"name": "type"}, ] self.count = 0 set_attrs = set() for key_value in content["CityObjects"].items(): value = key_value[1] if "attributes" in value: for kv in value["attributes"].items(): if not kv[0] in set_attrs: set_attrs.add(kv[0]) v = kv[1] t = "String" if isinstance(v, int): t = "Integer" elif isinstance(v, float): t = "Real" elif "Date" in kv[0]: t = "Date" self.fields.append({"name": kv[0], "type": t}) for geom in value["geometry"]: if geom["type"] in ("Solid", "MultiSurface"): self.count += 1 md = content["metadata"] srs = md["referenceSystem"] if "referenceSystem" in md else None if not srs and "crs" in md: if isinstance(md["crs"], dict) and "epsg" in md["crs"]: srs = "EPSG:" + str(md["crs"]["epsg"]) self.geometry_fields = [{"type": "MultiPolygonZ", "srs": srs}] if "transform" in content: self.translate = content["transform"]["translate"] self.scale = content["transform"]["scale"] else: self.translate = [0, 0, 0] self.scale = [1, 1, 1] self.vertices = content["vertices"] def extent(self, force): md = self.content["metadata"] if "geographicalExtent" in md: bbox = md["geographicalExtent"] else: bbox = md["bbox"] return [bbox[0], bbox[1], bbox[3], bbox[4]] def feature_count(self, force): return self.count def __iter__(self): fid = 1 vertices = self.vertices scale = self.scale translate = self.translate for key_value in self.content["CityObjects"].items(): value = key_value[1] for geom in value["geometry"]: geom_in = None boundaries = geom["boundaries"] if geom["type"] == "Solid": assert len(boundaries) == 1, (len(boundaries), key_value) geom_in = boundaries[0] elif geom["type"] == "MultiSurface": geom_in = boundaries if geom_in: out_mpoly = [] for poly in geom_in: out_poly = [] for ring in poly: out_ring = [] for vertex_idx in ring: v = vertices[vertex_idx] c = [v * s + t for v, s, t in zip(v, scale, translate)] out_ring.append(c) out_ring.append(out_ring[0]) out_poly.append(out_ring) out_mpoly.append(out_poly) properties = {"id": key_value[0], "type": value["type"]} if "attributes" in value: properties.update(value["attributes"]) wkt = "MULTIPOLYGON Z(" first_poly = True for poly in out_mpoly: if not first_poly: wkt += "," else: first_poly = False wkt += "(" first_ring = True for ring in poly: if not first_ring: wkt += "," else: first_ring = False wkt += "(" first_vertex = True for v in ring: if not first_vertex: wkt += "," else: first_vertex = False wkt += "%.18g %.18g %.18g" % (v[0], v[1], v[2]) wkt += ")" wkt += ")" wkt += ")" f = { "type": "OGRFeature", "id": fid, "fields": properties, "geometry_fields": {"": wkt}, } yield f fid += 1 class Dataset(BaseDataset): def __init__(self, filename, content): self.layers = [Layer(filename, content)] class Driver(BaseDriver): def identify(self, filename, first_bytes, open_flags, open_options={}): return ( b'"type":"CityJSON"' in first_bytes or b'"type": "CityJSON"' in first_bytes ) def open(self, filename, first_bytes, open_flags, open_options={}): if not self.identify(filename, first_bytes, open_flags): return None if filename.startswith("/vsi"): from osgeo import gdal f = gdal.VSIFOpenL(filename, "rb") if not f: return None gdal.VSIFSeekL(f, 0, 2) length = gdal.VSIFTellL(f) gdal.VSIFSeekL(f, 0, 0) try: content = json.loads(gdal.VSIFReadL(1, length, f)) finally: gdal.VSIFCloseL(f) else: with open(filename, "rt") as f: content = json.loads(f.read()) return Dataset(filename, content)