diff --git a/NAMESPACE b/NAMESPACE index 9c9f593b412bc9409a914d835c57d4d9d0967625..72481a1d7375a173b8f25c13d23630404684a49f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,7 +1,16 @@ # Generated by roxygen2: do not edit by hand export(convert_sic_params) +export(create_section_txt) +export(create_uniform_reach_txt) +export(dem_to_reach) +export(dem_to_reach_txt) +export(dem_to_section) +export(extract_reach) +export(get_section_centers) export(loadConfig) +export(merge_reaches) export(set_initial_conditions) export(sic_run_export) export(sic_run_fortran) +export(split_reach) diff --git a/R/create_section_txt.R b/R/create_section_txt.R new file mode 100644 index 0000000000000000000000000000000000000000..8f38a251869d5d95b1c4e5bff00a623df577e3d3 --- /dev/null +++ b/R/create_section_txt.R @@ -0,0 +1,43 @@ +#' Export a section in importation SIC format +#' +#' @param section_name [character], name of the section +#' @param abscissa [numeric], abscissa of the section in the reach +#' @param section_type 1-length [character], type of the section: "T" for Trapezoidal, "A" for Abscissa/Elevation, "L" for Width/Elevation +#' @param profile [list] or [matrix], profile of the section (See details) +#' @param distance_majeur [logical] or [numeric], `FALSE` for a minor bed section +#' +#' @return [character], section description in SIC text import format. +#' @export +#' @examples +#' # Trapezoidal section +#' export_section_txt("Trapeze", 0, "T", list(B = 2, S = 1.5, ZF = 100, ZB = 102)) +create_section_txt <- function(section_name, abscissa, section_type, profile, distance_majeur = FALSE) { + + if (section_type == "T") { + if (!is.list(profile) || !all(c("B", "S", "ZF", "ZB") %in% names(profile))) { + stop("With `section_type = \"T\", `profile` should be a list with the items B, S, ZF and ZB") + } + sic_profile <- c(paste(profile$B, profile$S, sep = "\t"), paste(profile$ZB, profile$ZF, sep = "\t")) + } else if (section_type %in% c("A", "L")) { + if (!is.matrix(profile) || ncol(profile) != 2) { + stop("With `section_type = \"A\" or \"L\", `profile` should be a matrix with 2 columns") + } + sic_profile <- capture.output(write.table(profile, sep = "\t", quote = FALSE, col.names = FALSE, row.names = FALSE)) + } else { + stop("section_type ", section_type, " not allowed. Possible choices are: A, L and T") + } + + if (!(is.logical(distance_majeur) || is.numeric(distance_majeur)) || length(distance_majeur) != 1) { + stop("`distance_majeur` should be a `logical` or a `numeric` of length 1") + } + bMajorBed <- !is.logical(distance_majeur) + section_txt <- c(paste(section_name, + abscissa, + ifelse(bMajorBed, distance_majeur, ""), + ifelse(bMajorBed, "1", "0"), + section_type, + sep = " $ "), + sic_profile) + class(section_txt) <- c("SectionTxt", class(section_txt)) + return(section_txt) +} diff --git a/R/create_uniform_reach_txt.R b/R/create_uniform_reach_txt.R new file mode 100644 index 0000000000000000000000000000000000000000..2594e8cfe0b117f4a41886eb5cc350f2ff6ec269 --- /dev/null +++ b/R/create_uniform_reach_txt.R @@ -0,0 +1,55 @@ +#' Generate text export geometry for an uniform reach +#' +#' @param abscissas [numeric] vector of section abscissas +#' @param upstream_bed_elevation [numeric], upstream bed elevation (m) +#' @param slope [numeric], bed slope of the reach (m/m) +#' @param names [character] vector of section names +#' @inheritParams create_section_txt +#' +#' @return A [list] from which each item is a section exported by [create_section_txt]. +#' Names of the list are the abscissas with trailing zeros for character sorting. +#' @export +#' +#' @examples +#' # Trapezoidal section +#' profT <- list(B = 2,S = 1, ZF = 100, ZB = 100 + 2) +#' +#' # Generate a reach with 2 sections at x=1000 and 2000 +#' create_uniform_reach_txt(abscissas = c(1000, 2000), +#' upstream_bed_elevation = 100, +#' slope = 0.001, +#' section_type = "T", +#' profile = profT) +create_uniform_reach_txt <- function(abscissas, + upstream_bed_elevation, + slope, + section_type, + profile, + section_names = paste0("Section x=", abscissas)) { + + sections <- lapply(seq_along(abscissas), function(i) { + x <- abscissas[i] + bed_elevation <- upstream_bed_elevation - (x - abscissas[1]) * slope + shifted_prof <- shift_profile(section_type, profile, bed_elevation) + create_section_txt(section_name = section_names[i], + abscissa = x, + section_type = section_type, + profile = shifted_prof) + }) + names(sections) <- sprintf("%08d", abscissas) + class(sections) <- c("ReachTxt", class(sections)) + return(sections) +} + +shift_profile <- function(section_type, profile, bed_elevation) { + shifted_prof <- profile + if (section_type == "T") { + shifted_prof$ZF <- bed_elevation + shifted_prof$ZB <- profile$ZB + bed_elevation - profile$ZF + } else if (section_type == "L") { + shifted_prof[,2] <- shifted_prof[,2] + bed_elevation - min(profile[, 2]) + } else { + stop("section_type ", section_type, " not allowed. Possible choices are: T and L") + } + return(shifted_prof) +} diff --git a/R/dem_to_reach.R b/R/dem_to_reach.R new file mode 100644 index 0000000000000000000000000000000000000000..173139f94bad9145c9664c17b3907d4375ce2ea6 --- /dev/null +++ b/R/dem_to_reach.R @@ -0,0 +1,129 @@ +#' Create SIC geometry sections from a DEM +#' +#' The coordinate system of `dem` should be a metric orthonormal coordinate system. +#' +#' @param dem [terra::SpatRaster] object with altitude data (m). +#' @param node_coords [matrix] with coordinates of starting and ending points of the line to browse +#' @param space_step 1-length [numeric], distance between each section (m) +#' @param section_width 1-length [numeric], width of the sections to create +#' @param nb_points 1-length [numeric], number of points to describe cross-section geometries +#' @param start 1-length [numeric], starting value for the chainage (i.e. section abscissa) along the reach +#' @param major_bed [logical], `TRUE` for major bed, `FALSE` for minor-medium bed +#' +#' @return +#' @rdname dem_to_reach +#' @export +#' +#' @examples +#' ## Inputs preparation +#' data("floodam_ead_dem") +#' dem <- terra::rast(floodam_ead_dem) +#' node_coords <- matrix(c(102550, 102550, 110000, 108000), ncol = 2) +#' space_step = 100 +#' section_width = 5000 +#' +#' ## Get section positions +#' reach_length <- as.numeric(dist(node_coords[1:2,])) +#' nb_sections <- ceiling(reach_length / space_step) + 1 +#' section_centers <- get_section_centers(node_coords, nb_sections) +#' +#' ## Create a section profile and plot it! +#' profile <- dem_to_section(dem, node_coords, section_centers[5,], section_width) +#' plot(profile) +#' +#' ## Generate cross section profiles for a reach +#' reach <- dem_to_reach(dem, node_coords, section_centers, section_width) +#' plot(reach[[5]]) +#' +#' ## Generate section profiles in SIC import text format +#' reach_txt <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) +#' reach_txt[[5]] +#' +dem_to_reach_txt <- function(dem, node_coords, space_step, section_width, nb_points = 50, start = 0, major_bed = FALSE) { + reach_length <- as.numeric(dist(node_coords[1:2,])) + nb_sections <- ceiling(reach_length / space_step) + 1 + section_centers <- get_section_centers(node_coords, nb_sections) + abscissas <- start + sapply(seq_len(nrow(section_centers)), function(i) {as.numeric(dist(section_centers[c(1,i), ]))}) + reach <- dem_to_reach(dem, node_coords, section_centers, section_width, nb_points) + reach_txt <- lapply(seq_along(reach), function(i) { + section_name <- paste0(ifelse(major_bed, "Major", "Minor"), " x=", abscissas[i]) + profile <- reach[[i]] + if (major_bed != FALSE) { + distance_majeur <- abscissas[i] - abscissas[max(i - 1, 1)] + } else { + distance_majeur <- FALSE + } + create_section_txt(section_name, abscissas[i], "A", profile, distance_majeur) + }) + names(reach_txt) <- sprintf("%08d", abscissas) + class(reach_txt) <- c("ReachTxt", class(reach_txt)) + return(reach_txt) +} + +#' @rdname dem_to_reach +#' @export +dem_to_reach <- function(dem, node_coords, section_centers, section_width, nb_points = 50) { + lapply(seq_len(nrow(section_centers)), function(i) { + section_center <- section_centers[i,] + dem_to_section(dem, node_coords, section_center, section_width, nb_points) + }) +} + + +#' Create a section cross profile from a DEM +#' +#' @inheritParams dem_to_reach +#' @param section_center 2-lenght [numeric], coordinates of the section center +#' +#' @return +#' @export +#' +#' @inherit dem_to_reach return examples +dem_to_section <- function(dem, node_coords, section_center, section_width, nb_points = 50) { + if (length(section_center) != 2) { + stop("`section_center` should be of lenght 2") + } + section_bounds <- + t(sapply(c(-1, 1), function(direction) { + get_section_bound(node_coords, section_center, section_width, direction) + })) + section_points <- + apply(section_bounds, 2, function(x) { + seq(x[1], x[2], length.out = nb_points) + }) + z <- terra::extract(dem, section_points, method = "bilinear")$lyr.1 + x_points <- sapply(seq_len(nrow(section_points)), function(i) {as.numeric(dist(section_points[c(1,i), ]))}) + m <- matrix(c(x_points, z), ncol = 2) + colnames(m) <- c("x", "z") + return(m) +} + +get_section_bound <- function(node_coords, section_center, section_width, direction = 1) { + theta_reach <- atan2(diff(node_coords[,2]), diff(node_coords[,1])) + theta <- theta_reach + direction * pi / 2 + r <- section_width / 2 + c(x = section_center[1] + r * cos(theta), y = section_center[2] + r * sin(theta)) +} + +#' Get coordinates of section centers +#' +#' From the coordinates of a segment, this function returns the coordinates +#' of equidistant points along the segment given a number of points +#' +#' @param node_coords 2-rows [matrix], with coordinates of segment boundaries +#' @param nb_sections 1-length [numeric], number of points to extract +#' +#' @return A [matrix] with the coordinates of the points along the segment +#' @export +#' +#' @examples +#' node_coords <- matrix(c(102550, 102550, 110000, 108000), ncol = 2) +#' reach_length <- as.numeric(dist(node_coords[1:2,])) +#' space_step <- 100 +#' nb_sections <- ceiling(reach_length / space_step) + 1 +#' get_section_centers(node_coords, nb_sections) +get_section_centers <- function(node_coords, nb_sections) { + apply(node_coords, 2, function(x) { + seq(x[1], x[2], length.out = nb_sections) + }) +} diff --git a/R/floodam_ead_dem.R b/R/floodam_ead_dem.R new file mode 100644 index 0000000000000000000000000000000000000000..d3bcc109ebf9d7e5d002ab1560c99a0eb1c80f24 --- /dev/null +++ b/R/floodam_ead_dem.R @@ -0,0 +1,13 @@ +#' Default DEM from the floodam-ead package +#' +#' This is the default DEM generated by [floodam.ead::generate_DEM] function with default arguments. +#' +#' @format A `PackedSpatRaster` object, use 'terra::rast()' to unpack it. +#' @source https://gitlab.irstea.fr/flood-appraisal/floodam/floodam-ead +#' +#' @examples +#' data("floodam_ead_dem", package = "rsic2") +#' floodam_ead_dem +#' terra::rast(floodam_ead_dem) +"floodam_ead_dem" + diff --git a/R/merge_reaches.R b/R/merge_reaches.R new file mode 100644 index 0000000000000000000000000000000000000000..dd3b31753f1de4b639ad7e1b8493ffa22c9ea667 --- /dev/null +++ b/R/merge_reaches.R @@ -0,0 +1,42 @@ +#' Merge several ReachTxt objects into one +#' +#' @param ... ReachTxt objects +#' +#' @return +#' @export +#' +#' @examples +#' # Minor bed generation +#' profT <- list( +#' B = 2, +#' S = (6 - 2) / 2 / 2, +#' ZF = 100, +#' ZB = 100 + 2 +#' ) +#' min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), +#' upstream_bed_elevation = 10 + 2000 * 0.002, +#' slope = 0.002, +#' section_type = "T", +#' profile = profT) +#' +#' # Major bed generation +#' data("floodam_ead_dem") +#' dem <- terra::rast(floodam_ead_dem) +#' node_coords <- matrix(c(102550, 102550, 110000, 100000), ncol = 2) +#' space_step = 100 +#' section_width = 5000 +#' maj_reach <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) +#' +#' # Merge reaches into one +#' reach <- merge_reaches(min_reach, maj_reach) +merge_reaches <- function(...) { + reaches <- list(...) + lapply(reaches, function(reach) { + if (!inherits(reach, "ReachTxt")) stop("Parameters must be of class ReachTxt") + }) + if (length(reaches) == 1) return(reaches[[1]]) + merged_reach <- do.call(c, reaches) + merged_reach <- merged_reach[order(names(merged_reach))] + class(merged_reach) <- c("ReachTxt", class(merged_reach)) + return(merged_reach) +} diff --git a/R/sic_import_reaches.R b/R/sic_import_reaches.R new file mode 100644 index 0000000000000000000000000000000000000000..d363252820e53d2b2f70e24e7f20f7145467d20f --- /dev/null +++ b/R/sic_import_reaches.R @@ -0,0 +1,33 @@ +sic_import_reaches <- function(reaches, import_mode = "ImportXml_UPDATE", cfg = loadConfig()) { + # Create reach files + import_path <- dirname(cfg$project$path) + reach_files <- sapply(seq_along(reaches), function(i) { + file <- file.path(import_path, sprintf("reach_%04d.txt", i)) + write_reach_txt(file, reaches[[i]]) + return(shQuote(basename(file), type = "cmd")) + }) + + # Import files with EdiSic https://sic.g-eau.fr/Import-sections-in-text-format?lang=en + project_name <- basename(tools::file_path_sans_ext(cfg$project$path)) + cmd_line <- shQuote( + paste( + file.path(cfg$sic$path, cfg$sic$edisic), + shQuote(dirname(cfg$project$path), type = "cmd"), + import_mode, + shQuote(project_name, type = "cmd"), + paste(reach_files, collapse = " ") + ), + type = "cmd2" + ) + logger::log_debug(cmd_line) + shell( + cmd_line, + wait = T, + translate = T + ) +} + +write_reach_txt <- function(file, reach) { + s <- unlist(reach) + writeLines(s, file) +} diff --git a/R/split_reach.R b/R/split_reach.R new file mode 100644 index 0000000000000000000000000000000000000000..e6850990aa4c84384dd99acf7b0ab4e05e33e76b --- /dev/null +++ b/R/split_reach.R @@ -0,0 +1,32 @@ +#' Split a reach into a list of reaches +#' +#' @param reach A ReachTxt object +#' @param x_limits Chainage (i.e. abscissas along the reach) of splitting points +#' +#' @return A [list] of ReachTxt objects +#' @export +#' +#' @examples +split_reach <- function(reach, x_limits) { + if(length(x_limits) < 2) { + stop("`x_limits` length must be greater or equal to 2") + } + lapply(seq_len(length(x_limits) - 1), function(i) extract_reach(reach, x_limits[i:(i+1)])) +} + +#' Select portion of a reach between two chainages +#' +#' @param reach A ReachTxt object +#' @param x_limits 2-length [numeric], min and max chainage +#' +#' @return A ReachTxt object containing the selected sections. +#' @export +#' +#' @examples +extract_reach <- function(reach, x_limits) { + reach_names <- names(reach) + x_limits <- sprintf("%08d", x_limits) + sel_reach <- reach[reach_names >= min(x_limits) & reach_names <= max(x_limits)] + class(sel_reach) <- c("ReachTxt", class(sel_reach)) + return(sel_reach) +} diff --git a/data-raw/floodam_ead_DEM.R b/data-raw/floodam_ead_DEM.R new file mode 100644 index 0000000000000000000000000000000000000000..e86c2c920fb722f586825fa03719ee7c713a163d --- /dev/null +++ b/data-raw/floodam_ead_DEM.R @@ -0,0 +1,138 @@ +## code to prepare `floodam-ead_DEM` + +# Prior to run this script, install all suggested packages with `remotes::install_deps(dependencies = TRUE)` + +# Function exported from floodam-ead 2022-02-03 +#' @title Generate floodplain DEM including streams +#' +#' @param min_elevation numeric, floodplain minimum elevation +#' @param D numeric, lateral floodplain elevation range (across floodplain) +#' @param rad numeric, x and y extent for a cell in pixels +#' @param n_y numeric, number of cells along y axis +#' @param r numeric, raster resolution in m +#' @param slope numeric, upstream-downstream slope +#' @param epsg, integer, EPSG Geodetic Parameter. Default value corresponding +#' Chaco area in Argentina (Def :UTM Zone 19S) +#' @param lower_left numeric, lower left, i.e. floodplain origin in metric +#' geodesic system +#' @param s_depth numeric, stream depth in m +#' @param s_up_width numeric, upper riverbed width in m (assuming symetric and +#' trapezoidal crosssection) +#' @param s_down_width numeric, lower riverbed width in m (assuming symetric +#' and trapezoidal crosssection) +#' @param ridge_ratio numeric, width ratio for ridge position +#' @param lambda numeric, sigmoid exponent controling bank slope +#' @param bank_width numeric, sigmoid width in pixels +#' @param random_microtopo logical, addition of random microtopo +#' @param seed integer, seed for generation of random microtopo +#' @param microtopo_sill numeric, square of macrotopo / microtopo ratio +#' (default 0.05) +#' @param microtopo_range numeric, microtopo spatial range in pixels +#' +#' @return Raster object +#' +#' @export +#' +#' @examples +#' DEM <- generate_DEM() +#' +#' @encoding UTF-8 +#' @author Jean-Stéphane BAILLY, \email{bailly@agroparistech.fr} +#' @author Frédéric Grelot, \email{frederic.grelot@inrae.fr} + +generate_DEM = function( + min_elevation = 10, + D = 20, + rad = 100, + n_y = 4, + r = 25, + slope = 0.002, + epsg = 32719, + lower_left = c(100000, 100000), + s_depth = 2, + s_up_width = 100, + s_down_width = 75, + ridge_ratio = 1 / 3, + lambda = 1, + bank_width = rad / 5, + random_microtopo = "FALSE", + seed = 123, + microtopo_sill = (0.02)^2, + microtopo_range = rad / 10 +) { + + # Cell design + x <- matrix(1:rad, nrow = rad, ncol = rad, byrow = TRUE) + y <- matrix(rad:1, nrow = rad, ncol = rad) + # ridge shoulder location along x axis (in pixels) + ridge_x <- (1 + sin(seq(pi / 2, 5 * pi / 2, length.out = rad))) / 2 + ridge_x <- round(1 + (ridge_x * ((ridge_ratio * rad)-1))) + sigm <- function(xi, lambda, bank_width) { + 1 / (1 + exp(lambda * (xi[-length(xi)] - xi[length(xi)]) / bank_width)) + } + prov <- cbind(x, ridge_x) + z0 <- t(apply(prov, 1, sigm, lambda = lambda, bank_width = bank_width)) + + # Cell replication + z0 <- cbind(z0, z0[ , ncol(z0):1]) + z0 <- do.call(rbind, mget(paste0("z", rep(0, n_y)))) + + # Adding local microtopo from spatial random realization + if (random_microtopo == "TRUE") { + set.seed(seed) + mt <- as.matrix( + RandomFields::RFsimulate( + x = 1:(rad * n_y), + y = 1:(rad * 2), + model = RandomFields::RMgauss( + var = microtopo_sill, + scale = microtopo_range + ), + grid = TRUE + ) + ) + z0 <- z0 + mt + } + + # Elevation scaling + z0 <- (z0 * D) + min_elevation + + # Stream concatenation + s <- seq(0, (s_up_width / 2), by = r) + zs <- rep(0, length(s)) + zs[s <= s_down_width / 2] <- - s_depth + zs[s > s_down_width / 2] <- + - s_depth + + 2 * s_depth * + (s[s > s_down_width / 2] - s_down_width / 2) / + (s_up_width - s_down_width) + zs <- c(rev(zs[-1]), zs) + zs <- as.matrix(do.call(rbind, mget(rep("zs", nrow(z0))))) + zs <- zs + min_elevation + z0 <- cbind( + z0[ , 1:(ncol(z0) / 2)], + zs, + z0[ , (1 + ncol(z0) / 2):ncol(z0)] + ) + + # Upstream-downstream slopping + z0 <- z0 + (nrow(z0) + 1 - row(z0)) * r * slope + + # Shape result + result <- terra::rast( + z0, + extent = terra::ext( + lower_left[1], + lower_left[1] + (r * (rad * 2 + ncol(zs))), + lower_left[2], + lower_left[2] + (r * rad * n_y) + ), + crs = sprintf("epsg:%s", epsg) + ) + + invisible(return(result)) +} + +floodam_ead_dem <- terra::wrap(generate_DEM()) + +usethis::use_data(floodam_ead_dem, overwrite = TRUE) diff --git a/data/floodam_ead_dem.rda b/data/floodam_ead_dem.rda new file mode 100644 index 0000000000000000000000000000000000000000..4b38496243e9d81dd3e1fd211b7839ddc6b3f80b Binary files /dev/null and b/data/floodam_ead_dem.rda differ diff --git a/man/create_section_txt.Rd b/man/create_section_txt.Rd new file mode 100644 index 0000000000000000000000000000000000000000..cd02ce1b2d7cc5467c1f085432653c477936510d --- /dev/null +++ b/man/create_section_txt.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/create_section_txt.R +\name{create_section_txt} +\alias{create_section_txt} +\title{Export a section in importation SIC format} +\usage{ +create_section_txt( + section_name, + abscissa, + section_type, + profile, + distance_majeur = FALSE +) +} +\arguments{ +\item{section_name}{\link{character}, name of the section} + +\item{abscissa}{\link{numeric}, abscissa of the section in the reach} + +\item{section_type}{1-length \link{character}, type of the section: "T" for Trapezoidal, "A" for Abscissa/Elevation, "L" for Width/Elevation} + +\item{profile}{\link{list} or \link{matrix}, profile of the section (See details)} + +\item{distance_majeur}{\link{logical} or \link{numeric}, \code{FALSE} for a minor bed section} +} +\value{ +\link{character}, section description in SIC text import format. +} +\description{ +Export a section in importation SIC format +} +\examples{ +# Trapezoidal section +export_section_txt("Trapeze", 0, "T", list(B = 2, S = 1.5, ZF = 100, ZB = 102)) +} diff --git a/man/create_uniform_reach_txt.Rd b/man/create_uniform_reach_txt.Rd new file mode 100644 index 0000000000000000000000000000000000000000..f47a6405f510c7dafba40ecfe942730cfaccfcda --- /dev/null +++ b/man/create_uniform_reach_txt.Rd @@ -0,0 +1,46 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/create_uniform_reach_txt.R +\name{create_uniform_reach_txt} +\alias{create_uniform_reach_txt} +\title{Generate text export geometry for an uniform reach} +\usage{ +create_uniform_reach_txt( + abscissas, + upstream_bed_elevation, + slope, + section_type, + profile, + section_names = paste0("Section x=", abscissas) +) +} +\arguments{ +\item{abscissas}{\link{numeric} vector of section abscissas} + +\item{upstream_bed_elevation}{\link{numeric}, upstream bed elevation (m)} + +\item{slope}{\link{numeric}, bed slope of the reach (m/m)} + +\item{section_type}{1-length \link{character}, type of the section: "T" for Trapezoidal, "A" for Abscissa/Elevation, "L" for Width/Elevation} + +\item{profile}{\link{list} or \link{matrix}, profile of the section (See details)} + +\item{names}{\link{character} vector of section names} +} +\value{ +A \link{list} from which each item is a section exported by \link{create_section_txt}. +Names of the list are the abscissas with trailing zeros for character sorting. +} +\description{ +Generate text export geometry for an uniform reach +} +\examples{ +# Trapezoidal section +profT <- list(B = 2,S = 1, ZF = 100, ZB = 100 + 2) + +# Generate a reach with 2 sections at x=1000 and 2000 +create_uniform_reach_txt(abscissas = c(1000, 2000), + upstream_bed_elevation = 100, + slope = 0.001, + section_type = "T", + profile = profT) +} diff --git a/man/dem_to_reach.Rd b/man/dem_to_reach.Rd new file mode 100644 index 0000000000000000000000000000000000000000..0cbfdd85a49f2bf752e49380292b36390bfd6247 --- /dev/null +++ b/man/dem_to_reach.Rd @@ -0,0 +1,66 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dem_to_reach.R +\name{dem_to_reach_txt} +\alias{dem_to_reach_txt} +\alias{dem_to_reach} +\title{Create SIC geometry sections from a DEM} +\usage{ +dem_to_reach_txt( + dem, + node_coords, + space_step, + section_width, + nb_points = 50, + start = 0, + major_bed = FALSE +) + +dem_to_reach(dem, node_coords, section_centers, section_width, nb_points = 50) +} +\arguments{ +\item{dem}{\link[terra:SpatRaster-class]{terra::SpatRaster} object with altitude data (m).} + +\item{node_coords}{\link{matrix} with coordinates of starting and ending points of the line to browse} + +\item{space_step}{1-length \link{numeric}, distance between each section (m)} + +\item{section_width}{1-length \link{numeric}, width of the sections to create} + +\item{nb_points}{1-length \link{numeric}, number of points to describe cross-section geometries} + +\item{start}{1-length \link{numeric}, starting value for the chainage (i.e. section abscissa) along the reach} + +\item{major_bed}{\link{logical}, \code{TRUE} for major bed, \code{FALSE} for minor-medium bed} +} +\value{ + +} +\description{ +The coordinate system of \code{dem} should be a metric orthonormal coordinate system. +} +\examples{ +## Inputs preparation +data("floodam_ead_dem") +dem <- terra::rast(floodam_ead_dem) +node_coords <- matrix(c(102550, 102550, 110000, 108000), ncol = 2) +space_step = 100 +section_width = 5000 + +## Get section positions +reach_length <- as.numeric(dist(node_coords[1:2,])) +nb_sections <- ceiling(reach_length / space_step) + 1 +section_centers <- get_section_centers(node_coords, nb_sections) + +## Create a section profile and plot it! +profile <- dem_to_section(dem, node_coords, section_centers[5,], section_width) +plot(profile) + +## Generate cross section profiles for a reach +reach <- dem_to_reach(dem, node_coords, section_centers, section_width) +plot(reach[[5]]) + +## Generate section profiles in SIC import text format +reach_txt <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) +reach_txt[[5]] + +} diff --git a/man/dem_to_section.Rd b/man/dem_to_section.Rd new file mode 100644 index 0000000000000000000000000000000000000000..f3dcd7ec1de8c5b77e0e323a239874e3da02572e --- /dev/null +++ b/man/dem_to_section.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dem_to_reach.R +\name{dem_to_section} +\alias{dem_to_section} +\title{Create a section cross profile from a DEM} +\usage{ +dem_to_section(dem, node_coords, section_center, section_width, nb_points = 50) +} +\arguments{ +\item{dem}{\link[terra:SpatRaster-class]{terra::SpatRaster} object with altitude data (m).} + +\item{node_coords}{\link{matrix} with coordinates of starting and ending points of the line to browse} + +\item{section_center}{2-lenght \link{numeric}, coordinates of the section center} + +\item{section_width}{1-length \link{numeric}, width of the sections to create} + +\item{nb_points}{1-length \link{numeric}, number of points to describe cross-section geometries} +} +\value{ + +} +\description{ +Create a section cross profile from a DEM +} +\examples{ +## Inputs preparation +data("floodam_ead_dem") +dem <- terra::rast(floodam_ead_dem) +node_coords <- matrix(c(102550, 102550, 110000, 108000), ncol = 2) +space_step = 100 +section_width = 5000 + +## Get section positions +reach_length <- as.numeric(dist(node_coords[1:2,])) +nb_sections <- ceiling(reach_length / space_step) + 1 +section_centers <- get_section_centers(node_coords, nb_sections) + +## Create a section profile and plot it! +profile <- dem_to_section(dem, node_coords, section_centers[5,], section_width) +plot(profile) + +## Generate cross section profiles for a reach +reach <- dem_to_reach(dem, node_coords, section_centers, section_width) +plot(reach[[5]]) + +## Generate section profiles in SIC import text format +reach_txt <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) +reach_txt[[5]] + +} diff --git a/man/extract_reach.Rd b/man/extract_reach.Rd new file mode 100644 index 0000000000000000000000000000000000000000..7b62ec508798c4676e1518f47bf286a30f3408a2 --- /dev/null +++ b/man/extract_reach.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/split_reach.R +\name{extract_reach} +\alias{extract_reach} +\title{Select portion of a reach between two chainages} +\usage{ +extract_reach(reach, x_limits) +} +\arguments{ +\item{reach}{A ReachTxt object} + +\item{x_limits}{2-length \link{numeric}, min and max chainage} +} +\value{ +A ReachTxt object containing the selected sections. +} +\description{ +Select portion of a reach between two chainages +} diff --git a/man/floodam_ead_dem.Rd b/man/floodam_ead_dem.Rd new file mode 100644 index 0000000000000000000000000000000000000000..5e30065d2936919adc5cb78657a132c7ffdb9acc --- /dev/null +++ b/man/floodam_ead_dem.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/floodam_ead_dem.R +\docType{data} +\name{floodam_ead_dem} +\alias{floodam_ead_dem} +\title{Default DEM from the floodam-ead package} +\format{ +A \code{PackedSpatRaster} object, use 'terra::rast()' to unpack it. +} +\source{ +https://gitlab.irstea.fr/flood-appraisal/floodam/floodam-ead +} +\usage{ +floodam_ead_dem +} +\description{ +This is the default DEM generated by \link[floodam.ead:generate_DEM]{floodam.ead::generate_DEM} function with default arguments. +} +\examples{ +data("floodam_ead_dem", package = "rsic2") +floodam_ead_dem +terra::rast(floodam_ead_dem) +} +\keyword{datasets} diff --git a/man/get_section_centers.Rd b/man/get_section_centers.Rd new file mode 100644 index 0000000000000000000000000000000000000000..584d51d15681b739bed06775c2557a18d5c6727d --- /dev/null +++ b/man/get_section_centers.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dem_to_reach.R +\name{get_section_centers} +\alias{get_section_centers} +\title{Get coordinates of section centers} +\usage{ +get_section_centers(node_coords, nb_sections) +} +\arguments{ +\item{node_coords}{2-rows \link{matrix}, with coordinates of segment boundaries} + +\item{nb_sections}{1-length \link{numeric}, number of points to extract} +} +\value{ +A \link{matrix} with the coordinates of the points along the segment +} +\description{ +From the coordinates of a segment, this function returns the coordinates +of equidistant points along the segment given a number of points +} +\examples{ +node_coords <- matrix(c(102550, 102550, 110000, 108000), ncol = 2) +reach_length <- as.numeric(dist(node_coords[1:2,])) +space_step <- 100 +nb_sections <- ceiling(reach_length / space_step) + 1 +get_section_centers(node_coords, nb_sections) +} diff --git a/man/merge_reaches.Rd b/man/merge_reaches.Rd new file mode 100644 index 0000000000000000000000000000000000000000..fa030e3e039bb72077bbf7428055310b2a90845e --- /dev/null +++ b/man/merge_reaches.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/merge_reaches.R +\name{merge_reaches} +\alias{merge_reaches} +\title{Merge several ReachTxt objects into one} +\usage{ +merge_reaches(...) +} +\arguments{ +\item{...}{ReachTxt objects} +} +\value{ + +} +\description{ +Merge several ReachTxt objects into one +} +\examples{ +# Minor bed generation +profT <- list( + B = 2, + S = (6 - 2) / 2 / 2, + ZF = 100, + ZB = 100 + 2 +) +min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), + upstream_bed_elevation = 10 + 2000 * 0.002, + slope = 0.002, + section_type = "T", + profile = profT) + +# Major bed generation +data("floodam_ead_dem") +dem <- terra::rast(floodam_ead_dem) +node_coords <- matrix(c(102550, 102550, 110000, 100000), ncol = 2) +space_step = 100 +section_width = 5000 +maj_reach <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) + +# Merge reaches into one +reach <- merge_reaches(min_reach, maj_reach) +} diff --git a/man/split_reach.Rd b/man/split_reach.Rd new file mode 100644 index 0000000000000000000000000000000000000000..fe18d7eca4062dbf16e9f0e11ea5d9b5f5c1f6f6 --- /dev/null +++ b/man/split_reach.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/split_reach.R +\name{split_reach} +\alias{split_reach} +\title{Split a reach into a list of reaches} +\usage{ +split_reach(reach, x_limits) +} +\arguments{ +\item{reach}{A ReachTxt object} + +\item{x_limits}{Chainage (i.e. abscissas along the reach) of splitting points} +} +\value{ +A \link{list} of ReachTxt objects +} +\description{ +Split a reach into a list of reaches +} diff --git a/tests/testthat/test-create_section_txt.R b/tests/testthat/test-create_section_txt.R new file mode 100644 index 0000000000000000000000000000000000000000..b014043213a0a9080cb58aceac6030351968d646 --- /dev/null +++ b/tests/testthat/test-create_section_txt.R @@ -0,0 +1,41 @@ +profT <- list( + B = 2, + S = (6 - 2) / 2 / 2, + ZF = 100, + ZB = 100 + 2 +) + +test_that("Trapezoidale minor section", { + expect_equal( + create_section_txt("toto", 1000, "T", profT, distance_majeur = FALSE), + c("toto $ 1000 $ $ 0 $ T", "2\t1", "102\t100") + ) + profT_wrong <- profT + profT_wrong$ZB <- NULL + expect_error( + create_section_txt("toto", 1000, "T", profT_wrong, distance_majeur = FALSE) + ) +}) + +test_that("X/Z minor section", { + expect_error( + create_section_txt("toto", 1000, "A", profT, distance_majeur = FALSE) + ) + profA <- matrix(c(0, 2, 4, 6, 102, 100, 100, 102), ncol = 2) + expect_equal( + create_section_txt("toto", 1000, "A", profA, distance_majeur = FALSE), + c("toto $ 1000 $ $ 0 $ A", "0\t102", "2\t100", "4\t100", "6\t102") + ) +}) + + +test_that("Major section", { + expect_equal( + create_section_txt("toto", 1000, "T", profT, distance_majeur = 100)[1], + "toto $ 1000 $ 100 $ 1 $ T" + ) + expect_equal( + create_section_txt("toto", 1000, "T", profT, distance_majeur = 0)[1], + "toto $ 1000 $ 0 $ 1 $ T" + ) +}) diff --git a/tests/testthat/test-create_uniform_reach_txt.R b/tests/testthat/test-create_uniform_reach_txt.R new file mode 100644 index 0000000000000000000000000000000000000000..c6eaa3442adbd6652aab15dba749d70d8e739500 --- /dev/null +++ b/tests/testthat/test-create_uniform_reach_txt.R @@ -0,0 +1,20 @@ +profT <- list( + B = 2, + S = (6 - 2) / 2 / 2, + ZF = 100, + ZB = 100 + 2 +) + +test_that("Trapezoidal minor section", { + expect_equal( + create_uniform_reach_txt(abscissas = c(1000, 2000), + upstream_bed_elevation = 100, + slope = 0.001, + section_type = "T", + profile = profT), + list( + "00001000" = c("Section x=1000 $ 1000 $ $ 0 $ T", "2\t1", "102\t100"), + "00002000" = c("Section x=2000 $ 2000 $ $ 0 $ T", "2\t1", "101\t99") + ) + ) +}) diff --git a/tests/testthat/test-dem_to_reach.R b/tests/testthat/test-dem_to_reach.R new file mode 100644 index 0000000000000000000000000000000000000000..8b137891791fe96927ad78e64b0aad7bded08bdc --- /dev/null +++ b/tests/testthat/test-dem_to_reach.R @@ -0,0 +1 @@ + diff --git a/tests/testthat/test-merge_reaches.R b/tests/testthat/test-merge_reaches.R new file mode 100644 index 0000000000000000000000000000000000000000..ac7a2a1fd6efe93bee6e855702b69b91dd1d57e8 --- /dev/null +++ b/tests/testthat/test-merge_reaches.R @@ -0,0 +1,26 @@ +# Minor bed generation +profT <- list( + B = 2, + S = (6 - 2) / 2 / 2, + ZF = 100, + ZB = 100 + 2 +) +min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), + upstream_bed_elevation = 10 + 2000 * 0.002, + slope = 0.002, + section_type = "T", + profile = profT) + +# Major bed generation +data("floodam_ead_dem") +dem <- terra::rast(floodam_ead_dem) +node_coords <- matrix(c(102550, 102550, 110000, 100000), ncol = 2) +space_step = 100 +section_width = 5000 +maj_reach <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) + +test_that("merge_reaches should work", { + reach <- merge_reaches(min_reach, maj_reach) + expect_s3_class(reach, "ReachTxt") + expect_length(reach, length(min_reach) + length(maj_reach)) +}) diff --git a/tests/testthat/test-sic_import_reaches.R b/tests/testthat/test-sic_import_reaches.R new file mode 100644 index 0000000000000000000000000000000000000000..4ac55c863d71a2a63c5ae50f29152c00da6c29d4 --- /dev/null +++ b/tests/testthat/test-sic_import_reaches.R @@ -0,0 +1,27 @@ +skip_on_ci() + +cfg <- loadLocalConfig() + +# Minor bed generation +profT <- matrix(c(2, 6, 0, 2), ncol = 2) +min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), + upstream_bed_elevation = 8 + 2000 * 0.002, + slope = 0.002, + section_type = "L", + profile = profT) + +# Major bed generation +data("floodam_ead_dem") +dem <- terra::rast(floodam_ead_dem) +node_coords <- matrix(c(102550, 102550, 110000, 100000), ncol = 2) +space_step = 100 +section_width = 5000 +maj_reach <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) + +# Merge minor and major beds and split into 2 reaches +reach <- merge_reaches(min_reach, maj_reach) +reaches <- split_reach(reach, seq(0, 10000, 5000)) + +test_that("Geometry Import works", { + sic_import_reaches(reaches, cfg = cfg) +}) diff --git a/tests/testthat/test-sic_run_export.R b/tests/testthat/test-sic_run_export.R index d2b413da70ac8de7a524058162390344d60e3fd1..45e6a50d82e1d3097ece670d28e54d6d7260e69c 100644 --- a/tests/testthat/test-sic_run_export.R +++ b/tests/testthat/test-sic_run_export.R @@ -2,7 +2,7 @@ skip_on_ci() cfg <- loadLocalConfig() -test_that("multiplication works", { +test_that("RunExport on Fluvia run works", { sic_run_fortran("fluvia", list(SCE = 1), cfg = cfg) m <- sic_run_export(scenario = 1, params = list(t = 0), cfg = cfg) expect_type(m, "double") diff --git a/tests/testthat/test-split_reach.R b/tests/testthat/test-split_reach.R new file mode 100644 index 0000000000000000000000000000000000000000..baea982190d17f9806c8ed725e7020797e2ef624 --- /dev/null +++ b/tests/testthat/test-split_reach.R @@ -0,0 +1,24 @@ +# Minor bed generation +profT <- list( + B = 2, + S = (6 - 2) / 2 / 2, + ZF = 100, + ZB = 100 + 2 +) +min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), + upstream_bed_elevation = 10 + 2000 * 0.002, + slope = 0.002, + section_type = "T", + profile = profT) + +test_that("extract_reach works", { + sel_reach <- extract_reach(min_reach, c(2000, 4000)) + expect_s3_class(sel_reach, "ReachTxt") + expect_s3_class(sel_reach[[1]], "SectionTxt") + expect_length(sel_reach, 21) +}) + +test_that("split_reach works", { + reaches <- split_reach(min_reach, seq(0, 10000, 2000)) + expect_length(reaches, 5) +})