Convert Wflow staticmaps netcdf to raster files#

In order to inspect or (manually) modify Wflow staticmaps it is convenient to export the maps to a raster format. Here we show how to read the model maps and save to a so-called mapstack (i.e.: a set of raster files with identical grid) using HydroMT.

Load dependencies#

[1]:

import xarray as xr
import shutil
from os.path import join, exists
from os import listdir, rename
import hydromt
from hydromt_wflow import WflowModel

Read wflow staticmaps#

HydroMT provides an easy method to read the model schematization through the Model API.

[2]:

root = "wflow_piave_subbasin"
mod = WflowModel(root, mode="r")
ds = mod.grid  # here the staticmaps netcdf is loaded
print(ds)

<xarray.Dataset> Size: 1MB
Dimensions:               (latitude: 53, longitude: 58, time: 12, layer: 4)
Coordinates:
  * latitude              (latitude) float64 424B 46.68 46.67 ... 45.83 45.82
  * longitude             (longitude) float64 464B 11.78 11.8 ... 12.72 12.73
  * time                  (time) int32 48B 1 2 3 4 5 6 7 8 9 10 11 12
  * layer                 (layer) int32 16B 0 1 2 3
    spatial_ref           int32 4B 0
Data variables: (12/77)
    x_out                 (latitude, longitude) float64 25kB nan nan ... nan nan
    y_out                 (latitude, longitude) float64 25kB nan nan ... nan nan
    idx_out               (latitude, longitude) int32 12kB -1 -1 -1 ... -1 -1 -1
    wflow_ldd             (latitude, longitude) uint8 3kB 255 255 ... 255 255
    wflow_subcatch        (latitude, longitude) int32 12kB 0 0 0 0 0 ... 0 0 0 0
    wflow_uparea          (latitude, longitude) float32 12kB -9.999e+03 ... -...
    ...                    ...
    M_original            (latitude, longitude) float32 12kB -9.999e+03 ... -...
    M                     (latitude, longitude) float32 12kB -9.999e+03 ... -...
    f                     (latitude, longitude) float32 12kB -9.999e+03 ... -...
    wflow_soil            (latitude, longitude) int32 12kB -9999 -9999 ... -9999
    wflow_gauges          (latitude, longitude) int32 12kB 0 0 0 0 0 ... 0 0 0 0
    wflow_gauges_grdc     (latitude, longitude) int32 12kB 0 0 0 0 0 ... 0 0 0 0

Write netcdf to mapstack#

The raster module provides many raster GIS methods throught the raster Dataset accessor. To write a Dataset to a mapstack one line with code is sufficient. We only need to provide the output folder in which all raster files are saved. The default output format is GeoTIFF, but this can be changed with the driver argument. To write to PCRaster map-files it is recommended to have PCRaster python installed.

[3]:

updated_staticmaps = "updated_staticmaps"
if exists(updated_staticmaps):
    shutil.rmtree(updated_staticmaps)

ds.raster.to_mapstack(updated_staticmaps)
listdir(updated_staticmaps)

[3]:

['M_original_.tif',
 'wflow_reservoirareas.tif',
 'wflow_riverwidth.tif',
 'RiverSlope.tif',
 'wflow_landuse.tif',
 'WaterFrac.tif',
 'f.tif',
 'thetaS.tif',
 'dem_subgrid.tif',
 'kc.tif',
 'KsatVer_200.0cm.tif',
 'SoilMinThickness.tif',
 'wflow_soil.tif',
 'thetaR.tif',
 'LakeStorFunc.tif',
 'idx_out.tif',
 'h3_low.tif',
 'c.tif',
 'ResSimpleArea.tif',
 'M_.tif',
 'wflow_subcatch.tif',
 'KsatVer_100.0cm.tif',
 'LakeThreshold.tif',
 'KsatVer_30.0cm.tif',
 'wflow_lakeareas.tif',
 'Sl.tif',
 'RiverDepth.tif',
 'y_out.tif',
 'Lake_b.tif',
 'wflow_reservoirlocs.tif',
 'SoilThickness.tif',
 'wflow_ldd.tif',
 'wflow_river.tif',
 'wflow_glacierfrac.tif',
 'alpha_h1.tif',
 'KsatVer.tif',
 'KsatVer_0.0cm.tif',
 'LinkedLakeLocs.tif',
 'glacareas.tif',
 'wflow_gauges_grdc.tif',
 'wflow_streamorder.tif',
 'ResTargetFullFrac.tif',
 'ResTargetMinFrac.tif',
 'Slope.tif',
 'LAI.tif',
 'ResDemand.tif',
 'wflow_gauges.tif',
 'Swood.tif',
 'subare.tif',
 'h4.tif',
 'wflow_lakelocs.tif',
 'wflow_dem.tif',
 'h1.tif',
 'M.tif',
 'LakeArea.tif',
 'RootingDepth.tif',
 'LakeAvgOut.tif',
 'LakeAvgLevel.tif',
 'h3_high.tif',
 'wflow_uparea.tif',
 'ResMaxVolume.tif',
 'ResMaxRelease.tif',
 'h2.tif',
 'PathFrac.tif',
 'Lake_e.tif',
 'x_out.tif',
 'N_River.tif',
 'f_.tif',
 'M_original.tif',
 'Kext.tif',
 'KsatVer_60.0cm.tif',
 'KsatVer_5.0cm.tif',
 'wflow_glacierstore.tif',
 'LakeOutflowFunc.tif',
 'KsatVer_15.0cm.tif',
 'N.tif',
 'wflow_riverlength.tif']

Now the model files can easily be inspected and modified e.g. QGIS.

Note: in QGIS, you can also visualize netcdf files but direct modification is not (yet) possible.

Create staticmaps netcdf files based on mapstack#

If you want to update the staticmaps after modification the maps can be read into a Dataset by hydromt. We recommend the following workflow:

read the original model
read the updated mapstack
change the model root to write the updated model to a new directory
update the staticmaps of the model
write the model

Note: We do not read the forcing as it is faster to just copy the file instead of loading it into python and writing it back to netcdf.

Note: The staticgeoms might be changed because of manual changes in the wflow_river, lakes, reservoir or glacier staticmaps and are therefore not read here. To change these maps we recommend using the hydromt update method to keep the staticgeoms and maps aligned.

[4]:

# Read the model without the updated staticmaps
root = "wflow_piave_subbasin"
mod = WflowModel(root, mode="r")
mod.read_grid()
mod.read_config()

Adding updated maps#

Given some GeoTiff files, there are 2 ways to add their contents to your model: with or without renaming of the files.

The simplest is to overwrite the entire grid by reading all .tif files produced using the steps above.

Note: If any maps have been manually updated (using QGIS), do not rename the updated files. WflowModel uses the filenames to generate variable names automatically. Renaming files can result in a WflowModel with staticmaps containing both the old and new versions of a map, which may lead to inconsistencies in the model.

[5]:

# read the entire updated mapstack
ds_updated = hydromt.open_mfraster(join(updated_staticmaps, "*.tif"))

# Reinitialize and update the new model staticmaps
mod._grid = xr.Dataset()
mod.set_grid(ds_updated)

# Change root to not overwrite the original model
updated_root = "wflow_piave_subbasin_updated"
mod.set_root(updated_root)
mod.write()

Update static maps separately using `setup_grid_from_raster`#

Lets say you have updated and renamed the KsatVer.tif static map to KsatVer_updated.tif, and that you would like to include it in your model.

You will need to find the wflow variable that corresponds to the map you updated, which in the case of KsatVer.tif is soil_surface_water__vertical_saturated_hydraulic_conductivity. For information on the available wflow variables, see the Wflow documentation

[6]:

mod = WflowModel(root, mode="r")
mod.read_grid()
mod.read_config()

# `KsatVer` corresponds to the `soil_surface_water__vertical_saturated_hydraulic_conductivity` variable
wflow_variable = "soil_surface_water__vertical_saturated_hydraulic_conductivity"
print(f"Before setup_grid_from_raster: {mod._config['input']['static'][wflow_variable]}")

# For this example we just copy and rename `KsatVer.tif`, but in practice this would be a different file
old_map = join(updated_staticmaps, "KsatVer.tif")
updated_map = join(updated_staticmaps, "KsatVer_updated.tif")
shutil.copyfile(old_map, updated_map)

mod.setup_grid_from_raster(raster_fn=updated_map, reproject_method='nearest', wflow_variables=[wflow_variable])
print(f"After setup_grid_from_raster: {mod._config['input']['static'][wflow_variable]}")

Before setup_grid_from_raster: KsatVer
After setup_grid_from_raster: KsatVer_updated