LoD1 terrain

[AdrianKriger]

I am sorry to bother you with this trivial matter but I need your help.

I would like to create a basic LoD1 cityjson terrain model but am unable too.
Vertices and geometry boundaries (indices / simplices) are via a Constrained Delaunay with python bindings of Shewchuck's; Triangle [option 'p' with holes - concavities removed].

Wavefront.obj is successfully created but I fail with cityjson. The code is:

def output_citysjon(extent, minz, maxz, T, pts):
    """
    basic function to produce LoD1 cityjson terrain
    """
    #-- create the JSON data structure for the City Model
    cm = {}
    cm["type"] = "CityJSON"
    cm["version"] = "0.9"
    cm["CityObjects"] = {}
    cm["vertices"] = []
    #-- Metadata is added manually
    cm["metadata"] = {
    "datasetTitle": "LoD1 terrain model of CPUT (Bellville) campus",
    "datasetReferenceDate": "2021-07-31",
    "geographicLocation": "Cape Town, South Africa",
    "referenceSystem": "urn:ogc:def:crs:EPSG::32733",
    "geographicalExtent": [
        extent[0],
        extent[1],
        minz ,
        extent[1],
        extent[1],
        maxz
      ],
    "datasetPointOfContact": {
        "contactName": "arkriger",
        "linkedin": "www.linkedin.com/in/adrian-kriger",
        "contactType": "private",
        "github": "https://github.com/AdrianKriger/osm_LoD1_3Dbuildings"
        },
    "metadataStandard": "ISO 19115 - Geographic Information - Metadata",
    "metadataStandardVersion": "ISO 19115:2014(E)"
    }
    
    add_terrain_v(pts, cm)
    grd = {}
    grd['type'] = 'TINRelief'
    grd['geometry'] = [] #-- a cityobject can have >1 
     #-- the geometry
    g = {} 
    g['type'] = 'CompositeSurface'
    g['lod'] = 1
    allsurfaces = [] #-- list of surfaces
    add_terrain_b(T, allsurfaces)
    g['boundaries'] = []
    g['boundaries'].append(allsurfaces)
      #-- add the geom 
    grd['geometry'].append(g)
      #-- insert the terrain as one new city object
    cm['CityObjects']['terrain01'] = grd

    return cm

def add_terrain_v(pts, cm):
    for p in pts:
        cm['vertices'].append([[p[0], p[1], p[2]]])
    
def add_terrain_b(T, allsurfaces):
    for i in T:
        allsurfaces.append([i[0] + 1, i[1] + 1, i[2] + 1])

where T and pts are the simplices and points respectively.
The complete workflow is here.
Your guidance is appreciated.

[liberostelios]

Hi Adrian. Glad to see you are interested in CityJSON.

I think your issue is in add_terrain_v, where there seems to be one extra square bracket wrapped around the vertex coordinates. The right code should be:

def add_terrain_v(pts, cm):
    for p in pts:
        cm['vertices'].append([p[0], p[1], p[2]]) # <- There was an extra square bracket here

Haven't checked thoroughly the rest of the code, but this for sure is an issue.

PS: It's very nice to see that you have some nice metadata in your model! 👏

[liberostelios]

That formula is awesome! 😅

Looking at the code again, I think you are also missing a bracket from add_terrain_b:

def add_terrain_b(T, allsurfaces):
    for i in T:
        allsurfaces.append([[i[0] + 1, i[1] + 1, i[2] + 1]]) # <- Every polygon is a list of rings (one outer + holes)

If that doesn't work, i'd suggest that you:

• validate your file with cjio to get some extra hints, or
• upload the CityJSON file (or a small sample) here so we can take a look at the result.

[AdrianKriger]

We appear to be getting closer to the solution. With:

from cjio import cityjson
cm = cityjson.load('citjsn_cput3d.json')

the result is:

TypeError: Boundary definition does not correspond to MultiSurface or CompositeSurface

I do not understand # <- Every polygon is a list of rings (one outer + holes).
The geometry I want to define is a TINRelief. A polygon with 3 vertices, and no interior ring. Would the boundaries not be the faces (indices/simplices) of the triangulation? Must the constraints be defined as well? or am I defining the CityObject incorrectly?

iether way I can:

cjio citjsn_cput3d.json upgrade_version save citjsnV1_cput3d.json

and have placed the result here.
then:

cjio citjsnV1_cput3d.json validate

the result:

-- Validating the syntax of the file
        (using the schemas 1.0.1)
=====
File is invalid
File has no warnings
--- ERRORS (total = 2) ---
1 ==> {'contactName': 'arkriger', 'linkedin': 'www.linkedin.com/in/adrian-kriger', 'contactType': 'private', 'github': 'https://github.com/AdrianKriger/osm_LoD1_3Dbuildings'} is not valid under any of the given schemas

2 ==> CityObject is not schema-valid: terrain01

Lastly: how do I call upgrade_version with the cjio API?

You are being extremely helpful. Thank you.

[liberostelios]

the result is:

TypeError: Boundary definition does not correspond to MultiSurface or CompositeSurface

Okay, that's indicative that there is indeed a problem with the depth of the lists in your boundaries. I just checked your data and realised that you are actually adding one extra list in the structure. This is because you do this:

g['boundaries'] = []
g['boundaries'].append(allsurfaces)

Instead, you should just do this:

g['boundaries'] = allsurfaces

After removing the outer brackets from the result I managed to visualise it in azul (although the output is messy) but not in ninja. I get an error that one of the vertices is out of bounds. I suspect that this has to do with you adding a + 1 for every index, you might want to do this:

def add_terrain_b(T, allsurfaces):
    for i in T:
        allsurfaces.append([[i[0], i[1], i[2]]])

CityJSON uses a 0-based array notation.

I do not understand # <- Every polygon is a list of rings (one outer + holes).
The geometry I want to define is a TINRelief. A polygon with 3 vertices, and no interior ring. Would the boundaries not be the faces (indices/simplices) of the triangulation? Must the constraints be defined as well? or am I defining the CityObject incorrectly?

So, for a MultiSurface or CompositeSurface you define the boundary as a list of faces. Every face, is composed by a list of rings; and every ring is a list of points (i.e. indices). The means, that for simple triangles you need to have something like this:

boundaries = [ # This is the list of faces
    [ # -- This is the list of rings
        [0, 1, 2] # -- Your triangle indices
    ]
]

So, you might not have holes in your dataset, but you still need to represent every triangle as a list of rings with only one ring. In your code you were representing a face as a list of indices, ignoring the list of rings. I hope that this makes sense.

iether way I can:

cjio citjsn_cput3d.json upgrade_version save citjsnV1_cput3d.json

and have placed the result here.

Lastly: how do I call upgrade_version with the cjio API?

Afair, you just call cm.upgrade_version("1.0"). You may want to check the source code here. Unfortunately, the API is not very thoroughly documented, so we should look into that in the future.

[AdrianKriger]

You have been extremely kind and helpful. result. Thank you.

[hugoledoux]

to upgrade, you can just call cm.upgrade_version() and it will upgrade to the next version.

but it would be better to just create v1.0 file directly! Just change that line and I think all is okay:
cm["version"] = "0.9"
=>
cm["version"] = "1.0"