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