Using hydromt.flw.d8_from_dem by using rasterio data

[AquaPore]

Dear Wflow users,

I am trying to use hydromt to process D8 data by using my own DEM data abstracted from rasterio instead of using the one provided by hydromt.

I could off couse use Pyflwdir, but I believe that the advantage of using hydromt.flw.d8_from_dem is that it can use the riv_burn_method for which the pyflwdir.from_dem does not have this functionality.

Below is a code that I am attempting to develop which does not work, and any help to improve it will be welcomed.

Many thanks for any help you may provide,
Joseph

Code that does not works by using hydromt

def DEM_2_FLOWDIRECTION(Path_Input, Path_Output, Max_Depth=-1.0, Outlets="edge"):
   
   # Read elevation data using rasterio
   with rasterio.open(Path_Input, "r") as src:
      elevtn = src.read(1)
      nodata = src.nodata
      transform = src.transform
      crs = src.crs
      latlon = src.crs.is_geographic
      prof = src.profile
      
      Dem = hydromt.io.open_raster(elevtn)

   FlowDirection_Pyflwdir = hydromt.flw.d8_from_dem(Dem ,max_depth=-1, outlets="edge", idxs_pit=None)
   
   FlowDirection_Array = FlowDirection_Pyflwdir.to_array(ftype="ldd")
      
   return FlowDirection_Array, FlowDirection_Pyflwdir

Code that works but by using pyflwdir.from_dem

import rasterio
import numpy as np
import pyflwdir

def DEM_2_FLOWDIRECTION2(Path_Input, Path_Output, Max_Depth=-100.0, Connectivity=8):
   
	# Read elevation data using rasterio
	with rasterio.open(Path_Input, "r") as src:
		elevtn = src.read(1)
		nodata = src.nodata
		transform = src.transform
		crs = src.crs
		#  extent = np.array(src.bounds)[[0, 2, 1, 3]]
		latlon = src.crs.is_geographic
		prof = src.profile

	# returns FlwDirRaster object
	FlowDirection_Pyflwdir = pyflwdir.from_dem(data=elevtn, nodata=src.nodata, transform=transform, latlon=latlon, outlets="min", max_depth=Max_Depth, connectivity=Connectivity)

	FlowDirection_Array = FlowDirection_Pyflwdir.to_array(ftype="ldd")

	# Write to tiff file
	prof.update(dtype=FlowDirection_Array.dtype, nodata=-999)
	with rasterio.open(Path_Output, "w", **prof) as src:
		src.write(FlowDirection_Array, 1)

	return FlowDirection_Array, FlowDirection_Pyflwdir

[shartgring]

Thank you for your question @AquaPore. hydromt.io.open_raster takes the path to the elevation data as an argument. You are now passing the raw elevation data itself. hydromt.flw.d8_from_dem requires the elevation data as a DataArray object (from xarray), which is ensured by hydromt.io.open_raster.

Could you try if this works:

Dem = hydromt.io.open_raster(Path_Input)
FlowDirection_Pyflwdir = hydromt.flw.d8_from_dem(Dem ,max_depth=-1, outlets="edge", idxs_pit=None)

As a bit of background info: the main difference between the hydromt function and the pyflwdir functions is that pyflwdir uses numpy arrays, and HydroMT wraps the DataArray from xarray around that using the functionality from pyflwdir

[shartgring]

Is the goal to create a flow direction map to build a model using hydromt_wflow? If that is the case, then working with d8 will suffice. As the first step in model building (setup_basemaps()) can handle flow direction maps in both d8 and ldd format. hydromt_wflow will recognize that it is in d8 format and accordingly transform it into ldd format using core_conversion.d8_to_ldd.

Otherwise, this may also work:

import xarray
from pyflwdir.core_conversion import d8_to_ldd

Flwdir_D8 = hydromt.flw.d8_from_dem(Dem ,max_depth=Max_Depth, outlets=Outlets, idxs_pit=None)
Flwdir_Ldd = xarray.apply_ufunc(d8_to_ldd, Flwdir_D8.load())

As d8_to_ldd is a function that works with numpy arrays, we can use xarray.apply_ufunc (https://docs.xarray.dev/en/latest/generated/[xarray.apply_ufunc](https://docs.xarray.dev/en/latest/generated/xarray.apply_ufunc.html).html) with the function as an argument to perform it over the values of the DataArray

[shartgring]

Note: This quick solution does require loading the entire flow direction map into memory (Flwdir_D8.load()) which can be a problem for large maps. The implementation inside hydromt_wflow works differently and makes use of options in xarray to perform computations in a chucked way (so loading only chuncks of the maps into memory to perform calculations instead of the whole map)

[AquaPore]

A great thanks for your amazing kindness. I am preparing files to run WFLOW model for small catchements in Ireland, and we have our own lidar and river data, so therefore I am not using hdromt data to parameterize the model. I am therefore writting my own code with Python/ Julia to prepare the files for WFLOW model. My code will be made publically available.

Your solution is very clean by using xarray.apply_ufunc and the output is an xarray which is amazing..

Nevertheless, I got the message, nodata value missing I was not sucesfull to trie to put a nondata value. In Julia we can use NaN but I do not think in Python. How to put a nondata value and what value would you recommend?

"""

import hydromt
import pyflwdir
import xarray

def DEM_2_LDD_HYDROMT(Path_Input, Path_Output, Max_Depth=-1.0, Outlets="edge"):

# Read elevation data using hydromt which will transform the data into xarray

Dem = hydromt.io.open_raster(Path_Input)

Flwdir_D8= hydromt.flw.d8_from_dem(Dem ,max_depth=Max_Depth, outlets=Outlets, idxs_pit=None)

Flwdir_Ldd_Pyflwdir = xarray.apply_ufunc(pyflwdir.core_conversion.d8_to_ldd, Flwdir_D8.load())

Flwdir_Ldd_Pyflwdir.raster.to_raster(Path_Output)

return Flwdir_Ldd_Pyflwdir

import hydromt
import pyflwdir
import xarray

def DEM_2_LDD_HYDROMT(Path_Input, Path_Output, Max_Depth=-1.0, Outlets="edge"):

	# Read elevation data using hydromt which will transform the data into xarray
 
	Dem = hydromt.io.open_raster(Path_Input)
  
	Flwdir_D8= hydromt.flw.d8_from_dem(Dem ,max_depth=Max_Depth, outlets=Outlets, idxs_pit=None)
 
	Flwdir_Ldd_Pyflwdir = xarray.apply_ufunc(pyflwdir.core_conversion.d8_to_ldd, Flwdir_D8.load())
 
	Flwdir_Ldd_Pyflwdir.raster.to_raster(Path_Output)
   
	return Flwdir_Ldd_Pyflwdir

[shartgring]

That is a good question. The map with the flow directions uses 8 bit integer values. The nodata value has to also be an integer. In this case following LDD convention, it is 255.

I think that using a numpy NaN value np.nan could result in the map ending up as floats instead of 8 bit integers.

[AquaPore]

A great thanks, this explains very well why np.nan does not work .

[AquaPore]

Harmonizing SLOPE with LDD

The question is that we now created a new ldd , by buring the river and by filling the depressions, which has modified the DEM. Now I would need to use the modified DEM to compute the SLOPE.

The question, is there any ways of retrieving the corrected DEM?

This can be performed in pyflwdir.dem.fill_depressions but would not correct for the burning of the river?

Many thanks for any help you may provide,
Joseph

import geopandas
import os as os

Path_Input_River_Geojson = os.path.join(Input_Root, River_Input)

River_Geopandas = geopandas.read_file(Path_Input_River_Geojson)

River_GeoDataFrame = geopandas.GeoDataFrame(River_Geopandas)

import hydromt
from pyflwdir.core_conversion import d8_to_ldd
import xarray
import numpy as np


def DEM_2_LDD_HYDROMT(
    Path_Input,
    Path_Output,
    Ldd_MaxDepth=Ldd_MaxDepth,
    Ldd_RiverDepth=Ldd_RiverDepth,
    NoData=NoData,
    Outlets="edge",
    River=[],
):
    # Read elevation data using hydromt which will transform the data into xarray
    Dem = hydromt.io.open_raster(Path_Input)

    # Deriving the D8 from DEM
    Flwdir_D8 = hydromt.flw.d8_from_dem(
        da_elv=Dem,
        max_depth=Ldd_MaxDepth,
        gdf_riv=River,
        outlets=Outlets,
        idxs_pit=None,
        riv_burn_method="fixed",
        riv_depth=Ldd_RiverDepth,
    )

    # Converting D8-> Ldd as required by WFLOW model
    Flwdir_Ldd_Pyflwdir = xarray.apply_ufunc(d8_to_ldd, Flwdir_D8.load())

    # Save output
    Flwdir_Ldd_Pyflwdir.raster.to_raster(Path_Output, nodata=NoData)

    return Flwdir_Ldd_Pyflwdir

[AquaPore]

I think that I found a solution based on flwdir.dem_adjust
Flwdir