slope generation

.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

@@ -6,6 +6,13 @@ The generator incorporates steps from [the TSMP1 static file generator](https://
 Instead of including raw input data in this repository, the repo is instead kept small.
 The static-files input data can be found [here](https://icg4geo.icg.kfa-juelich.de/ExternalReposPublic/tsmp2-static-files/grids_parflow_cordex-eur-11u).
 
+We depend on the [SLOTH helper scripts](https://hpscterrsys.github.io/SLOTH/README.html), which can be installed through
+
+```
+cd $HOME/.local/share/
+git clone --recurse-submodules https://github.com/HPSCTerrSys/SLOTH.git
+```
+
 If you are running this generator on a [JSC](https://www.fz-juelich.de/en/ias/jsc) machine, sourcing the provided environment file
 
 ```
@@ -14,6 +21,7 @@ source jsc.2025.gnu.psmpi
 
 makes the necessary utilities and libraries available.
 Otherwise you have to make sure that the respective software is installed or made available on your system.
+On OpenBSD you need the packages `curl cmake eccodes gcc hdf5 ncview netcdf openmpi py3-{gdal,geopandas,h5py,matplotlib,netcdf4,scipy,xarray} tcl-8.6`.
 
 To create all static files needed to run ParFlow, you first need to download static input files:
 
@@ -40,24 +48,110 @@ However, this classification process introduces some artifacts along the coastli
 To address these artifacts, we run an optimization loop that checks if lake pixels are neighbors of sea pixels.
 If they are, we treat them as sea pixels as well, ensuring a more accurate representation of the coastline.
 
-It is important to note that COSMO also uses the `FR_LAND` variable as a land-sea mask with also a threshold of 0.5.
+It is important to note that COSMO and ICON also use the `FR_LAND` variable as a land-sea mask with also a threshold of 0.5.
 This ensures that all components of TSMP share the same land-sea mask.
 
 Furthermore, special treatment is required due to the coarse horizontal resolution of our target grid.
 Some topographic formations, like the Bosporus connecting the Black Sea and the Mediterranean Sea, cannot be adequately resolved.
-To address this, we modify the original `FR_LAND` variable by setting the value to 0, based on a predefined shape-files containing the coordinates of the Bosporus break through, effectively designating those areas as total water pixels.
+To address this, we modify the original `FR_LAND` variable by setting the value to 0, based on a predefined shape-files containing the coordinates of the Bosporus breakthrough, effectively designating those areas as total water pixels.
 
 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
 ```
 
 ## Creation of the flow direction and slopes
 
-## Creation of the mask solids mask
+Using DEMs straight forward to calculate slopes for ParFlow could lead to smaller or bigger problems in river corridor placement.
+In particular for coarse spatial resolutions this is easy to imagine as tight canyons are smoothed out.
+For example a 12km resolution as for the EU11 domain does not see the breakthrough valley [`irongate` for donau river](https://de.wikipedia.org/wiki/Eisernes_Tor) leading to the result, that the donau is flowing around the related mountain range.
+A further example are the Netherlands, where major parts of the land area are below sea level and rivers do not follow the ‘natural’ river-corridor but are forced to follow artificial canals.
+To fix those and other issues / problems the real river-corridors are ‘burned’ to the DEM within this approach.
+
+Burning the correct river corridors is achieved in two steps:
+
+1) `burnShape2Topo.py`  
+The correct river positions are mapped to the target grid (in this case hydroSHEDS data were used) and then river pixels (and neighboring ones) are pushed down within the DEM.
+
+2) `modTopo.py`  
+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’.
+
+The Python package *priority\_flow* 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.
+ParFlow author Reed Maxwell [recommended](https://github.com/parflow/parflow/issues/696#issuecomment-3920959452) to use [PriorityFlow](https://github.com/lecondon/PriorityFlow), which was since ported to Python as [priority\_flow](https://github.com/hydroframe/priority_flow), which we will use here.
+
+#### Usage
+First extract and adjust the HydroSHEDS river data:
+
+```
+cd ../mkslopes
+wget https://data.hydrosheds.org/file/HydroRIVERS/HydroRIVERS_v10_eu_shp.zip
+wget https://data.hydrosheds.org/file/HydroRIVERS/HydroRIVERS_v10_af_shp.zip
+unzip -u HydroRIVERS_v10_eu_shp.zip
+unzip -u HydroRIVERS_v10_af_shp.zip
+./burnShape2Topo.py
+./modTopo.py
+```
+
+#### Optional / to be tested / unsure whether pysheds is needed
+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.
+
+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.
+
+#### Usage of priority\_flow for calculating slopes
+
+```
+mkdir data
+pip3 install priority_flow
+./create_pfl_slopes.py
+```
+
+**N.B. Even though this can create D4 slopes, they are untested and may not be useful with ParFlow.
+Further development and testing of the slopes generation is needed.
+If you want to run TSPM2 with ParFlow, you may create the land mask, solids and textures with this generator, but for slopes it is probably best to rely on 
+[the existing slopes](https://gitlab.jsc.fz-juelich.de/detect/detect_z03_z04/constant_fields/TSMP_EUR-11/-/tree/main/static/parflow).**
+
+## Creation of the mask and solids files
+
+ParFlow requires a "solid file" to distinguish between active and inactive cells in the simulation.
+In this case, active cells represent land areas, while inactive cells represent ocean areas.
+To generate this solid file, we use the land-sea mask that was previously created.
+
+However, it is important to note that we want to preserve lakes in the simulation, as rivers should not disappear when they flow into a lake.
+Therefore, the solid file is created based solely on the presence or absence of ocean areas, without masking out the lakes.
+
+The purpose of the solid file is to inform ParFlow that ocean areas should be treated as inactive cells, meaning they are not considered for computation.
+On the other hand, land areas (including lakes) should be treated as active cells, allowing for proper simulation of water flow and interactions.
+
+#### Usage
+
+Set `PARFLOW_DIR` to your ParFlow installation prefix (the directory that contains ParFlow's `bin`, `lib` and such).
+In the mask and solid files creation below, we exemplified this by assuming you are using the TSMP2 framework on JURECA:
+
+```
+cd ../mksolids
+./createPfbMask.py
+export PARFLOW_DIR="$TSMP2_DIR/bin/JURECADC_ParFlow" # modify to install prefix
+$PARFLOW_DIR/bin/pfmask-to-pfsol --z-top 30 --z-bottom 0 --mask PfbMask4SolidFile.pfb --pfsol PfbMask4SolidFile.pfsol >PfbMask4SolidFile.log
+```
+
+You have to choose the value for `--z-top` according to `ComputationalGrid.DZ * ComputationalGrid.NZ`.
+In the case of EUR-11 this comes down to $$z_\textrm{bottom} = 2 \cdot 15$$.
 
 ## Creation of the texture indicator