How to get Latitude/Longitude

[Benyuel]

How do you get the latitude/longitude and height from a valid cityjson?

Example

Download NewYork.zip from here

Unzip

import json

blob = json.load('DA13_3D_Buildings_Merged.json')

blob['metadata']
# {'geographicalExtent': [985320.346, 219464.155, -1.5, 1009996.302, 259992.617, 677.862]}

blob['transform']
# {'scale': [0.001, 0.001, 0.001], 'translate': [985320.346, 219464.155, -1.5]}

example_building = 'gml_7KKGXJDCBAXCTELB70OMHK04I376NBA6NGHC'
example_city_obj = blob['CityObjects'][example_building]
# {'type': 'Building', 'attributes': {'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, 'geometry': [{'type': 'MultiSurface', 'boundaries': [[[327517, 327518, 327519, 327520, 327521, 327522, 327523, 327524, 327525, 327526, 327527, 327528, 327529, 327530]], [[327531, 327532, 327533, 327534, 327535, 327536, 327537, 327538, 327539]], [[327540, 327541, 327542, 327543, 327544, 327545, 327546]], [[327544, 327543, 327532, 327531]], [[327530, 327531, 327539, 327517]], [[327529, 327545, 327544, 327530]], [[327528, 327546, 327545, 327529]], [[327527, 327540, 327546, 327528]], [[327526, 327541, 327540, 327527]], [[327525, 327542, 327541, 327526]], [[327524, 327543, 327542, 327525]], [[327523, 327533, 327532, 327524]], [[327522, 327534, 327533, 327523]], [[327521, 327535, 327534, 327522]], [[327520, 327536, 327535, 327521]], [[327519, 327537, 327536, 327520]], [[327518, 327538, 327537, 327519]], [[327517, 327539, 327538, 327518]]], 'semantics': {'values': [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17], 'surfaces': [{'type': 'GroundSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'RoofSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'RoofSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}, {'type': 'WallSurface', 'SOURCE_ID': '13210006270', 'BIN': '1061273', 'DOITT_ID': '168309'}]}, 'lod': 2}]}

From my understanding, each of the geometry.boundaries is an index for the blob.vertices, so to get the vertex for the first element of the first boundary above 327517, we can do

blob['vertices'][327517]
# [13360760, 19999727, 122068]

Based on the cityjson spec, we also now need to use the blob.transform on that.

x_in, y_in, z_in = blob['vertices'][327517]
x_out = (x_in * blob['transform']['scale'][0]) + blob['transform']['translate'][0]
y_out = (y_in * blob['transform']['scale'][1]) + blob['transform']['translate'][1]
z_out = (z_in * blob['transform']['scale'][2]) + blob['transform']['translate'][2]
x_out, y_out, z_out
# 998681.106, 239463.88199999998, 120.568

This is not a lat/lon/height, so am I doing something wrong or are there additional steps in order to get the latitude and longitude?

Thanks in advance for any help!

[kenohori]

This is a projected CRS. My guess would be NAD83. So, you'd need to convert the coordinates using a library like PROJ.4.

[Benyuel]

ok thanks @kenohori!, I'll try that

[hugoledoux]

to do so you can use cjio with the functions 'assign_epsg' (to assign one to the file; I reckon it's that one https://epsg.io/2263) and then you can use the function 'reproject' to epgs:4326.