3DBAG filteren om 3D geometrie van dakvlakken/muurvlakken te verkrijgen

We hebben de geopackage, cityjson, WFS en PostGreSQL dump van 3DBAG bekeken, maar konden nog niet eenvoudig de 3D geometrie van dakvlakken en muurvlakken daaruit krijgen. We lazen dat 3D geometrie van dakvlak/muurvlak semantisch al aangeduid zou moeten zijn bij LOD 2.2 in cityjson als ‘RoofSurface’/‘WallSurface’ en in viewers (zoals ninja: CityJSON Ninja met file: https://data.3dbag.nl/cityjson/v210908_fd2cee53/3dbag_v210908_fd2cee53_5907.json) is dakvlak heel duidelijk rood aangeduid. Dus we hopen dat er eenvoudige en standaard bepaling is om dakvlak en muurvlak uit 3DBAG te krijgen?

Weet iemand wat makkelijke manier (FME/Python/ArcGIS/PostGIS/?) is om 3D geometrie van dakvlak/muurvlak uit 3DBAG te krijgen?

Met nieuwe 3DBAG van oktober 2023 is het zeer eenvoudig om dak/muur/vloerdelen te filteren door nieuw labels kenmerk (Data Attributes - 3DBAG):
(0: GroundSurface , 1: RoofSurface , 2: OuterWallSurface , 3: InnerWallSurface )

Dit is voorbeeld Python code o.b.v. 3D geometrie in gpkg:

from shapely import wkt
pand_geom_3d_labels = "(14:1,2,2,2,2,2,2,2,2,2,2,2,0,0)"
pand_geom_3d = "MULTIPOLYGON Z (((70033.875 379772.46875 1.195, 70022.757812 379758.0625 1.195, 69994.867188 379779.53125 1.195, 70005.890625 379793.90625 1.195, 70022.101562 379815 1.195, 70030.234375 379808.75 1.195, 70027.65625 379805.4375 1.195, 70035.5625 379799.21875 1.195, 70033.171875 379796.0625 1.195, 70041.007812 379789.90625 1.195, 70041.890625 379782.90625 1.195, 70033.875 379772.46875 1.195)), ((70041.007812 379789.90625 6.226, 70041.007812 379789.90625 1.195, 70033.171875 379796.0625 1.195, 70033.171875 379796.0625 6.203, 70041.007812 379789.90625 6.226)), ((70041.890625 379782.90625 7.655, 70041.890625 379782.90625 1.195, 70041.007812 379789.90625 1.195, 70041.007812 379789.90625 6.226, 70041.890625 379782.90625 7.655)), ((70005.890625 379793.90625 11.447, 70005.890625 379793.90625 1.195, 69994.867188 379779.53125 1.195, 69994.867188 379779.53125 6.206, 70005.890625 379793.90625 11.447)), ((70022.101562 379815 3.864, 70022.101562 379815 1.195, 70005.890625 379793.90625 1.195, 70005.890625 379793.90625 11.447, 70022.101562 379815 3.864)), ((69994.867188 379779.53125 6.206, 69994.867188 379779.53125 1.195, 70022.757812 379758.0625 1.195, 70022.757812 379758.0625 6.124, 69994.867188 379779.53125 6.206)), ((70027.65625 379805.4375 5.052, 70027.65625 379805.4375 1.195, 70030.234375 379808.75 1.195, 70030.234375 379808.75 3.854, 70027.65625 379805.4375 5.052)), ((70035.5625 379799.21875 5.077, 70035.5625 379799.21875 1.195, 70027.65625 379805.4375 1.195, 70027.65625 379805.4375 5.052, 70035.5625 379799.21875 5.077)), ((70022.757812 379758.0625 6.124, 70022.757812 379758.0625 1.195, 70033.875 379772.46875 1.195, 70033.875 379772.46875 11.398, 70022.757812 379758.0625 6.124)), ((70033.171875 379796.0625 6.203, 70033.171875 379796.0625 1.195, 70035.5625 379799.21875 1.195, 70035.5625 379799.21875 5.077, 70033.171875 379796.0625 6.203)), ((70030.234375 379808.75 3.854, 70030.234375 379808.75 1.195, 70022.101562 379815 1.195, 70022.101562 379815 3.864, 70030.234375 379808.75 3.854)), ((70033.875 379772.46875 11.398, 70033.875 379772.46875 1.195, 70041.890625 379782.90625 1.195, 70041.890625 379782.90625 7.655, 70033.875 379772.46875 11.398)), ((70022.757812 379758.0625 6.124, 70033.875 379772.46875 11.398, 70005.890625 379793.90625 11.447, 69994.867188 379779.53125 6.206, 70022.757812 379758.0625 6.124)), ((70033.875 379772.46875 11.398, 70041.890625 379782.90625 7.655, 70041.007812 379789.90625 6.226, 70033.171875 379796.0625 6.203, 70035.5625 379799.21875 5.077, 70027.65625 379805.4375 5.052, 70030.234375 379808.75 3.854, 70022.101562 379815 3.864, 70005.890625 379793.90625 11.447, 70033.875 379772.46875 11.398)))"

pand_polygon = wkt.loads(pand_geom_3d)
pand_polygon_geometries = pand_polygon.geoms
pand_geom_3d_labels_list = pand_geom_3d_labels.split(":")[-1][:-1].split(",")
for pand_polygon_index in range(0, len(pand_polygon_geometries)):
    if pand_geom_3d_labels_list[pand_polygon_index] == "0": # check if pand element is ground surface
        print("Ground surface:", pand_polygon_geometries[pand_polygon_index])
    elif pand_geom_3d_labels_list[pand_polygon_index] == "1": # check if pand element is roof surface
        print("Roof surface:", pand_polygon_geometries[pand_polygon_index])
    elif pand_geom_3d_labels_list[pand_polygon_index] == "2": # check if pand element is outer wall surface
        print("Outer wall surface:", pand_polygon_geometries[pand_polygon_index])