Support of 3D layered meshes

[PeterPetrik]

QGIS Enhancement: Support of 3D layered meshes

Date 2019/10/16

Author Peter Petrik (@PeterPetrik)

Contact peter.petrik@lutraconsulting.co.uk

maintainer @PeterPetrik

Version QGIS 3.12

Summary

The initial implementation of mesh layer in QGIS supports only basic type of 2D irregular meshes with data defined on faces and vertices. The standard for mesh layers UGRID contains other mesh types that are not yet supported. This proposal describes the required steps to allow usage of 3D Layered Meshes for QGIS users.

Description

3D Layered Mesh consists of multiple stacked 2D unstructured meshes each extruded in the vertical direction (levels) by means of a vertical coordinate. The vertices and faces have the same topology in each vertical level. The mesh definition (vertical level extrusion) could in general change in time. The data is usually defined in volume centres or by some parametric function. Here is an example of regular mesh with 4 levels spaced by 1 meter.

Here is an another example from TUFLOW FV that demonstrates one possible type of vertical layering. When water level increases or decreases the layering stretches or compresses, modifying the absolute elevation in time

The number of vertical levels can't change in time, but it can vary in spatial dimension

Commonly in meteorological/hydrology/geology numerical results, you may get tens of vertical levels and to be able to visualise in on 2D canvas, various averaging methods are used to derive 2D datasets. Also there should be possibility to display just one particular vertical level (see the green quantity on following figure (© Met.3D contributors 2017. CC Attribution 4.0 International License. met3d.wavestoweather.de))

Example Data

There are data defined for each time, vertical level and volume. Volumes are contructed from base 2D mesh and vertical levels.

For example lets have simple 2D mesh with 2 quad elements F1 and F2 and 8 vertices V1 - V8.

Imagine we measure air temperature above the ground. We define 3 levels in height 2500m, 1500m and 0m. Levels are stacked for each mesh element from from top to bottom. The height is absolute to datum.

LEVELS = [2500m, 1500m, 0m]

With 3 levels, we have 4 volumes, 2 volumes C1 and C2 above face F1 and two C3 and C4 above face F2

For each volume we measure air temperature in the volume's center during day (12:00) and night (24:00).

So we have in total (size(LEVELS)-1) * size(FACES) * size(TIMES) = (3-1) * 2 * 2 = 8 measured temperatures

Volume #	Time 12:00
C1 (2000m)	-2°C
C2 (750m)	10°C
C3 (2000m)	0°C
C4 (750m)	12°C

Volume #	Time 24:00
C1 (2000m)	-2°C
C2 (750m)	8 °C
C3 (2000m)	0°C
C4 (750m)	6°C

Changes in MDAL

We implement support for 3D Layered Meshes in MDAL. For formats that do not support meshes with vertical dimension, the API will be backward compatible. There will be new API on datasets to determine number of vertical levels for a given face on 2D base mesh frame. Also one would frequently need to determine first index of 3D volume from 2D face from base mesh frame (e.g. for identify tool to list vertical levels)

MDAL_EXPORT int MDAL_D_volumesCount( DatasetH dataset );

Also there will be 5 new enum options for retrieving data from MDAL_D_data

• For VERTICAL_LEVEL_COUNT_INTEGER, the returning array will be number of vertical levels for particular face, with size faceCount * size_of(int)
• For VERTICAL_LEVEL_DOUBLE, the returning array will be extrusion for given vertical level for particular face, with size (faceCount + volumesCount) * size_of(double)
• For FACE_INDEX_TO_VOLUME_INDEX_INTEGER, the returning array will be first 3D index for particular face, with size faceCount * size_of(int)
• For SCALAR_VOLUMES_DOUBLE, the returning array will be scalar values for given 2D base face for all stacked volumes above, with size volumesCount * size_of(double)
• For VECTOR_2D_VOLUMES_DOUBLE, the returning array will be vector values for given face for all stacked volumes, with size 2 * volumesCount * size_of(double)

Data for volumes will be ordered from bottom to top in resulting array.

//! Populates buffer with values from the dataset
//! for nodata, returned is numeric_limits<double>::quiet_NaN
//!
//! \param dataset handle to dataset
//! \param indexStart index of face/vertex to start reading of values to the buffer
//! \param count number of values to be written to the buffer
//! \param dataType type of values to be written to the buffer
//! \returns number of values written to buffer. If return value != count requested, see MDAL_LastStatus() for error type
MDAL_EXPORT int MDAL_D_data( DatasetH dataset, int indexStart, int count, MDAL_DataType dataType, void *buffer );

To return to our example data from previous section the returns will be (note that dataset already defines which time you selected):

MDAL_M_volumesCount(ds) == 4
MDAL_D_data(ds, 0, 2, VERTICAL_LEVEL_COUNT_INTEGER, ...) == [3, 3]
MDAL_D_data(ds, 0, 6, VERTICAL_LEVEL_DOUBLE, ...) == [2500, 1500.0, 0.0, 2500, 1500.0, 0.0]
MDAL_D_data(ds, 0, 2, FACE_INDEX_TO_VOLUME_INDEX_INTEGER, ...) == [0, 2]
MDAL_D_data(ds, 0, 4, SCALAR_VOLUMES_DOUBLE, ...) == [-2.0, 10.0, 0.0, 12.0]

Also we will introduce new option so data can be defined on vertices, faces and volumes.

Changes in QGIS

On QGIS side, the QgsMeshLayer dataset's group tree will recognise that the dataset group consists of various vertical levels and group the stacked levels in subgroup. The root of the subgroup can be selected as an active contour/vector group. User can also select one particular level in subgroup and in that case the rendering is done same way as for other datasets

To be able to visualise the subgroup with with various levels, we need to do averaging of the values in the levels to one single value to be displayed. We will use common methods described in this link. For example sigma method where you can calculate average with some levels excluded

and this would be corresponding QGIS widget

Affected Files

most of the mesh related files and MDAL library

Performance Implications

None

Further Considerations/Improvements

• The support of 3D Layered Mesh in QGIS 3D view is considered for future.
• Integrate the averaging (aggregation) methods to mesh layer calculator, so user would be able to persist the calculated datasets.

Backwards Compatibility

For all formats that do not support layered mesh, the API will be backward compatible.

Issue Tracking ID(s)

• https://github.com/lutraconsulting/MDAL/issues/154
• https://github.com/qgis/QGIS/issues/31574

[nyalldawson]

Sounds great to me!

[pcav]

Same here, cool stuff.

[DelazJ]

I don't fully understand what it's about but wanted to say it's great/pleasant to have handwritten items in the description

[timlinux]

Great stuff! Does the averaging take into account only the checked / enabled levels?

[PeterPetrik]

@timlinux there are multiple different averaging methods and we can add more on demand. I think basic methods just take average or you can exclude edge levels or so...

[wonder-sk]

@DelazJ hmm if you have problems understanding the QEP, then probably other people will have troubles as well. Anything we can improve in to proposal to make it easier to understand? Maybe expand the intro?

As a TL;DR summary, if you are familiar with 3D rasters, this QEP is about adding support for a similar thing for meshes, just more generic (extras: 1. data can change over time. 2. elevations can change over time).

[NathanW2]

This sounds pretty bloody great.

[saberraz]

For those interested in 3D raster, how GRASS GIS supports it: https://grass.osgeo.org/grass76/manuals/raster3dintro.html

It would be good also to support 3D rasters (if it is supported by GDAL) in addition to the MDAL unstructured mesh formats.

[andreasneumann]

This sounds pretty bloody great.

agreed - looking forward to it!

[PeterPetrik]

** update 22nd October 2019: Generalize the QEP to allow "The number of vertical levels can't change in time, but it can vary in spatial dimension"