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* diff --git a/README.md b/README.md index 9ccb01f..2b8fc33 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 ``` @@ -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. diff --git a/jsc.2025.gnu.psmpi b/jsc.2025.gnu.psmpi index d1766f3..c6ea025 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,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 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..bed5555 100755 --- a/mkslopes/create_pfl_slopes +++ b/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 diff --git a/mkslopes/create_pfl_slopes.R b/mkslopes/create_pfl_slopes.R new file mode 100755 index 0000000..26f2796 --- /dev/null +++ b/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") 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') + 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 diff --git a/mkslopes/netCDF2simpleASCII.py b/mkslopes/netCDF2simpleASCII.py new file mode 100755 index 0000000..97a7e91 --- /dev/null +++ b/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") diff --git a/mkslopes/simpleASCII2netCDF.py b/mkslopes/simpleASCII2netCDF.py new file mode 100755 index 0000000..52956c8 --- /dev/null +++ b/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