Load and Style a Mesh Layer in PyQGIS

A mesh layer is what a numerical model actually produces: values on an unstructured grid of triangles or quadrilaterals, usually several variables, usually many timesteps. QGIS reads them through MDAL and renders them natively, and the whole thing is scriptable — but the addressing scheme is unlike anything else in the API, and the renderer settings follow a copy-modify-write pattern that silently does nothing if you get it wrong.

This recipe belongs to Temporal & 3D Visualization in PyQGIS. It covers loading, enumerating dataset groups and timesteps, sampling values at a point, and configuring the scalar and vector renderers.

One grid, several variables, many timestepsThe mesh geometry is fixed: a set of vertices joined into triangular or quadrilateral faces. On top of it sit dataset groups, each one a variable such as depth or velocity. Each group holds one dataset per timestep, and each dataset holds one value per vertex or per face.QgsMeshDatasetIndex(group, dataset) addresses everythingthe grid — fixedvertices and faces, never changedataset groupsgroup 0 · "Bed Elevation" · scalar · 1 datasetgroup 1 · "Water Depth" · scalar · 72 datasetsgroup 2 · "Velocity" · vector · 72 datasetsgroup indices are file order — always enumerate, never hard-codea re-run of the model with an extra output variable renumbers everything after it

Prerequisites

  • QGIS 3.34 LTR or newer. Mesh support arrived in 3.4 and matured through 3.14.
  • A mesh file — NetCDF, GRIB, UGRID, XMDF, 2DM, or another format MDAL reads.

Load and enumerate

from qgis.core import QgsMeshLayer, QgsMeshDatasetIndex, QgsProject

mesh = QgsMeshLayer("/data/model/flood.nc", "flood", "mdal")
if not mesh.isValid():
    raise SystemExit("MDAL could not read this file")

print(mesh.crs().authid(), mesh.extent().toString(0))

for group in range(mesh.datasetGroupCount()):
    index = QgsMeshDatasetIndex(group, 0)
    meta = mesh.datasetGroupMetadata(index)
    steps = mesh.datasetCount(index)
    kind = "vector" if meta.isVector() else "scalar"
    print(f"{group}: {meta.name()!r} {kind}, {steps} timestep(s)")

QgsProject.instance().addMapLayer(mesh)

Breakdown: The provider key is always "mdal" — the specific driver is chosen from the file, so there is no per-format constructor to remember. QgsMeshDatasetIndex(group, dataset) is the addressing pair used everywhere in this API; group metadata is the same for every timestep, so passing dataset 0 to fetch it is the convention. A group with one dataset is static — bed elevation, roughness, a mask — and one with many is a time series. Enumerating rather than hard-coding group numbers is not defensive programming, it is necessary: the indices are file order, and a model re-run with one extra output variable shifts every group after it.

Many meshes carry no CRS, in which case authid() comes back empty and the layer lands wherever the project puts it. Set it explicitly:

from qgis.core import QgsCoordinateReferenceSystem

if not mesh.crs().isValid():
    mesh.setCrs(QgsCoordinateReferenceSystem("EPSG:27700"))

Breakdown: Unlike a point cloud, where the fix belongs in the file, setting the CRS on a mesh layer is the normal remedy — most model formats have no place to store one, so QGIS expects you to declare it, and it is saved with the project. Get it wrong and the mesh is in the right shape in the wrong place, which is obvious against a basemap and invisible on its own.

Timesteps and reference time

group = 1
count = mesh.datasetCount(QgsMeshDatasetIndex(group, 0))

for step in range(min(count, 5)):
    meta = mesh.datasetMetadata(QgsMeshDatasetIndex(group, step))
    print(step, meta.time(), "hours after reference")

print("reference time:", mesh.temporalProperties().referenceTime())

Breakdown: meta.time() returns hours since the layer's reference time, not an absolute datetime — that is the crucial detail, and it is why setting the reference time correctly matters. Model files usually carry one; where they do not, QGIS defaults to something arbitrary and every timestep is offset by the difference. mesh.temporalProperties().setReferenceTime() fixes it, and doing so shifts the whole animation onto the real calendar so the mesh lines up with rainfall records, tide gauges or anything else time-aware in the project.

Sample a value

import math
from qgis.core import QgsPointXY

point = QgsPointXY(445120, 187330)
series = []
for step in range(count):
    value = mesh.datasetValue(QgsMeshDatasetIndex(group, step), point)
    scalar = value.scalar()
    series.append(None if math.isnan(scalar) else scalar)

print("peak depth:", max(v for v in series if v is not None))

Breakdown: datasetValue interpolates within the face containing the point, so it works anywhere inside the mesh rather than only at vertices. Outside the mesh it returns NaN, which is why the isnan check is required — NaN compares false against everything, so a max() over raw values silently returns a real number while a sum() returns NaN, and neither tells you the point was outside. For a vector group, value.x() and value.y() give the components and value.scalar() gives the magnitude. Sampling a whole series at one point like this is the fastest route from a model output to a chart, and it is a loop rather than a bulk call.

Style the scalar surface

Settings are copies until you put them backCalling rendererSettings returns a copy of the layer's settings. Modifying that copy changes nothing on the layer. Only calling setRendererSettings with the modified copy applies the change, and the same pattern applies one level down to the per-group scalar settings.Forget the last step and nothing happens, quietlyrendererSettings()returns a copyscalarSettings(g)a copy of a copymodify the rampstill no effectwrite both backnow it appliesboth writes are needed, innermost firstsettings.setScalarSettings(g, scalar)mesh.setRendererSettings(settings)the same pattern governs QgsRasterLayer renderers and layout item properties

from qgis.core import QgsColorRampShader, QgsStyle

depth_group = 1

settings = mesh.rendererSettings()
settings.setActiveScalarDatasetGroup(depth_group)

scalar = settings.scalarSettings(depth_group)
shader = scalar.colorRampShader()
shader.setColorRampType(QgsColorRampShader.Interpolated)
shader.setSourceColorRamp(QgsStyle.defaultStyle().colorRamp("Blues"))
shader.setMinimumValue(0.0)
shader.setMaximumValue(3.5)
shader.classifyColorRamp(classes=8, band=-1)
scalar.setColorRampShader(shader)

settings.setScalarSettings(depth_group, scalar)
mesh.setRendererSettings(settings)
mesh.triggerRepaint()

Breakdown: Both write-backs are required and in that order — modifying the scalar settings alone changes a temporary object. Setting the minimum and maximum explicitly rather than letting QGIS classify from the current timestep is what keeps colours comparable across an animation; classified per frame, a rising flood stays the same colour while the legend silently changes underneath it. A minimum of exactly zero on a depth variable renders dry ground as the ramp's lightest colour rather than as nothing, which is usually not what you want — setting the minimum slightly above zero and enabling the "below minimum" transparency is the conventional fix.

Style the vector field

vector_group = 2

settings = mesh.rendererSettings()
settings.setActiveVectorDatasetGroup(vector_group)

vector = settings.vectorSettings(vector_group)
vector.setLineWidth(0.6)
vector.setColor(QColor("#1e3a5f"))
vector.setShaftLengthMethod(0)          # 0 = min/max scaling
settings.setVectorSettings(vector_group, vector)
mesh.setRendererSettings(settings)

Breakdown: A scalar group and a vector group can be active at once, which is exactly the flood map convention: depth as a colour surface, velocity as arrows over it. The shaft length method decides how magnitude becomes arrow length — scaled between the data's minimum and maximum, fixed, or scaled by a factor — and min/max scaling is the readable default until the magnitudes span orders of magnitude, at which point a fixed length coloured by magnitude reads better. Arrow density is a separate setting and matters more than any of this: a mesh with fifty thousand faces drawn with one arrow each is a solid block of ink.

Finding the range across the whole run

Fixing the colour range needs the run's true minimum and maximum, and the group metadata carries them without a scan.

meta = mesh.datasetGroupMetadata(QgsMeshDatasetIndex(depth_group, 0))
print("group range:", meta.minimum(), "to", meta.maximum())

peaks = []
for step in range(count):
    step_meta = mesh.datasetMetadata(QgsMeshDatasetIndex(depth_group, step))
    peaks.append((step_meta.time(), step_meta.maximum()))

worst_time, worst_value = max(peaks, key=lambda pair: pair[1])
print(f"peak of {worst_value:.2f} at {worst_time:.1f} h")

Breakdown: Group metadata gives the extremes across every timestep, which is exactly what the colour ramp should be classified against. Per-dataset metadata gives them per timestep, which is how you find the peak of a flood without sampling anything — and knowing which hour that was tells you which frame to use as the still image in a report. Both are read from the file's own summary rather than computed, so they are instant even on a large run. A group whose reported maximum is absurd usually means the model wrote a fill value into the data rather than a nodata flag.

QGIS version compatibility

Mesh layers arrived in QGIS 3.4 and the renderer settings API reached its current shape by 3.14. QgsMeshLayer.datasetValue and the dataset index have been stable since. Temporal properties for mesh layers, including the reference time, arrived with the wider temporal framework in 3.14. The available MDAL drivers grow with each release, so a file that will not open on 3.34 is worth retrying on a newer build.

Troubleshooting

  • isValid() is False. MDAL has no driver for this file, or the file is a NetCDF whose variables MDAL cannot interpret as a mesh.
  • The mesh is in the wrong place. No CRS in the file. Set it on the layer.
  • Every sampled value is NaN. The point is outside the mesh, or the group is a vector group being read with .scalar() when the components are null.
  • Styling changes have no effect. A missing setRendererSettings or setScalarSettings write-back.
  • Timesteps are decades off. The reference time is wrong; meta.time() is relative, not absolute.
  • The map is solid with arrows. Reduce arrow density, or switch to streamlines for a dense field.

Conclusion

Enumerate the groups rather than assuming their numbers, set the CRS and the reference time before anything else, check for NaN on every sample, and remember that renderer settings are copies until written back. Mesh support in QGIS is genuinely good; the API just asks you to be explicit about things other layer types infer.

Frequently Asked Questions

Can I convert a mesh to a raster? Yes — native:meshrasterize grids a chosen dataset group and timestep into a raster, after which everything in raster analysis workflows applies. There are matching algorithms for exporting contours and vectors from a mesh.

How do I animate a mesh? Set the reference time, activate the layer's temporal properties, and drive the controller as usual — the layer picks the dataset matching each frame automatically.

Can I edit mesh values from PyQGIS? QGIS supports editing the mesh frame — adding and removing vertices and faces — but dataset values are read-only. Changing values means regenerating the file from the model.

Why do two dataset groups have the same name? Some formats store maximum and instantaneous variants under similar names. Print the group index alongside the name and select by index once you have identified the right one.