From 883bfaa67621de8ef540641f4b5f8e3a667ae287 Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Thu, 26 Feb 2026 11:06:35 +0100
Subject: [PATCH 01/15] Added nc to sa conversion script

---
 mkslopes/netCDF2simpleASCII.py | 12 ++++++++++++
 1 file changed, 12 insertions(+)
 create mode 100755 mkslopes/netCDF2simpleASCII.py

diff --git a/mkslopes/netCDF2simpleASCII.py b/mkslopes/netCDF2simpleASCII.py
new file mode 100755
index 0000000..750e5a0
--- /dev/null
+++ b/mkslopes/netCDF2simpleASCII.py
@@ -0,0 +1,12 @@
+#!/usr/bin/env python3
+from netCDF4 import Dataset
+import numpy as np
+
+nc = Dataset('HSURFBurnedAndMod.nc', 'r')
+var = nc.variables['HSURFBurnedAndMod'][:]
+nc.close()
+var_fortran = np.array(var, order='F')
+with open('HSURFBurnedAndMod.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")

From 61ff7efe5e4cca0567a22ce1bd6544256ce9362a Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Thu, 26 Feb 2026 11:41:58 +0100
Subject: [PATCH 02/15] Use PriorityFlow and moved to Stages/2026

The goal is to create slope files for ParFlow.  This tools should be able
to do that.
---
 jsc.2025.gnu.psmpi                    |  5 ++--
 mkslopes/create_pfl_flow_direction.py | 34 ------------------------
 mkslopes/create_pfl_slopes            | 37 +++------------------------
 mkslopes/create_pfl_slopes.R          | 11 ++++++++
 4 files changed, 17 insertions(+), 70 deletions(-)
 delete mode 100755 mkslopes/create_pfl_flow_direction.py
 create mode 100644 mkslopes/create_pfl_slopes.R

diff --git a/jsc.2025.gnu.psmpi b/jsc.2025.gnu.psmpi
index d1766f3..0d8e5e6 100644
--- a/jsc.2025.gnu.psmpi
+++ b/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,14 @@ 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
 #
 # cd $HOME/.local/share/
 # git clone --recurse-submodules https://github.com/HPSCTerrSys/SLOTH.git
diff --git a/mkslopes/create_pfl_flow_direction.py b/mkslopes/create_pfl_flow_direction.py
deleted file mode 100755
index 9d414cf..0000000
--- a/mkslopes/create_pfl_flow_direction.py
+++ /dev/null
@@ -1,34 +0,0 @@
-#!/usr/bin/env python3
-# Calculation of flow direction grid from HSURFBurnedAndMod-Sphere.tiff
-# based on https://mattbartos.com/pysheds/
-
-from pysheds.grid import Grid
-
-# Read elevation raster
-grid = Grid.from_raster('HSURF-Sphere.tiff')
-dem = grid.read_raster('HSURF-Sphere.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)
-
-# Determine D8 flow directions from DEM
-dirmap = (64, 128, 1, 2, 4, 8, 16, 32)
-#fdir = grid.flowdir(inflated_dem, dirmap=dirmap)
-fdir = grid.flowdir(dem, dirmap=dirmap)
-
-# Calculate the slopes
-slopes = grid.cell_slopes(dem, fdir)
-
-# Write flow direction and slope fields to GeoTIFF
-grid.to_raster(data=fdir, file_name='flow_direction.tiff')
-grid.to_raster(data=slopes, file_name='slopes.tiff')
-
-# Determine accumulation and write to GeoTIFF
-#acc = grid.accumulation(fdir, dirmap=dirmap)
-#grid.to_raster(data=acc, file_name='accumulation.tiff')
-
-# Determine river network
-# N.B. 'branches' is not a Raster but a FeatureCollection containing tuples!
-#branches = grid.extract_river_network(fdir, acc > 10, dirmap=dirmap)
diff --git a/mkslopes/create_pfl_slopes b/mkslopes/create_pfl_slopes
index 7848b28..2ed5021 100755
--- a/mkslopes/create_pfl_slopes
+++ b/mkslopes/create_pfl_slopes
@@ -2,40 +2,9 @@
 # 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)"
+topo_nc="./HSURFBurnedAndMod.nc"
 
-# 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"
-
-# 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
-
-# Use pysheds to calculate the location and direction of the streams, akin to
+# Use PriorityFlow to calculate the slopes, akin to
 # GRASS' r.watershed (sva_static_pfl.ncl)
-module load Python matplotlib
-pip3 install pysheds
-python3 create_pfl_flow_direction.py
-
-# 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
+Rscript create_pfl_slopes.R
 
-# TODO: better names for output field variables
diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R
new file mode 100644
index 0000000..f42f9a3
--- /dev/null
+++ b/mkslopes/create_pfl_slopes.R
@@ -0,0 +1,11 @@
+#!/usr/bin/env R
+# Create D4 ParFlow slope files
+
+# The 'remotes' package is to grab stuff from GitHub &co, but at JSC it is in R
+#install.packages('remotes')
+library(remotes)
+install_github("lecondon/PriorityFlow", subdir="Rpkg")
+
+install.packages('fields')
+
+library('PriorityFlow')

From 8415a817f1012c97e821e148b0fcaa7750ce97cd Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Thu, 26 Feb 2026 15:51:30 +0100
Subject: [PATCH 03/15] Call PriorityFlow::SlopeCalcUP()

This should be consistent with the ParFlow overland flow boundary condition.
---
 mkslopes/create_pfl_slopes   |  2 +-
 mkslopes/create_pfl_slopes.R | 17 +++++++++++++++--
 2 files changed, 16 insertions(+), 3 deletions(-)

diff --git a/mkslopes/create_pfl_slopes b/mkslopes/create_pfl_slopes
index 2ed5021..ffa58f6 100755
--- a/mkslopes/create_pfl_slopes
+++ b/mkslopes/create_pfl_slopes
@@ -6,5 +6,5 @@ topo_nc="./HSURFBurnedAndMod.nc"
 
 # Use PriorityFlow to calculate the slopes, akin to
 # GRASS' r.watershed (sva_static_pfl.ncl)
-Rscript create_pfl_slopes.R
+Rscript create_pfl_slopes.R ${topo_nc}
 
diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R
index f42f9a3..a797590 100644
--- a/mkslopes/create_pfl_slopes.R
+++ b/mkslopes/create_pfl_slopes.R
@@ -1,11 +1,24 @@
 #!/usr/bin/env R
 # Create D4 ParFlow slope files
 
-# The 'remotes' package is to grab stuff from GitHub &co, but at JSC it is in R
-#install.packages('remotes')
+topo_nc <- paste0(args[1])
+
+# 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")
 
+# Install packages that are dependencies but not described as such by PriorityFlow
 install.packages('fields')
 
+# Use the PriorityFlow package to calculate slopes
+# TODO: complete and correct arguments!
 library('PriorityFlow')
+PriorityFlow::SlopeCalcUP(
+            dem = topo_nc,
+            direction = , # e.g. flow_direction.{tiff,nc} from pysheds routine
+            dx = 0.11, dy = 0.11,
+            mask = "../mklandmask/EUR-11_TSMP_FZJ-IBG3_444x432_LAND-LAKE-SEA-MASK.nc",
+            borders = ,
+            rivermask = # e.g. RiverMask.nc created by burnShape2Topo.py
+)

From a1eb957bd7524f37ea0acfa454fa45ebd82dc209 Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Fri, 6 Mar 2026 17:27:32 +0100
Subject: [PATCH 04/15] Completed arguments for slopes calculation

This is calls PriorityFlow::SlopeCalStan() with the right arguments, but
ParFlow needs PriorityFlow::SlopeCalcUP().
---
 jsc.2025.gnu.psmpi           |  1 +
 mkslopes/create_pfl_slopes   |  5 +++--
 mkslopes/create_pfl_slopes.R | 40 ++++++++++++++++++++++++++----------
 3 files changed, 33 insertions(+), 13 deletions(-)

diff --git a/jsc.2025.gnu.psmpi b/jsc.2025.gnu.psmpi
index 0d8e5e6..f10123c 100644
--- a/jsc.2025.gnu.psmpi
+++ b/jsc.2025.gnu.psmpi
@@ -28,6 +28,7 @@ module load SciPy-Stack
 module load Cartopy
 module load h5py
 module load R
+module load R-bundle-CRAN/2025.11
 #
 # cd $HOME/.local/share/
 # git clone --recurse-submodules https://github.com/HPSCTerrSys/SLOTH.git
diff --git a/mkslopes/create_pfl_slopes b/mkslopes/create_pfl_slopes
index ffa58f6..3cd2787 100755
--- a/mkslopes/create_pfl_slopes
+++ b/mkslopes/create_pfl_slopes
@@ -2,9 +2,10 @@
 # Create D4 ParFlow slope files
 set -e
 
-topo_nc="./HSURFBurnedAndMod.nc"
+topo_nc="./HSURFBurned.nc" #AndMod
+mask_nc="/p/scratch/cdetect/vanhulten1/testing/TSMP_EUR-11.25A07/static.resource/03_LandLakeSeaMask/EUR-11_TSMP_FZJ-IBG3_444x432_LAND-LAKE-SEA-MASK.nc"
 
 # Use PriorityFlow to calculate the slopes, akin to
 # GRASS' r.watershed (sva_static_pfl.ncl)
-Rscript create_pfl_slopes.R ${topo_nc}
+Rscript create_pfl_slopes.R ${topo_nc} ${mask_nc}
 
diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R
index a797590..01860d2 100644
--- a/mkslopes/create_pfl_slopes.R
+++ b/mkslopes/create_pfl_slopes.R
@@ -1,24 +1,42 @@
 #!/usr/bin/env R
 # Create D4 ParFlow slope files
 
-topo_nc <- paste0(args[1])
+topo_nc <- args
+mask_nc <- args
 
+# 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")
 
 # Install packages that are dependencies but not described as such by PriorityFlow
-install.packages('fields')
+install.packages('fields', repos="https://ftp.fau.de/cran/")
 
 # Use the PriorityFlow package to calculate slopes
-# TODO: complete and correct arguments!
 library('PriorityFlow')
-PriorityFlow::SlopeCalcUP(
-            dem = topo_nc,
-            direction = , # e.g. flow_direction.{tiff,nc} from pysheds routine
-            dx = 0.11, dy = 0.11,
-            mask = "../mklandmask/EUR-11_TSMP_FZJ-IBG3_444x432_LAND-LAKE-SEA-MASK.nc",
-            borders = ,
-            rivermask = # e.g. RiverMask.nc created by burnShape2Topo.py
-)
+
+# 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)
+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 )
+#image(slope$slopex)
+#dev.off()
+
+slopex_rast <- rast( slope$slopex )
+slopey_rast <- rast( slope$slopey )
+slope_dataset <- sds(slopex_rast, slopey_rast)
+writeCDF(slope_rast, filename="slopes-out.nc", varname=slope_dataset)

From 6cb828954144efcd3229eded13c49ea4aedd7ea6 Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Wed, 11 Mar 2026 11:07:07 +0100
Subject: [PATCH 05/15] Write slope_dataset to netCDF

Slopes in x- and y-direction appear okay, except that the right half of both
variables are sort of a ghost mirror image of the left half.  Figure out at
which step this happens.
---
 README.md                    | 1 +
 mkslopes/create_pfl_slopes.R | 5 +++--
 2 files changed, 4 insertions(+), 2 deletions(-)
 mode change 100644 => 100755 mkslopes/create_pfl_slopes.R

diff --git a/README.md b/README.md
index 9ccb01f..8c2e423 100644
--- a/README.md
+++ b/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
 ```
diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R
old mode 100644
new mode 100755
index 01860d2..d1046bf
--- a/mkslopes/create_pfl_slopes.R
+++ b/mkslopes/create_pfl_slopes.R
@@ -11,7 +11,7 @@ library(remotes)
 install_github("lecondon/PriorityFlow", subdir="Rpkg")
 
 # Install packages that are dependencies but not described as such by PriorityFlow
-install.packages('fields', repos="https://ftp.fau.de/cran/")
+if (!requireNamespace("fields", quietly=TRUE)) install.packages('fields', repos="https://ftp.fau.de/cran/")
 
 # Use the PriorityFlow package to calculate slopes
 library('PriorityFlow')
@@ -39,4 +39,5 @@ slope = PriorityFlow::SlopeCalStan( dem=hsurf_shaped, direction=zero_matrix,
 slopex_rast <- rast( slope$slopex )
 slopey_rast <- rast( slope$slopey )
 slope_dataset <- sds(slopex_rast, slopey_rast)
-writeCDF(slope_rast, filename="slopes-out.nc", varname=slope_dataset)
+varnames(slope_dataset) <- c("slopex", "slopey")
+writeCDF(slope_dataset, filename="slopes-out.nc")

From b85fa33db107e61bac0ddd1161bd1f058984aecd Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Wed, 11 Mar 2026 16:47:00 +0100
Subject: [PATCH 06/15] Generalised conversion script and call from
 create_pfl_slopes

Handling of arguments in R and Python
---
 mkslopes/create_pfl_slopes     |  2 ++
 mkslopes/create_pfl_slopes.R   |  9 ++++++---
 mkslopes/netCDF2simpleASCII.py | 12 +++++++++---
 3 files changed, 17 insertions(+), 6 deletions(-)

diff --git a/mkslopes/create_pfl_slopes b/mkslopes/create_pfl_slopes
index 3cd2787..319d661 100755
--- a/mkslopes/create_pfl_slopes
+++ b/mkslopes/create_pfl_slopes
@@ -9,3 +9,5 @@ mask_nc="/p/scratch/cdetect/vanhulten1/testing/TSMP_EUR-11.25A07/static.resource
 # GRASS' r.watershed (sva_static_pfl.ncl)
 Rscript create_pfl_slopes.R ${topo_nc} ${mask_nc}
 
+# Convert to simple ASCII for ParFlow
+./netCDF2simpleASCII.py slopes-out
diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R
index d1046bf..57a5002 100755
--- a/mkslopes/create_pfl_slopes.R
+++ b/mkslopes/create_pfl_slopes.R
@@ -1,14 +1,17 @@
 #!/usr/bin/env R
 # Create D4 ParFlow slope files
 
-topo_nc <- args
-mask_nc <- args
+if (!interactive()) {
+    args <- commandArgs(trailingOnly = TRUE)
+    topo_nc <- args
+    mask_nc <- args
+}
 
 # 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")
+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/")
diff --git a/mkslopes/netCDF2simpleASCII.py b/mkslopes/netCDF2simpleASCII.py
index 750e5a0..4da77f2 100755
--- a/mkslopes/netCDF2simpleASCII.py
+++ b/mkslopes/netCDF2simpleASCII.py
@@ -1,12 +1,18 @@
 #!/usr/bin/env python3
 from netCDF4 import Dataset
 import numpy as np
+import sys
 
-nc = Dataset('HSURFBurnedAndMod.nc', 'r')
-var = nc.variables['HSURFBurnedAndMod'][:]
+try:
+    name = sys.argv[1]
+except IndexError:
+    sys.exit(f"Argument needed: '{sys.argv[0]} myvar' with myvar a variable in myvar.nc.")
+
+nc = Dataset(f"{name}.nc", 'r')
+var = nc.variables[f"{name}"][:]
 nc.close()
 var_fortran = np.array(var, order='F')
-with open('HSURFBurnedAndMod.sa', 'w') as f:
+with open(f"{name}.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")

From 14aa2d20c84a19965bade481d69474d5c0bf2cc5 Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Thu, 26 Mar 2026 16:19:31 +0100
Subject: [PATCH 07/15] Modified to resolve issues with the grid by using SLOTH

I am not certain if this is relevant, as I've found out that the script
(also in its previous form) does nothing to HSURFBurned.nc, i.e. the
values of HSURFBurnedAndMod are identical to HSURFBurned.
---
 mkslopes/modTopo.py | 37 ++++++++++++++++++++++++-------------
 1 file changed, 24 insertions(+), 13 deletions(-)

diff --git a/mkslopes/modTopo.py b/mkslopes/modTopo.py
index 318065d..591e16c 100755
--- a/mkslopes/modTopo.py
+++ b/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

From 2e5edc13a9d39c474810a24953876faa082a2a4f Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Fri, 27 Mar 2026 11:23:33 +0100
Subject: [PATCH 08/15] Pass arguments correctly and transpose the slopes

N.B. This still uses PriorityFlow::SlopeCalStan() and should use
PriorityFlow::SlopeCalcUP() for ParFlow.  Also direction based on non-burnt DEM
should be used in the calculation.

Resolves: #10
---
 README.md                    |  4 ++--
 mkslopes/create_pfl_slopes.R | 12 +++++-------
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/README.md b/README.md
index 8c2e423..c26f99d 100644
--- a/README.md
+++ b/README.md
@@ -80,8 +80,8 @@ 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.
+3) `create_pfl_slopes`
+The R package *PriorityFlow* is used 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.
 
 This way the slope values are in line with the origin DEM, but flow direction is according to the correct river-corridors.
diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R
index 57a5002..667dd8d 100755
--- a/mkslopes/create_pfl_slopes.R
+++ b/mkslopes/create_pfl_slopes.R
@@ -3,8 +3,8 @@
 
 if (!interactive()) {
     args <- commandArgs(trailingOnly = TRUE)
-    topo_nc <- args
-    mask_nc <- args
+    topo_nc <- args[1]
+    mask_nc <- args[2]
 }
 
 # Install PriorityFlow
@@ -36,11 +36,9 @@ 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 )
-#image(slope$slopex)
-#dev.off()
 
-slopex_rast <- rast( slope$slopex )
-slopey_rast <- rast( slope$slopey )
+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-out.nc")
+writeCDF(slope_dataset, filename="slopes.nc")

From 7f3ac124626846715f7faefb284f3da55f1ca2fc Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Wed, 1 Apr 2026 11:48:31 +0200
Subject: [PATCH 09/15] Added an sa-to-netCDF conversion script

This should be the opposite of netCDF2simpleASCII.py
---
 mkslopes/simpleASCII2netCDF.py | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100755 mkslopes/simpleASCII2netCDF.py

diff --git a/mkslopes/simpleASCII2netCDF.py b/mkslopes/simpleASCII2netCDF.py
new file mode 100755
index 0000000..f3e1e59
--- /dev/null
+++ b/mkslopes/simpleASCII2netCDF.py
@@ -0,0 +1,32 @@
+#!/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'")
+
+print(sa_file)
+
+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

From 166a6d409a43a34ef9ba56a2307a54d9fa3cd291 Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Wed, 1 Apr 2026 14:56:30 +0200
Subject: [PATCH 10/15] Completed slopes workflow

N.B. This does not generate the expected slopes!

Does not resolve #10

+ Use pysheds to fill pits
+ Use GDAL for conversion between GeoTIFF and netCDF
---
 README.md                      | 14 +++++++++-----
 jsc.2025.gnu.psmpi             |  1 +
 mkslopes/create_pfl_slopes     | 19 ++++++++++++++++---
 mkslopes/netCDF2simpleASCII.py | 15 ++++++++++-----
 mkslopes/simpleASCII2netCDF.py |  2 --
 5 files changed, 36 insertions(+), 15 deletions(-)

diff --git a/README.md b/README.md
index c26f99d..2b8fc33 100644
--- a/README.md
+++ b/README.md
@@ -82,8 +82,10 @@ Those need manual adjustment and are corrected ‘pixel-by-pixel’.
 
 3) `create_pfl_slopes`
 The R package *PriorityFlow* is used 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.
+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
@@ -99,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.
diff --git a/jsc.2025.gnu.psmpi b/jsc.2025.gnu.psmpi
index f10123c..c6ea025 100644
--- a/jsc.2025.gnu.psmpi
+++ b/jsc.2025.gnu.psmpi
@@ -29,6 +29,7 @@ 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
diff --git a/mkslopes/create_pfl_slopes b/mkslopes/create_pfl_slopes
index 319d661..3a2466e 100755
--- a/mkslopes/create_pfl_slopes
+++ b/mkslopes/create_pfl_slopes
@@ -2,12 +2,25 @@
 # Create D4 ParFlow slope files
 set -e
 
-topo_nc="./HSURFBurned.nc" #AndMod
+topo_nc="./HSURFBurnedAndMod.nc"
 mask_nc="/p/scratch/cdetect/vanhulten1/testing/TSMP_EUR-11.25A07/static.resource/03_LandLakeSeaMask/EUR-11_TSMP_FZJ-IBG3_444x432_LAND-LAKE-SEA-MASK.nc"
 
+# Convert netCDF to GeoTIFF
+gdal_translate -a_srs ./WKT-Sphere.txt HSURF.nc HSURF-Sphere.tiff
+gdal_translate -a_srs ./WKT-Sphere.txt HSURFBurnedAndMod.nc HSURFBurnedAndMod-Sphere.tiff
+
+# Use pysheds to fill pitts
+#module load Python matplotlib
+pip3 install pysheds
+python3 fill_pits.py
+
+# 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 ${topo_nc} ${mask_nc}
+Rscript create_pfl_slopes.R "pit-filled_dem.nc" ${mask_nc}
 
 # Convert to simple ASCII for ParFlow
-./netCDF2simpleASCII.py slopes-out
+./netCDF2simpleASCII.py slopes slopex
diff --git a/mkslopes/netCDF2simpleASCII.py b/mkslopes/netCDF2simpleASCII.py
index 4da77f2..97a7e91 100755
--- a/mkslopes/netCDF2simpleASCII.py
+++ b/mkslopes/netCDF2simpleASCII.py
@@ -4,15 +4,20 @@
 import sys
 
 try:
-    name = sys.argv[1]
+    fname = sys.argv[1]
 except IndexError:
-    sys.exit(f"Argument needed: '{sys.argv[0]} myvar' with myvar a variable in myvar.nc.")
+    sys.exit(f"File and variable names needed: '{sys.argv[0]} FILE VARIABLE'")
 
-nc = Dataset(f"{name}.nc", 'r')
-var = nc.variables[f"{name}"][:]
+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"{name}.sa", 'w') as 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")
diff --git a/mkslopes/simpleASCII2netCDF.py b/mkslopes/simpleASCII2netCDF.py
index f3e1e59..52956c8 100755
--- a/mkslopes/simpleASCII2netCDF.py
+++ b/mkslopes/simpleASCII2netCDF.py
@@ -8,8 +8,6 @@
 except IndexError:
     sys.exit(f"File and variable names needed: '{sys.argv[0]} FILE VARIABLE'")
 
-print(sa_file)
-
 try:
     vname = sys.argv[2]
 except IndexError:

From 9d7f0c5c7e415634306338b662321544c9828e0f Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Thu, 2 Apr 2026 17:12:10 +0200
Subject: [PATCH 11/15] Use the generated land-sea-mask

---
 mkslopes/create_pfl_slopes | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/mkslopes/create_pfl_slopes b/mkslopes/create_pfl_slopes
index 3a2466e..5a5c5b1 100755
--- a/mkslopes/create_pfl_slopes
+++ b/mkslopes/create_pfl_slopes
@@ -3,7 +3,7 @@
 set -e
 
 topo_nc="./HSURFBurnedAndMod.nc"
-mask_nc="/p/scratch/cdetect/vanhulten1/testing/TSMP_EUR-11.25A07/static.resource/03_LandLakeSeaMask/EUR-11_TSMP_FZJ-IBG3_444x432_LAND-LAKE-SEA-MASK.nc"
+mask_nc="../mklandmask/EUR-11_TSMP_FZJ-IBG3_444x432_LAND-LAKE-SEA-MASK.nc"
 
 # Convert netCDF to GeoTIFF
 gdal_translate -a_srs ./WKT-Sphere.txt HSURF.nc HSURF-Sphere.tiff

From 96ea2d6034c867618b49062400fd6b6044f6cdda Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Fri, 10 Apr 2026 11:45:59 +0200
Subject: [PATCH 12/15] Reversed indexing in latitude of mask

Apparently the mask is upside down; with this patch the mask is applied
correctly.

This commit resolves Problem 2 of issue #10.
---
 mkslopes/create_pfl_slopes.R | 1 +
 1 file changed, 1 insertion(+)

diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R
index 667dd8d..26f2796 100755
--- a/mkslopes/create_pfl_slopes.R
+++ b/mkslopes/create_pfl_slopes.R
@@ -27,6 +27,7 @@ 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)

From 15cacc7f7bda20ce05d2e8e7d73003a83d6b605c Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Thu, 16 Apr 2026 17:14:05 +0200
Subject: [PATCH 13/15] Pass variables as filenames

---
 mkslopes/create_pfl_slopes | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/mkslopes/create_pfl_slopes b/mkslopes/create_pfl_slopes
index 5a5c5b1..bed5555 100755
--- a/mkslopes/create_pfl_slopes
+++ b/mkslopes/create_pfl_slopes
@@ -4,15 +4,15 @@ set -e
 
 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
-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 fill pitts
 #module load Python matplotlib
 pip3 install pysheds
-python3 fill_pits.py
+python3 fill_pits.py $topo_tiff
 
 # Convert from GeoTIFF to netCDF
 #gdal_translate -a_srs ./WKT-Sphere.txt flow_direction.tiff flow_direction.nc

From 006687d1fea4f6fdc31798596ab34631e67810fc Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Fri, 17 Apr 2026 09:54:41 +0200
Subject: [PATCH 14/15] Ignore files that should not be tracked

---
 .gitignore | 40 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 40 insertions(+)
 create mode 100644 .gitignore

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..9d73acd
--- /dev/null
+++ b/.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*

From 21eac3a91ee8e49822ac32f4d39490d6c16ea27b Mon Sep 17 00:00:00 2001
From: Marco van Hulten <Marco.van.Hulten@uni-bonn.de>
Date: Fri, 17 Apr 2026 11:39:33 +0200
Subject: [PATCH 15/15] Fill pits and depressions and resolve flats

This uses pysheds.
---
 mkslopes/fill_pits.py | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)
 create mode 100755 mkslopes/fill_pits.py

diff --git a/mkslopes/fill_pits.py b/mkslopes/fill_pits.py
new file mode 100755
index 0000000..39b367d
--- /dev/null
+++ b/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')
+