Deltares / imod-qgis

🗺️🧭 QGIS plugin for iMOD
GNU General Public License v2.0
7 stars 1 forks source link

How can we use the QGIS-plugin for refined grids (disv) #70

Closed rubencalje closed 5 months ago

rubencalje commented 5 months ago

Hi developers,

The QGIS-plugin looks great! I am interested in using the plugin to plot cross-sections and time-series graphs in QGIS, for models we generate with our Python-package nlmod, see https://github.com/gwmod/nlmod. This package has some similarities to iMod-Python, as we both use xarray Datasets to store model data.

The link on the main page of your GitHUB (https://github.com/Deltares/imod-qgis), that should contain an example of preparing such a dataset in python is broken: https://deltares.github.io/iMOD-Documentation/workflow_wq.html#convert-output-data.

So my question is: how can we best generate input for the QGIS-plugin from a xarray-Dataset, also for models that use grid-refinement (disv-package)?

Thank you for your time!

Huite commented 5 months ago

Hi @rubencalje,

For structured data (DIS), we use these two utilities to convert to a UGRID netCDF format, which QGIS can read and treats as a mesh:

https://deltares.github.io/imod-python/api/generated/util/imod.util.to_ugrid2d.html#imod.util.to_ugrid2d https://deltares.github.io/imod-python/api/generated/util/imod.util.mdal_compliant_ugrid2d.html#imod.util.mdal_compliant_ugrid2d

The reason is structured netCDF data doesn't work so well in QGIS due to GDAL, since it treats layer and time as a band. The Mesh Layers actually have support for the time axis in QGIS, which makes it easier to work with. I've been thinking about a "netCDF manager" in QGIS to make working with multidimensional structured netCDFs a bit easier, but I never got past a very rudimentary prototype.

For DISV grids, we use xugrid: https://github.com/deltares/xugrid

Anyway, concretely, if you have some output you want to visualize in QGIS:

import imod

head = imod.mf6.open_hds("results.hds", "disv.grb")  # returns a xugrid.UgridDataset
head_dataset = head.ugrid.to_dataset()  # convert to ordinary xarray.Dataset
mdal_dataset = imod.util.to_mdal_compliant_ugrid2d(head_dataset)
mdal_dataset.to_netcdf("ready-for-qgis.nc")

There's some logic in the open_hds which determines whether to return an xarray Datase (DIS) or an xugrid UgridDataset (DISV). The xugrid Dataset has a separate grid object inside of it, storing the unstructured topology. This page hopefully gives a bit of an overview: https://deltares.github.io/xugrid/examples/quick_overview.html

With .ugrid.to_dataset() this data is converted into xarray variables, according to the UGRID conventions: https://ugrid-conventions.github.io/ugrid-conventions/

Finally, this works well for variables with a single layer, not for multi-layer variables. So we end up using the to_mdal_compliant_ugrid2d since it unstacks the layers into separate variables. This isn't ideal, it's on our to-do list to expand QGIS and MDAL support for layered UGRID netCDFs.

If you want to visualize input and the like, you will need to construct a xugrid.Ugrid2d topology and join it with your dataset into a UgridDataset or UgridDataArray. That should be straightforward, you essentially need the data from this block of DISV input:

BEGIN CELL2D
<icell2d> <xc> <yc> <ncvert> <icvert(ncvert)>
<icell2d> <xc> <yc> <ncvert> <icvert(ncvert)>
...
END CELL2D

@JoerivanEngelen should probably take a look at why the link doesn't work.

rubencalje commented 5 months ago

Thanks for the quick response, and for sending me in the right direction!

JoerivanEngelen commented 5 months ago

I fixed the broken link in the documentation. I think this makes a good Q&A discussion, as I think more people are probably wondering how to prepare their data, so I'm converting this to a discussion.