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 shutil
from os.path import join, exists
from os import listdir
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"
wflow = WflowModel(root, mode="r")
ds = wflow.grid  # here the staticmaps netcdf is loaded
print(ds)

<xarray.Dataset> Size: 1MB
Dimensions:                               (latitude: 53, longitude: 58,
                                           time: 12, layer: 4)
Coordinates:
  * time                                  (time) int32 48B 1 2 3 4 ... 10 11 12
  * layer                                 (layer) int64 32B 0 1 2 3
    spatial_ref                           int64 8B 0
  * longitude                             (longitude) float64 464B 11.78 ... ...
  * latitude                              (latitude) float64 424B 46.68 ... 4...
Data variables: (12/72)
    meta_subgrid_outlet_x                 (latitude, longitude) float64 25kB ...
    meta_subgrid_outlet_y                 (latitude, longitude) float64 25kB ...
    meta_subgrid_outlet_idx               (latitude, longitude) int32 12kB -1...
    local_drain_direction                 (latitude, longitude) uint8 3kB 255...
    subcatchment                          (latitude, longitude) int32 12kB 0 ...
    meta_upstream_area                    (latitude, longitude) float32 12kB ...
    ...                                    ...
    meta_soilgrids_ksat_vertical_200.0cm  (latitude, longitude) float32 12kB ...
    soil_f_                               (latitude, longitude) float32 12kB ...
    soil_f                                (latitude, longitude) float32 12kB ...
    meta_soil_texture                     (latitude, longitude) int32 12kB -9...
    outlets                               (latitude, longitude) int32 12kB 0 ...
    gauges_grdc                           (latitude, longitude) int32 12kB 0 ...

Write netcdf to mapstack#

The raster module provides many raster GIS methods through the raster Dataset accessor. To write a Dataset into several raster files (eg one geotiff file per map), one line with code is sufficient using raster.to_mapstack. 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.

[3]:

# Name of the folder where the tif files will be saved
updated_staticmaps = "updated_staticmaps"

# First remove the folder if it already exists
if exists(updated_staticmaps):
    shutil.rmtree(updated_staticmaps)

# Save the tif files
ds.raster.to_mapstack(updated_staticmaps)
# Print the content of the folder
listdir(updated_staticmaps)

[3]:

['vegetation_feddes_alpha_h1.tif',
 'land_water_fraction.tif',
 'reservoir_max_volume.tif',
 'meta_soilgrids_ksat_vertical_15.0cm.tif',
 'gauges_grdc.tif',
 'meta_subgrid_outlet_y.tif',
 'soil_f.tif',
 'meta_soil_texture.tif',
 'reservoir_area.tif',
 'lake_storage_curve.tif',
 'river_depth.tif',
 'meta_soilgrids_ksat_vertical_200.0cm.tif',
 'land_elevation.tif',
 'subcatchment.tif',
 'local_drain_direction.tif',
 'lake_initial_depth.tif',
 'glacier_fraction.tif',
 'reservoir_target_min_fraction.tif',
 'lake_lower_id.tif',
 'river_width.tif',
 'vegetation_feddes_h1.tif',
 'meta_upstream_area.tif',
 'reservoir_target_full_fraction.tif',
 'soil_compacted_fraction.tif',
 'lake_area_id.tif',
 'meta_subgrid_outlet_idx.tif',
 'river_mask.tif',
 'vegetation_crop_factor.tif',
 'vegetation_feddes_h2.tif',
 'reservoir_demand.tif',
 'meta_soilgrids_ksat_vertical_5.0cm.tif',
 'river_manning_n.tif',
 'lake_outflow_threshold.tif',
 'reservoir_max_release.tif',
 'meta_subgrid_outlet_x.tif',
 'lake_b.tif',
 'lake_e.tif',
 'meta_soilgrids_ksat_vertical_30.0cm.tif',
 'lake_rating_curve.tif',
 'river_length.tif',
 'vegetation_kext.tif',
 'vegetation_feddes_h3_high.tif',
 'vegetation_feddes_h3_low.tif',
 'lake_area.tif',
 'meta_streamorder.tif',
 'soil_brooks_corey_c.tif',
 'vegetation_feddes_h4.tif',
 'soil_ksat_vertical.tif',
 'vegetation_leaf_area_index.tif',
 'soil_f_.tif',
 'reservoir_area_id.tif',
 'meta_soilgrids_ksat_vertical_100.0cm.tif',
 'soil_theta_r.tif',
 'soil_thickness.tif',
 'vegetation_leaf_storage.tif',
 'land_slope.tif',
 'river_slope.tif',
 'meta_soilgrids_ksat_vertical_0.0cm.tif',
 'meta_soilgrids_ksat_vertical_60.0cm.tif',
 'land_manning_n.tif',
 'outlets.tif',
 'meta_glacier_area_id.tif',
 'soil_theta_s.tif',
 'vegetation_root_depth.tif',
 'meta_landuse.tif',
 'reservoir_outlet_id.tif',
 'vegetation_wood_storage.tif',
 'glacier_initial_leq_depth.tif',
 'meta_lake_mean_outflow.tif',
 'meta_subgrid_elevation.tif',
 'lake_outlet_id.tif',
 'meta_subgrid_area.tif']

Now the model files can easily be inspected and modified for example using QGIS (e.g. with the Serval plugin).

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

Update a specific map using `setup_grid_from_raster`#

Let’s say you have updated (manually or not) the soil_ksat_vertical.tif static map and that you would like to include it in your model. For demonstration purpose, we will make a copy and name it soil_ksat_vertical_updated.tif.

[4]:

old_map = join(updated_staticmaps, "soil_ksat_vertical.tif")
updated_map = join(updated_staticmaps, "soil_ksat_vertical_updated.tif")
shutil.copyfile(old_map, updated_map)

[4]:

'updated_staticmaps/soil_ksat_vertical_updated.tif'

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

And then you can simply use setup_grid_from_raster.

In the below example, we will use python to update our model. If you wish to update via command line, the steps are:

Create a data catalog entry for the maps in the “updated_staticmaps” folder
Prepare a hydromt config file with the steps: setup_grid_from_raster, write_grid and write_config
Call the hydromt update command line re-using the catalog and configuration file you created.

[5]:

# Instantiate a new WflowModel object in read/write mode (r+) to be able to update it
wflow = WflowModel(root, mode="r+")

# `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: {wflow._config['input']['static'][wflow_variable]}")

# Path to the (updated) `soil_ksat_vertical_updated.tif`
updated_map = join(updated_staticmaps, "soil_ksat_vertical_updated.tif")

# Update the model KsatVer map with setup_grid_from_raster
wflow.setup_grid_from_raster(
    raster_fn=updated_map,
    reproject_method='nearest',
    variables=["soil_ksat_vertical_updated"],
    wflow_variables=[wflow_variable],

)
print(f"After setup_grid_from_raster: {wflow._config['input']['static'][wflow_variable]}")

# Write the updated grid and config files to new files
wflow.write_grid(fn_out="staticmaps_updated.nc")
wflow.write_config(config_name="wflow_sbm_updated.toml")

Before setup_grid_from_raster: soil_ksat_vertical
After setup_grid_from_raster: soil_ksat_vertical_updated