Calculate slopes with PriorityFlow

.gitignore

@@ -0,0 +1,40 @@
+# ignore big data formats as
+*.tif
+*.tiff
+*.nc
+*.hdf
+*.sa
+*.pfb
+*.pfsol
+*.gdb
+*.zip
+*.tar
+
+# ignore shape files
+*.shp
+*.cpg
+*.dbf
+*.prj
+*.shx
+
+# ignore images
+*.png
+*.jpeg
+*.pdf
+
+# ignore swap files
+*.swp
+
+# ignore error / log files
+*.log
+*.LOG
+*err.*
+*out.*
+**/nohup.out
+
+# ignore compiled files
+*.pyc
+*.so
+
+# ignore MANIFEST
+*MANIFEST*

README.md

@@ -58,6 +58,7 @@ Above steps are performed by two scripts located in this directory:
 
 ```
 cd mklandmask/
+pip3 install geopandas
 python3 Breakthrough-Bosporus.py
 python3 make_land_lake_sea_mask.py
 ```
@@ -79,10 +80,12 @@ The correct river positions are mapped to the target grid (in this case hydroSHE
 Some pixels do need extra treatment, as for example the Elbe river in this setup.
 Those need manual adjustment and are corrected ‘pixel-by-pixel’.
 
-3) `sva_static_pfl.ncl`  
-GRASS algorithmus is taken to calculate flow direction and main river streams based on the burned DEM.
-To keep the correct slope values, those are calculated based on the original DEM, but the flow direction, represented by the sign of the slope value, is taken from flow direction calculated by GRASS algorithm.
+3) `create_pfl_slopes`
+The R package *PriorityFlow* is used to calculate flow direction and main river streams based on the burned DEM.
+There are several other tools available, a couple of them also doing D4 slopes.
+For slopes we use PriorityFlow, as [recommended by Reed Maxwell](https://github.com/parflow/parflow/issues/696#issuecomment-3920959452) who authored ParFlow and this tool.
 
+To keep the correct slope values, those are calculated based on the original DEM, but the flow direction, represented by the sign of the slope value, is taken from flow direction calculated by pysheds.
 This way the slope values are in line with the origin DEM, but flow direction is according to the correct river-corridors.
 
 #### Usage
@@ -98,15 +101,17 @@ unzip -u HydroRIVERS_v10_af_shp.zip
 ./modTopo.py
 ```
 
-We use [pysheds](https://mattbartos.com/pysheds/) to calculate slopes and flow directions.
-The following script will install and run pysheds.
-pysheds uses XArray as an internal format and reads and writes to disk in GeoTIFF format.
-We wrap around that, converting between netCDF and GeoTIFF.
+We use [pysheds](https://mattbartos.com/pysheds/) for pit filling, flooding depressions and resolving flats.
+We use [PriorityFlow](https://github.com/lecondon/PriorityFlow) to calculate the D4 slopes.
+The following script will install and run pysheds and PriorityFlow:
 
 ```
 ./create_pfl_slopes
 ```
 
+pysheds uses XArray as an internal format and reads and writes to disk in GeoTIFF format.
+We wrap around that by converting between netCDF and GeoTIFF.
+
 ## Creation of the mask and solids files
 
 ParFlow requires a "solid file" to distinguish between active and inactive cells in the simulation.

jsc.2025.gnu.psmpi

@@ -3,7 +3,7 @@ module --force purge
 export PYTHONPATH=""
 #
 module use $OTHERSTAGES
-module load Stages/2025
+module load Stages/2026
 module load GCC
 module load ParaStationMPI
 #
@@ -20,13 +20,16 @@ module load cURL
 #module load Szip
 module load Python
 #module load NCO
-module load ncview
+#module load ncview
 module load CMake
 #
 module load netcdf4-python
 module load SciPy-Stack
 module load Cartopy
 module load h5py
+module load R
+module load R-bundle-CRAN/2025.11
+module load GDAL
 #
 # cd $HOME/.local/share/
 # git clone --recurse-submodules https://github.com/HPSCTerrSys/SLOTH.git

mkslopes/create_pfl_slopes

@@ -2,40 +2,25 @@
 # Create D4 ParFlow slope files
 set -e
 
-fil_topo="./HSURFBurnedAndMod.nc"
-gisbase="/p/project1/cdetect/easybuild/jurecadc/software/GRASS/8.4.1-foss-2024a/grass8"
-gisdbase="$(mktemp -d -t pfl-XXXXXX)"
-
-# GRASS GIS extraction setup
-topo_nc=$fil_topo
-location_name="d4_map"
-#grassrc6=(/"GISDBASE: " + gisdbase, "LOCATION_NAME: " + location_name , "MAPSET: PERMANENT"/)
-#...
-
-export PATH="$gisbase/bin:$gisbase/scripts:$PATH"
-export LD_LIBRARY_PATH="$gisbase/lib:$LD_LIBRARY_PATH"
-export PYTHONPATH="$gisbase/etc/python"
-export GRASS_VERSION="8.4.1"
-export GIS_LOCK="$$"
-export tmp="$(mktemp -d -t grass6-XXXXXX)"
-export GISRC="$tmp/gisrc"
+topo_nc="./HSURFBurnedAndMod.nc"
+mask_nc="../mklandmask/EUR-11_TSMP_FZJ-IBG3_444x432_LAND-LAKE-SEA-MASK.nc"
+topo_tiff="HSURF-Sphere.tiff"
 
 # Convert netCDF to GeoTIFF
-# FIXME: HSURFBurnedAndMod.nc does not have geotransform, gcps, or rpcs
-# information included.  This was not preserved in burnShape2Topo.py.
-gdal_translate -a_srs ./WKT-Sphere.txt HSURF.nc HSURF-Sphere.tiff
-gdal_translate -a_srs ./WKT-Sphere.txt HSURFBurnedAndMod.nc HSURFBurnedAndMod-Sphere.tiff
+gdal_translate -a_srs ./WKT-Sphere.txt $topo_nc $topo_tiff
 
-# Use pysheds to calculate the location and direction of the streams, akin to
-# GRASS' r.watershed (sva_static_pfl.ncl)
-module load Python matplotlib
+# Use pysheds to fill pitts
+#module load Python matplotlib
 pip3 install pysheds
-python3 create_pfl_flow_direction.py
+python3 fill_pits.py $topo_tiff
 
-# Convert flow direction from GeoTIFF to netCDF
-gdal_translate -a_srs ./WKT-Sphere.txt flow_direction.tiff flow_direction.nc
-gdal_translate -a_srs ./WKT-Sphere.txt slopes.tiff slopes.nc
-#gdal_translate -a_srs ./WKT-Sphere.txt accumulation.tiff accumulation.nc
-#gdal_translate -a_srs ./WKT-Sphere.txt branches.tiff branches.nc
+# Convert from GeoTIFF to netCDF
+#gdal_translate -a_srs ./WKT-Sphere.txt flow_direction.tiff flow_direction.nc
+gdal_translate -a_srs ./WKT-Sphere.txt pit-filled_dem.tiff pit-filled_dem.nc
+
+# Use PriorityFlow to calculate the slopes, akin to
+# GRASS' r.watershed (sva_static_pfl.ncl)
+Rscript create_pfl_slopes.R "pit-filled_dem.nc" ${mask_nc}
 
-# TODO: better names for output field variables
+# Convert to simple ASCII for ParFlow
+./netCDF2simpleASCII.py slopes slopex

mkslopes/create_pfl_slopes.R

@@ -0,0 +1,45 @@
+#!/usr/bin/env R
+# Create D4 ParFlow slope files
+
+if (!interactive()) {
+    args <- commandArgs(trailingOnly = TRUE)
+    topo_nc <- args[1]
+    mask_nc <- args[2]
+}
+
+# Install PriorityFlow
+# On JSC HPC 'remotes' is already included in the R module
+if (!requireNamespace("remotes", quietly=TRUE)) install.packages('remotes')
+library(remotes)
+install_github("lecondon/PriorityFlow", subdir="Rpkg", quiet=TRUE)
+
+# Install packages that are dependencies but not described as such by PriorityFlow
+if (!requireNamespace("fields", quietly=TRUE)) install.packages('fields', repos="https://ftp.fau.de/cran/")
+
+# Use the PriorityFlow package to calculate slopes
+library('PriorityFlow')
+
+# netCDF support; part of R-bundle-CRAN/2025.11
+library('terra')
+
+llsm <- rast(mask_nc)
+llsm[llsm==2] <- 1  # treat lakes as land
+lsm_df <- data.frame(llsm)
+lsm_matrix = data.matrix(lsm_df)
+lsm_shaped = matrix(lsm_matrix, nrow=444)
+lsm_shaped = lsm_shaped[,432:1]  # up and down must be flipped
+hsurf <- rast(topo_nc)
+hsurf_df <- data.frame(hsurf)
+hsurf_matrix <- data.matrix(hsurf_df)
+hsurf_shaped <- matrix(hsurf_matrix, nrow=444)
+zero_matrix <- array( 0, dim(hsurf_shaped) )
+
+# Calculate slopes; ParFlow needs SlopeCalcUP()
+slope = PriorityFlow::SlopeCalStan( dem=hsurf_shaped, direction=zero_matrix,
+                                    dx=12500, dy=12500, mask=lsm_shaped )
+
+slopex_rast <- rast( t( slope$slopex ) )
+slopey_rast <- rast( t( slope$slopey ) )
+slope_dataset <- sds(slopex_rast, slopey_rast)
+varnames(slope_dataset) <- c("slopex", "slopey")
+writeCDF(slope_dataset, filename="slopes.nc")

mkslopes/fill_pits.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+import sys
+import numpy as np
+#from netCDF4 import Dataset
+from pysheds.grid import Grid
+
+try:
+    topo_tiff = sys.argv[1]
+except IndexError:
+    sys.exit(f"File name needed: '{sys.argv[0]} FILE'")
+
+# Read elevation raster
+grid = Grid.from_raster(topo_tiff)
+dem = grid.read_raster(topo_tiff)
+
+## Fill pits and depressions and resolve flats in DEM
+pit_filled_dem = grid.fill_pits(dem)
+flooded_dem = grid.fill_depressions(pit_filled_dem)
+inflated_dem = grid.resolve_flats(flooded_dem)
+
+# Write pit-filled to GeoTIFF
+grid.to_raster(data=inflated_dem, file_name='pit-filled_dem.tiff')
+

mkslopes/modTopo.py

@@ -20,9 +20,9 @@
 import sys
 import os
 import csv
-import fiona
 from shapely.geometry import shape
 import sloth.mapper
+import datetime
 
 ###############################################################################
 #### START OF: stuff you need to adjust!
@@ -32,12 +32,20 @@
     HSURFin = nc_file.variables['HSURFBurned'][...]
     lon2D   = nc_file.variables['lon'][...]
     lat2D   = nc_file.variables['lat'][...]
-# file to save adjusted topo to
-nc_file      = nc.Dataset(f'./HSURFBurnedAndMod.nc', 'w')
-latDim       = nc_file.createDimension('lat',HSURFin.shape[0])
-lonDim       = nc_file.createDimension('lon',HSURFin.shape[1])
-out_HSURF    = nc_file.createVariable(f'HSURFBurnedAndMod','f8',('lat','lon'),
-                                           zlib=True)
+
+rawdata_dir = os.environ['RAWDATA_DIR']
+griddesFileName = f'{rawdata_dir}/EUR-11_TSMP_FZJ-IBG3_CLMPFLDomain_444x432_griddes.txt'
+
+netCDFFileName_HSURFBurnedAndMod = sloth.IO.createNetCDF('./HSURFBurnedAndMod.nc',
+        domain=griddesFileName,
+        calcLatLon=True,
+        author='Marco VAN HULTEN',
+        contact='Marco.van.Hulten@uni-bonn.de',
+        institution='Uni Bonn',
+        history=f'Created: {datetime.datetime.now().strftime("%Y-%m-%d %H:%M")}',
+        description='HSURF with major rivers burnt and certain pixels modified',
+        source='')
+
 # pixel to adjust
 # Read in vertices of regions to adjust
 with open("./BurnRiversAndCanyons_Vertices_EUR11.csv", "r") as f:
@@ -55,9 +63,6 @@
         ToAdjust[regionID]['lat'].append(float(line[3]))
 
 # Initialize mapper from SLOTH
-#Mapper = sloth.mapper.mapper(SimLons=lon2D, SimLats=lat2D,
-#                            ObsLons=pointLons, ObsLats=pointLats,
-#                            ObsIDs=pointIDs)
 Mapper = sloth.mapper.mapper(SimLons=lon2D, SimLats=lat2D)
 for regionID in ToAdjust.keys():
     Mapper.ObsLons = np.array(ToAdjust[regionID]['lon'])
@@ -91,6 +96,12 @@
         print(f'handling: {Plat} | {Plon}')
         tmp_out_HSURF[Plat, Plon] = tmp_out_HSURF[LATref,LONref]
 
-out_HSURF[...] = tmp_out_HSURF[...]
-nc_file.close()
-
+# write results to netCDF object
+with nc.Dataset(netCDFFileName_HSURFBurnedAndMod, 'a') as nc_file:
+    nc_HSURFBurnedAndMod = nc_file.createVariable('HSURFBurnedAndMod', 'f8', ('rlat', 'rlon'), zlib=True)
+    nc_HSURFBurnedAndMod.standard_name = "HSURFBurnedAndMod"
+    nc_HSURFBurnedAndMod.long_name = "HSURFBurnedAndMod"
+    nc_HSURFBurnedAndMod.units = "m"
+    nc_HSURFBurnedAndMod.coordinates = "lon lat"
+    nc_HSURFBurnedAndMod.grid_mapping = "rotated_pole"
+    nc_HSURFBurnedAndMod[...] = tmp_out_HSURF

mkslopes/netCDF2simpleASCII.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+from netCDF4 import Dataset
+import numpy as np
+import sys
+
+try:
+    fname = sys.argv[1]
+except IndexError:
+    sys.exit(f"File and variable names needed: '{sys.argv[0]} FILE VARIABLE'")
+
+try:
+    vname = sys.argv[2]
+except IndexError:
+    sys.exit(f"File and variable names needed: '{sys.argv[0]} FILE VARIABLE'")
+
+nc = Dataset(f"{fname}.nc", 'r')
+var = nc.variables[f"{vname}"][:]
+nc.close()
+var_fortran = np.array(var, order='F')
+with open(f"{fname}.sa", 'w') as f:
+     f.write(f"{var.shape[1]} {var.shape[1]} {var.shape[0]}\n1\n")
+     for val in var_fortran.flatten():
+         f.write(f"{val}\n")

mkslopes/simpleASCII2netCDF.py

@@ -0,0 +1,30 @@
+#!/usr/bin/env python3
+from netCDF4 import Dataset
+import numpy as np
+import sys
+
+try:
+    sa_file = sys.argv[1]
+except IndexError:
+    sys.exit(f"File and variable names needed: '{sys.argv[0]} FILE VARIABLE'")
+
+try:
+    vname = sys.argv[2]
+except IndexError:
+    sys.exit(f"File and variable names needed: '{sys.argv[0]} FILE VARIABLE'")
+
+# Load ParFlow simple ASCII data
+with open(sa_file, 'r') as f:
+    nx, ny, nz = map(int, f.readline().split())
+var = np.loadtxt(sa_file, skiprows=1)
+myvar = var.reshape((nx, ny, nz), order='C')
+
+# Write to netCDF
+nc_file = f"{vname}.nc"
+with Dataset(nc_file, 'w', format='NETCDF4') as nc:
+    nc.createDimension('lon', nx)
+    nc.createDimension('lat', ny)
+    nc.createDimension('time', nz)
+    nc.description = 'Converted from ParFlow simple ASCII to netCDF'
+    var = nc.createVariable(f'{vname}', 'f4', ('time', 'lat', 'lon'))
+    var[:] = myvar