How to calculate building volume using 3D BAG and PostgreSQL?

Datasets BAG
#1

Hi, I’m trying to calculate the volumes of buildings. I have done this in the past by generating point clouds from 3D BAG and AHN data. Now, I want to use the 3D geometries provided by the 3D BAG PostgreSQL database dump. Building geometries are provided as MultiPolygonZ, and currently the highest LoD is 2.2 (source: General Concepts - 3D BAG).

I’m having trouble with calling ST_Volume() on the geometries (source: ST_Volume). I have the postgis_sfcgal extension installed. Here’s what I tried:

SELECT
    -- ST_Volume(geometrie),
    -- ST_Volume(ST_MakeValid(geometrie)),
    ST_Volume(ST_MakeSolid(geometrie)),
    -- ST_Volume(ST_Tesselate(geometrie)),
    -- ST_Volume(ST_Buffer(ST_Buffer(geometrie, -0.001), 0.001)),
    gid
FROM bag3d.lod22_3d LIMIT 1;

Doing ST_Volume(geometrie), ST_Volume(ST_MakeSolid(geometrie)), or ST_Volume(ST_Tesselate(geometrie)) results in:

[XX000] ERROR: MultiPolygon is invalid : intersection between Polygon 0 and 10

Doing ST_Volume(ST_MakeValid(geometrie)) results in:

[XX000] ERROR: GeometryCollection is invalid : MultiPolygon 0 is invalid: Polygon 1 is invalid: points don't lie in the same plane

Doing ST_Volume(ST_Buffer(ST_Buffer(geometrie, -0.001), 0.001)) does sometimes result in the error below, and otherwise always yields a volume of zero.

[XX000] ERROR: Intersection does not support polygon with connected rings

When I split the MultiPolygonZs into Polygons, I get similar results:

SELECT
    n,
    ST_Volume(ST_GeometryN(geometrie, n)),
    ST_Volume(ST_MakeValid(ST_GeometryN(geometrie, n))),
    ST_Volume(ST_MakeSolid(ST_GeometryN(geometrie, n))),
    ST_Volume(ST_Tesselate(ST_GeometryN(geometrie, n))),
    ST_Volume(ST_Buffer(ST_Buffer(geometrie, -0.001), 0.001)),
    gid
FROM (SELECT * FROM bag3d.lod22_3d LIMIT 1) a
CROSS JOIN generate_series(1, 100) n
WHERE n <= ST_NumGeometries(geometrie);

None of these combinations result in an error, but instead, they always yield a volume of zero.

In all cases, other postgis_sfcgal functions such as ST_3DArea() do return some useful output. Does this mean that I have to convert the MultiPolygonZs into POLYHEDRALSURFACE or POLYHEDRALSURFACEZ first, before I can call ST_MakeSolid() on them? Or am I missing something else?

Any help with this is greatly appreciated!

#2

It seems like there are some geometric errors in the source data. Have you tried making a small piece of testdata where you know the geometry is proper? In that way you can test the method of calculating the volume of certain 3D polygon.
Seems like you won’t find a solution when you are testing on “unusable” data.

#3

Thank you for your reply!

Unfortunately, none of the MultiPolygonZ geometries in the entire 3D BAG database return true for either ST_IsValid() or ST_IsSolid(), even if I do ST_MakeValid() or ST_MakeSolid() first.

However, if I work through some examples from the Postgis docs, it seems like the only way to calculate a volume is from a POLYHEDRALSURFACE or POLYHEDRALSURFACEZ.

I am starting to think that the proper way to phrase my question would be:

How can I convert a MultiPolygonZ to a POLYHEDRALSURFACEZ?

#4

I’m afraid I can’t help you there (don’t know that much about PostGreSQL). In Oracle you can convert geometry into some other definition and possibly with some modifications/corrections convert it back into the original geometry. Not in possession of FME by any chance?

Also, maybe the problem can be resolved using an application that doesn’t see the geometry as faulty? Perhaps using QGIS to calculate area, store it as an attribute or directly run a calculation on it?

#5

Thanks @ademeijer for the suggestions!

I actually tried loading the data into Python, but calculating volume from 3D geometries (e.g., using Shapely and SciPy) apparently is not straightforward at all!

However, I think I found a solution workaround. To answer my own question (‘How can I convert a MultiPolygonZ to a POLYHEDRALSURFACEZ?’), in a very hacky way:

SELECT
    ST_Volume(ST_MakeSolid(ST_GeomFromText(
        REPLACE(
            ST_AsText(bag3d.lod22_3d.geometrie),
            'MULTIPOLYGON Z', 'POLYHEDRALSURFACEZ'
        )))) AS volume,
    bag3d.pand.identificatie AS pand_id
FROM bag3d.lod22_3d
JOIN bag3d.pand
ON bag3d.lod22_3d.fid = bag3d.pand.fid
LIMIT 1;

Using this to calculate some known volumes of buildings actually seems to give accurate results.

Also thanks to @BalazsDukai and the 3D BAG team for providing this dataset!

1 like
#6

I want to add, for anyone else stumbling upon this, you’d probably want to add

WHERE cardinality(bag3d.lod22_3d.val3dity_codes) = 0

to your query.

#7

For anyone reading this, make sure you first read the relevant part of the 3DBAG documentation. Geometric validity in 3D is not as simple as in 2D, that is why val3dity exists. PostGIS’s MakeValid doesn’t work for the 3D models.
Anything below code 306 (including) will probably give incorrect results for volume calculations.

@petersaalbrink Shapely is 2D only. I guess SciPy could work, but I never tried that myself.