Skip to content

Commit

Permalink
changes to files
Browse files Browse the repository at this point in the history
  • Loading branch information
@Ash-Manoj
Ash-Manoj committed Aug 6, 2024
1 parent 4eb67a3 commit 63c8c51
Show file tree
Hide file tree
Showing 4 changed files with 397 additions and 270 deletions.
9 changes: 3 additions & 6 deletions Dockerfile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,12 +4,9 @@ FROM r-base:4.2.0
# the parameter parsing function always needs the rjson and yaml packages
RUN R -e "install.packages(c('jsonlite', 'yaml'))"

# install json2aRgs v0.3.0 to parse parameters from /in/parameters.json
RUN wget -q https://cran.r-project.org/src/contrib/json2aRgs_0.3.0.tar.gz
RUN R CMD INSTALL json2aRgs_0.3.0.tar.gz

# remove json2aRgs tarball
RUN rm json2aRgs_0.3.0.tar.gz
# install json2aRgs latest version from github
RUN R -e 'install.packages("remotes")'
RUN R -e 'remotes::install_github("VForWaTer/json2aRgs")'

# install Catflow-R-Package dependencies
RUN R -e "install.packages(c('deSolve', 'RColorBrewer', 'zoo', 'xts'))"
Expand Down
279 changes: 15 additions & 264 deletions src/catlib.R
View file
Original file line number Diff line number Diff line change
@@ -1,275 +1,26 @@
#--------------------------------------------------------------------------------------
# This file essentially maps the Catflow preprocessing functions from the
# Catflow-R-Pacakge and makes them available in the container.
# "Real" workflows which combine the different functions can be found in lib.R
# This file additionally contains a function based on the hillslope wizards developed
# This file contains a function based on the hillslope wizards developed
# by Ralf Loritz (2014). The function is used to infer a representative hillslope from
# provided .tif files.

# load Catflow package
library(Catflow)

catlib_make_geometry <- function(params) {
# output file is always saved to /out/CATFLOW/in/
params$out.file <- "geometry.geo"

# create project folder
project.path <- "/out/CATFLOW/in"
system("mkdir -p /out/CATFLOW/in/")
system("chmod 777 /out/CATFLOW")
system("chmod 777 /out/CATFLOW/in")

# create folder for plots
system("mkdir -p /out/plots")
system("chmod 777 /out/plots")

# run function make.geometry() with params as input parameters
pdf("/out/plots/geometry.pdf")
output <- make.geometry(params, plottitle = "Model Geometry", project.path = project.path)
dev.off()

# save output of make.geometry() for use in other tools
saveRDS(output, file = "/out/geom.Rds")
}

catlib_write_facmat <- function(params) {
# load geometry generated by the tool make_geometry and attach
geom <- readRDS(params$geometry)
attach(geom)

# create project folder
system("mkdir -p /out/CATFLOW/in")
system("chmod 777 /out/CATFLOW")
system("chmod 777 /out/CATFLOW/in")

# multipliers for scaling saturated hydraulic conductivity
write.facmat(output.file = "/out/CATFLOW/in/ksmult.dat")

# multipliers for scaling saturated water content / porosity
write.facmat(output.file = "/out/CATFLOW/in/thsmult.dat")

# write initial conditions for hydraulic conductivity
if (params$write_soilhyd_ini) {
write.facmat(output.file = "/out/CATFLOW/in/soilhyd.ini",
headr = paste("PSI ", 0, 1, length(eta), length(xsi), 1))
}

# write soiltype identifiers
if (params$write_soil_types) {
write.facmat(output.file = "/out/CATFLOW/in/soils.bod",
headr = paste("BODEN", length(eta), length(xsi), 1),
fac = matrix(c(rep(1, ceiling(length(eta) / 2)),
rep(2, floor(length(eta) / 2))),
nrow = length(eta), ncol = length(xsi)))
}
}

catlib_write_precip <- function(params) {
# create project folder
system("mkdir -p /out/CATFLOW/in")
system("chmod 777 /out/CATFLOW")
system("chmod 777 /out/CATFLOW/in")

# write precipitation data with params as input
write.precip(
raindat = params$raindat,
output.file = "/out/CATFLOW/in/rain.dat",
start.time = params$start.time,
time.unit = params$time.unit,
faktor.p = params$faktor.p
)

# create folder for plots
system("mkdir -p /out/plots")
system("chmod 777 /out/plots")

# plot rainfall data
pdf("/out/plots/raindat.pdf")
plot(params$raindat, t = "s", xlab = "time", ylab = "precipitation")
dev.off()
}

catlib_write_climate <- function(params) {
# create project folder
system("mkdir -p /out/CATFLOW/in")
system("chmod 777 /out/CATFLOW")
system("chmod 777 /out/CATFLOW/in")

# write climate data with params as input
write.climate(
climadat = params$climadat,
output.file = "/out/CATFLOW/in/clima.dat",
start.time = params$start.time,
time.unit = params$time.unit,
rBilart = params$rBilart,
ref.height = params$ref.height,
sw0 = params$sw0,
sw1 = params$sw1,
sw2 = params$sw2,
trueb = params$trueb,
truebf = params$truebf,
NA.flag = params$NA.flag
)
}

catlib_write_printout <- function(params) {
# create project folder
system("mkdir -p /out/CATFLOW/in")
system("chmod 777 /out/CATFLOW")
system("chmod 777 /out/CATFLOW/in")

# write printout times
write.printout(
output.file = "/out/CATFLOW/in/printout.prt",
start.time = params$start.time,
end.time = params$end.time,
intervall = params$intervall,
time.unit = params$time.unit,
flag = params$flag
)
}

catlib_write_surface_pob <- function(params) {
# load geometry generated by the tool make_geometry and attach
geom <- readRDS(params$geometry)
attach(geom)

# create project folder
system("mkdir -p /out/CATFLOW/in")
system("chmod 777 /out/CATFLOW")
system("chmod 777 /out/CATFLOW/in")

# drop geometry from params (unused in write.surface.pob)
params$geometry <- NULL

# write surface attributes
write.surface.pob(
output.file = "/out/CATFLOW/in/surface.pob",
xs = xsi,
lu = params$lu,
precid = params$precid,
climid = params$climid,
windid = params$windid
)
}

catlib_write_control <- function(params) {
# create the project folder and set permissions
system("mkdir -p /out/CATFLOW")
system("chmod 777 /out/CATFLOW")

# write project controll file
write.control(
output.file = "run_cat.in",
project.path = "/out/CATFLOW",
start.date = params$start.time,
end.date = params$end.time,
slope.in.list = list(
slope1 = list(
geo.file = "geometry.geo",
soil.file = "soils.bod",
ks.fac = "ksmult.dat",
ths.fac = "thsmult.dat",
macro.file = "profil.mak",
cv.file = "cont_vol.cv",
ini.file = "soilhyd.ini",
print.file = "printout.prt",
surf.file = "surface.pob",
bc.file = "boundary.rb"
)
)
)

# write main control file
write.CATFLOW.IN(
control.files = "run_cat.in",
project.path = "/out/CATFLOW"
)
}

catlib_complete_file_structure <- function(params) {
# create the project folder CATFLOW and the in/ folder and set permissions
system("mkdir -p /out/CATFLOW/in")
system("chmod 777 /out/CATFLOW")
system("chmod 777 /out/CATFLOW/in")

# macro.file: "profil.mak"
cat(paste("1 0 2", "ari", "0.00 1.00 0.00 1.00 1 1.00 1.00 ", sep = "\n"), file = "/out/CATFLOW/in/profil.mak")

# cv.file: "cont_vol.cv"
cat(paste("1", "0.8 0.9 0.98 1.0", sep = "\n"), file = "/out/CATFLOW/in/cont_vol.cv")

# bc.file: "boundary.rb"
cat(paste("L", "1 0", "0. 1. 0", " ",
"R", "1 0", "0. 1. -10", " ",
"T", "1 0", "0. 1. -99 ", " ",
"B", "1 0", "0. 1. 0", " ",
"S", "1 0", "0. 1. 0. 1. -99", " ",
"M", 0, sep = "\n"),
file = "/out/CATFLOW/in/boundary.rb")

# winddir.file: "winddir.def"
cat(paste("4", "240 0.81", " 50 0.78", " 80 0.97", "220 0.94", sep = "\n"), file = "/out/CATFLOW/in/winddir.def")

# soildef.file: soils.def
# writes a soil type definition for two soil types:
cat(paste("2", "1 Loamy Sand, porosity 0.55, bulk dens 1 g/cm3",
"1 800 1. 1. 1e-4 0.5 0.34 0.11 20. 0.70 0.050 1. 1. 1.",
"4.05e-5 0.55 0.06 12.40 2.28 -6.00 8.00 1000.00 0.80",
"0. 0. 0.", "0. 0. 0.", "0. 0. 0.",
"2 Sandy Clay Loam (30% S, 40 % U; 30 % T)",
"1 800 1. 1. 1e-4 0.5 0.34 0.11 20. 0.70 0.050 1. 1. 1.",
"3.42e-6 0.48 0.08 0.96 1.5 -6.00 8.00 1200.00 0.80",
"0. 0. 0.", "0. 0. 0.", "0. 0. 0.", sep = "\n"),
file = "/out/CATFLOW/in/soils.def")

# timeser.file: timeser.def
cat(paste("PREC", "1", "in/rain.dat", "",
"BC", "0", "", "SINKS", "0", "", "SOLUTE", "0", "",
"LAND-USE", "in/landuse/lu_ts.dat", "",
"CLIMATE", "1", "in/clima.dat", "", sep = "\n"),
file = "/out/CATFLOW/in/timeser.def")

# file related to landuse specification
# create landuse subdirectory and set permissions
system("mkdir -p /out/CATFLOW/in/landuse")
system("chmod 777 /out/CATFLOW/in/landuse")

# pointer to land-use parameters
cat(paste("3", "coniferous forest", "in/landuse/conif.par", sep = " "), file = "/out/CATFLOW/in/landuse/lu_file.def")

# time-series of land-use parameters
cat(paste("01.01.2004 00:00:00.00", "in/landuse/lu_set1.dat", "01.01.2005 00:00:00.00", sep = "\n"), file = "/out/CATFLOW/in/landuse/lu_ts.dat")

# land-use parameters
cat(paste(
paste("10", "KST", "MAK", "BFI", "BBG", "TWU", "PFH",
"PALB", "RSTMIN", "WP_BFW", "F_BFW", sep = " "),
"0. 3. 1. 5. 0.95 5.0 5.0 0.15 1. 1. 1.",
paste(c(" 1", "366"),
"2. 1. 1. 1.0 1.0 1.0 1.0 546. 0.05 30.",
sep = " ", collapse = "\n"), sep = "\n"),
file = "/out/CATFLOW/in/landuse/conif.par")

# pointer to surface node attributes
cat(paste(1, "33 3 %coniferous forest", sep = "\n"), file = "/out/CATFLOW/in/landuse/lu_set1.dat")
}

catlib_make_geometry_representative_hillslope <- function(params) {
make_geometry_representative_hillslope <- function(params,data_paths) {
# loader raster package to open .tif files
library("raster")

# load hillslope tool
source("hillslope_method.R")

#import spatial data
flow_accum <- raster(params$flow_accumulation) # flowaccumulation Attention should be in log scale (Edit: Not working with log - Ashish)
hillslopes <- raster(params$basin) # !!! entire basin as one hillslope
elev_2_river <- raster(params$elevation2river) # elevation to !!!elev2riv_mod !!! (because of failed calculation areas, gasp are filled with values of SAGA calculation)
dist_2_river <- raster(params$distance2river) # distance to river
dem <- raster(params$dem) # digital elevation modell
aspect <- raster(params$aspect) # aspect
river_id <- raster(params$river_id) # stream link id
flow_accum <- raster(data_paths$flow_accumulation) # flowaccumulation Attention should be in log scale (Edit: Not working with log - Ashish)
hillslopes <- raster(data_paths$hillslopes) # !!! entire basin as one hillslope
elev_2_river <- raster(data_paths$elev2river) # elevation to !!!elev2riv_mod !!! (because of failed calculation areas, gasp are filled with values of SAGA calculation)
dist_2_river <- raster(data_paths$dist2river) # distance to river
dem <- raster(data_paths$filled_dem) # digital elevation modell
aspect <- raster(data_paths$aspect) # aspect
river_id <- raster(data_paths$river_id) # stream link id

# plot spatial data
system("mkdir /out/plots")
Expand Down Expand Up @@ -302,14 +53,14 @@ catlib_make_geometry_representative_hillslope <- function(params) {
"stream_id" = river_id)

# Select a hillslope from the hillslope raster map (integer value) - halfbasins file
hillslope_nr <- 1
hillslope_nr <- params$hillslope_id

# Run hillslope function. no_rf= Percentage of no flow area for region with almost no slope
# min_dist = minimum number of cells within a hillslope; freedom= freedom of spline function
pdf("/out/plots/hillslope_plots.pdf")

hill <- hillslope_tool(hillslope_nr, li_spatial, plot_hillslope_width = TRUE, plot_2d_catena = TRUE,
no_rf = 0.37, min_dist = 10, min_area = 10000, freedom = 10)
no_rf = params$no_flow_area, min_dist = params$min_cells, min_area = 10000, freedom = 10)

dev.off()

Expand All @@ -333,11 +84,11 @@ catlib_make_geometry_representative_hillslope <- function(params) {
bh = hill$short_rep_hill$short_width_corr, # Width
dist = hill$short_rep_hill$short_dist,
tot.area = hill$area, # Area
htyp = 1, # Hillslope type
dyy = 1, # Thickness of profile [m] average of drillings 2.1 m
htyp = params$hill_type, # Hillslope type
dyy = params$depth, # Thickness of profile [m] average of drillings 2.1 m
xsi = seq(0, 1, length = max(hill$short_rep_hill$short_dist) + 1), # in order to get every 1m a node length of hillslope + 1 # discretitation
eta = c(seq(0, 0.60, length = 4), seq(0.80, 1, length = 11)), # eta starts at the bottom, upper 50 cm, dx=10cm, 50-400cm dx=25 cm
out.file = "geometry.geo" # outpath
eta = c(seq(0,0.625,length=6),seq(0.7,0.85, length=3),seq(0.875,1, length=6)), # eta starts at the bottom, upper 50 cm, dx=10cm, 50-400cm dx=25 cm
out.file = "rep_hill.geo" # outpath
)

# make geometry from hillslope parameters
Expand Down
Loading

0 comments on commit 63c8c51

Please sign in to comment.