Day 2: Input Data Preparation Tasks

[pallavhydro]

Session Preparation

Download:

• Link: https://nc.ufz.de/s/6fyD6wC24MgKSm4
• Password: on the screen.

This includes:

• Environment installation script (MacOS, Ubuntu, WSL) and environment yaml file
• Jupyter Notebook
• Raw data

Environment and Jupyter

• Go to the installation files:

cd environment_installation/

• Execute the installation file, it will create an environment called mhm_day2 (change name if wanted)

bash installation_ubuntu_mac.sh mhm_day2 requirements.yml
cd ..

• Activate the environment

conda activate mhm_day2

Note: If the environment doesn't get recognized just by name, try using the full path instead. Full path can be found using conda env list

• Open Jupyter lab from the conda environment

jupyter lab

• Navigate and open the Jupyter notebook data_prep_mhm_workshop_2026.ipynb

---

First Block

Exercise 1 [Morphology] Prepare DEM, slope, aspect, flow direction, and flow accummulation

• Fill the raw DEM and visualize.
  • Can you see any difference between the raw and filled DEMs?
  • Save filled DEM as ASCII
• Explore the pyflowdir package to generate flow direction from the filled DEM
  • flow direction can have different nomenclature, are you familiar to one required by mHM?
  • Visualize the flow direction
  • Save the fdir as an ASCII in python
• Now, explore the pyflowdir package to generate flow accumulation
  • Visualize the flow accumulation
  • Can you make smaller rivers appear in the plot as well?
  • Save the fdir as an ASCII in python
• Use GDAL to generate slope from the filled DEM
  • Visualize the slope. Should the slope be in percentage rise or angle of elevation in mHM?
  • Save the slope as an ASCII
• Use GDAL to generate aspect from the filled DEM
  • Visualize the aspect.
  • Save the aspect as an ASCII

---

Exercise 2 [Morphology] Prepare the catchment shape and gauge location data

• Delineate both catchments using mhm-tools create_catchment. As input either use your filled dem or your fdir file.
• Visualize the mask or uparea_grid from the created basins_ids.nc file
• Optional: To see where the gauge is located you can load your facc.asc file you created in Ex1, plot it again and plot the gauge you get from idgauges.ascon top of it. If you want you can also overlapp the uparea_gridfrom the basin_ids.ncfile.
• Compare the two catchments visually.
• Use gdalinfo to find the latlon extend that covers both catchments.
• Extend both idgauges.asc files to cover the whole domain and merge them into one file.
• Mask the dem.asc file(s) using the mask from catchment delineation and expand it/them to the whole domain. Merge if necessary.
• Visuallize the masked dem.
• Crop all files you created in Exercise 1 to to the domain extend. You can use CDO, GDAL or mhm-tools crop_mhm_setup.

---

Exercise 3 [Morphology] Crop the rest of the morphological input files.

• Use mhm-tools crop_mhm_setupfor to clip/crop the remaining morphological data. Try to find a smart way.
  • For a help page of mhm-tools run mhm-tools -h
  • If you need information about a specific tool and it's input run mhm-tools <TOOL_NAME> -h.

---

Second Block

Exercise 4 [latlon file] Generate the relationship that connects model resolution to morphology

• mhm-tools can helps you to make latlon files.
  • Explore the help of the latlon command:

  `mhm-tools latlon -h`
  

  • Prepare the necessary arguments for mhm-tools latlon command
    • --level0 LEVEL0
    • --level1 LEVEL1
    • --level11 LEVEL11 (optional, but necessary when routing streamflow)
  • Run the mhm-tools latlon command

---

Exercise 5 [Meteorology] Prepare the meteorological forcings

• Prepare the meterological data such that mHM can use it. Steps that might or might not be needed:

  1. merge the files together into one
  2. crop to the domain you want
  3. create header file
  4. The variable names are hard coded in mHM so you might have to change them.
  5. Any data provided for mHM should be of the data type DOUBLE.
  6. Time must be in Integers. e.g hours since …, days since …
  7. Units must be correct
  8. No data value must be -9999.0 (if float)
  9. Latitude must be decreasing

• If you need inspiration: We provide a bash and a python workflow for precipitation data

---

Third Block

Exercise 6 [Observations] Prepare the observed discharge data

• Format raw time series observations to mHM format
  • Use langauge of your choice (python, bash, R, ...) to process the format
  • Raw data is available at data_raw/observation/grdc/
  • Format of mHM can be referred to from the test domain input data from Day 1

Hint: ChatGPT or AI companion of your choice comes in handy here!

---

Exercise 7 [mHM test run] Run mHM with the input data prepared from this session

• Prepare the mHM namelists for the setup
  • Download the nml folder from the nextcloud
  • Save the path to PATH_NMLS in the preable
  • Adjust the paths to your data in mhm.nml
• Run mHM to test the data as follows...
  • use only the 6340200 catchment with the mhm_parameter.nml file provided.
    • Are you happy with the model performance?
    • plot the hydrograph from mHM output discharge.nc
  • use only the 6340200 catchment with the mhm_parameter_calibrated.nml file provided.
    • Any changes in the model performance?
    • plot the hydrograph from mHM output discharge.nc
  • use both the 6340200 and 6340220 catchments with the mhm_parameter_calibrated.nml file provided.
    • What can be done to improve the small catchment performance?
  • use both the 6340200 and 6340220 catchments this time with the resolution of 0.0625 degree.
    • Did this help the model performance? If so, what could be the reason?

---

Optional Exercise 8

• Nested simulations (Idea from Bibek + Stephan)
  • Run only the small catchment 6340220 using its corresponding DEM mask at 0.0625 degree.
  • Format the output time series at 6340220 to mHM gauge input time series format.
  • Run the larger catchment 6340200 using 6340220 as inflow gauge.
    • What conclusions do you make?

---

Optional Exercise 9

• Repeat the above experiment but calibrate 6340220 at 0.0625 degree resolution.
  • Transfer the parameter to larger catchment 6340200.
    • What conclusions do you make?

[MuellerSeb]

A good workaround for my_path in the beginning is:

my_path = Path("..").resolve()
print(my_path)

Assuming you went into the jupyter_notebook/ folder to start the notebook.

[silentin]

Solution Exercise1:

Paths:

# Set paths as environment so that it is available to all cells in this notebook

my_path = Path("..").resolve() # if you opened jupyter lab from the path where you also downloaded the data to 
os.environ["PATH_DATA_SRC"]   = os.path.join(my_path, "data/data_raw/")
os.environ["PATH_DATA_TGT"]   = os.path.join(my_path, "data/data_processed/")
os.environ["PATH_WORK"]       = os.path.join(my_path, "data/work_dir/")
os.environ["PATH_NMLS"]       = os.path.join(my_path, "nmls/")
os.environ["PATH_MHM_OUTPUT"] = os.path.join(my_path, "mhm_output/")
os.environ["PATH_SCRIPTS"]    = os.path.join(my_path, "scripts/")

If there is %%bash at top of a cell it is a bash cell. You can then call the names you set in the beginning with $PATH_DATA_TGT.
If you are in a python cell you need to call them using: os.environ["PATH_DATA_TGT"].

Solution Exercise 1.1

# Plot and visualize the raw DEM file

# Path to DEM
dem_raw_file = os.path.join(os.environ["PATH_DATA_SRC"], "static/morph/dem.asc")

# Load DEM as xarray
dem_xr = rxr.open_rasterio(dem_raw_file, masked=True).squeeze(drop=True)

# Visualize the DEM
dem_xr.plot()
plt.show()
# ====== FILL THE DEM FILE =========

# You can easily fill a DEM file using the pyflwdir command as below:

# Convert to NumPy array
dem_np = dem_xr.values.squeeze()  # 2D numpy array

# Fill depressions using pyflowdir
dem_filled_np, filled_mask = pfd.dem.fill_depressions(dem_np)

# convert back to xarray which preserves metadata
dem_filled_xr = xr.DataArray(dem_filled_np,
                           coords=dem_xr.coords,
                           dims=dem_xr.dims,
                           attrs=dem_xr.attrs)

# Visualize the filled DEM
dem_filled_xr.plot()
plt.show()

Solution Exercise 1.2

# >>>>>> Exercise 1.2 <<<<<<<<

# Visualize the difference between the raw and filled DEM

# Diff the two DEMs
dem_diff_xr = dem_filled_xr - dem_xr


dem_diff_xr.plot()
plt.show()
# Include nodata attribute
dem_filled_xr.rio.write_nodata(-9999, inplace=True)

# The DEM is currently still in the raw extent and needs to be clipped later. 
# For now save these intermediate files to PATH_WORK

# Save as ASCII
filled_dem_path = Path(os.environ["PATH_WORK"]) / "full_extent/static/morph/dem_filled.asc"
filled_dem_path.parent.mkdir(exist_ok=True, parents=True)
dem_filled_xr.rio.to_raster(filled_dem_path, driver="AAIGrid")

Solution Exercise 1.3:

# Compute flow direction (D8 nomenclature for mHM)
fdir_pfd = pfd.from_dem(dem_filled_np) # using np and not xr object as the function needs numpy array as input

# convert the pyflwdir object to numpy array
fdir_np = fdir_pfd.to_array()

# convert back to xarray which preserves metadata
fdir_xr = xr.DataArray(fdir_np,
                           coords=dem_xr.coords,
                           dims=dem_xr.dims,
                           attrs=dem_xr.attrs)

# Visualize the flow direction
fdir_xr.plot()
plt.show()

# Include nodata attribute
fdir_xr.rio.write_nodata(0, inplace=True) # 1 to 128 are the valid directions. Storing -9999 would require signed integer data type. So, 0 for nodata. 


# Save as ASCII
fdir_xr.rio.to_raster(os.path.join(os.environ["PATH_WORK"], "full_extent/static/morph/fdir.asc"), driver="AAIGrid")

Solution Exercise 1.4:

# Look for pyflwdir command to generate facc.asc.

# Compute flow accumulation
facc_pfd = fdir_pfd.upstream_area(unit = 'cell') # needs FlwdirRaster

# convert back to xarray which preserves metadata
facc_xr = xr.DataArray(facc_pfd,
                           coords=dem_xr.coords,
                           dims=dem_xr.dims,
                           attrs=dem_xr.attrs)

# Visualize the flow direction
facc_xr.plot()
plt.show()
# Include nodata attribute
facc_xr.rio.write_nodata(-9999, inplace=True)


# Save as ASCII
facc_xr.rio.to_raster(os.path.join(os.environ["PATH_WORK"], "full_extent/static/morph/facc.asc"), driver="AAIGrid")

Solution Exercise 1.5
bash cell:

gdaldem slope -p -s 111120 ${PATH_WORK}/full_extent/static/morph/dem_filled.asc ${PATH_WORK}/full_extent/static/morph/slope.asc && echo "slope file created and saved!"
# visualize the slope

python cell:

# Path to slope file
slope_file = os.path.join(os.environ["PATH_WORK"], "full_extent/static/morph/slope.asc")


# Load DEM as xarray
slope_xr = rxr.open_rasterio(slope_file, masked=True).squeeze(drop=True)

# plot
slope_xr.plot()
plt.show()

Solution Exercise 1.6:
bash cell:

gdaldem aspect ${PATH_WORK}/full_extent/static/morph/dem_filled.asc ${PATH_WORK}/full_extent/static/morph/aspect.asc && echo "aspect file created and saved!"
# visualize the aspect

python cell:

# Path to aspect file
aspect_file = os.path.join(os.environ["PATH_WORK"], "full_extent/static/morph/aspect.asc")

# Load DEM as xarray
aspect_xr = rxr.open_rasterio(aspect_file, masked=True).squeeze(drop=True)

# plot
aspect_xr.plot()
plt.show()

[silentin]

Extend you need to set in a python cell:

# Set the extents

%env XMIN = 10.125
%env XMAX = 11.875
%env YMIN = 50.625
%env YMAX = 51.875

[silentin]

Description of the era5 grid: https://confluence.ecmwf.int/display/CKB/ERA5%3A+What+is+the+spatial+reference

[Jetsonic]

Due to some exception some of use were having an error where the extent of the merged idgauge file changed when doing with gdalwarp for the work around suggest to use a python script shared here (Exercise 2.9) :

import rasterio
from rasterio.merge import merge

def merge_rasters_preserve_extent(input_folder, output_file, nodata_value=-9999,reference_file=None):
  """
  Merge all raster files (.tif) from a folder into one output raster.
  Skips the output file if it already exists in the folder.

  Parameters
  ----------
  input_folder : str or Path
      Folder containing raster files to merge.
  output_file : str
      Path to the merged raster output.
  nodata_value : int or float, optional
      NoData value for input and output rasters. Default is -9999.
  reference_file : str or Path, optional
      Reference raster to align extent/resolution to. Default is None.

  Returns
  -------
  str
      Path to the merged raster file.
  """
  
  input_folder = Path(input_folder)
  output_file = Path(output_file)

  # Find all .tif files but exclude the output file if present
  tif_files = sorted([f for f in input_folder.glob("*.tif") if f.resolve() != output_file.resolve()])
  
  if not tif_files:
      raise ValueError(f"No .tif files found to merge in folder: {input_folder}")

  print(f"Found {len(tif_files)} .tif files to merge (output file skipped if present).")

  # Open all source rasters
  srcs = [rasterio.open(f) for f in tif_files]

  if reference_file:
      with rasterio.open(reference_file) as ref:
          out_transform = ref.transform
          out_crs = ref.crs
          out_shape = (ref.height, ref.width)
      print(f"Aligning to reference grid: {reference_file}")
      merged_array, merged_transform = merge(
          srcs,
          bounds=ref.bounds,
          res=ref.res,
          nodata=nodata_value,
          method='first'
      )
  else:
      merged_array, merged_transform = merge(srcs, nodata=nodata_value, method='first')
      out_crs = srcs[0].crs

  # Prepare metadata
  out_meta = srcs[0].meta.copy()
  out_meta.update({
      "driver": "GTiff",
      "height": merged_array.shape[1],
      "width": merged_array.shape[2],
      "transform": merged_transform,
      "crs": out_crs,
      "nodata": nodata_value,
      "dtype": merged_array.dtype
  })

  # Write merged raster
  with rasterio.open(output_file, "w", **out_meta) as dest:
      dest.write(merged_array)

  for src in srcs:
      src.close()

  print(f"All idgauges merged and  raster saved as: {output_file}")
  return str(output_file)

merge_rasters_preserve_extent(
  input_folder=Path(os.environ["PATH_WORK"]) / "process",
  output_file=Path(os.environ["PATH_WORK"]) / "process" / "idgauges.tif",
  nodata_value=-9999,
  reference_file=None
)

[MuellerSeb]

This is great. Thanks for sharing! 🎉

[elojonen]

Thanks! However, it is important to double check if there are any other tifs in the directory. If yes, it will also merge these and the output gets corrupted. This can be fixed by filtering the tif_files list one level further for only files that e.g. contain "idgauges" in their name:

tif_files = sorted([
    f for f in input_folder.glob("*.tif")
    if f.resolve() != output_file.resolve()
    and "idgauges" in f.name.lower()
])