diff --git a/NAMESPACE b/NAMESPACE index f92b891b43a4b45ae97ce3d43fb2319299bb0eae..43ec1e978ce1374a6e9251df08ad31d7293097ec 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ S3method(SicInput,POSIXt) S3method(SicInput,data.frame) S3method(SicInput,matrix) S3method(SicInput,numeric) +S3method(merge,ReachTxt) S3method(merge,SicInput) export(SicInput) export(SicLocation) @@ -19,7 +20,6 @@ export(get_result) export(get_result_tree) export(get_section_centers) export(loadConfig) -export(merge_reaches) export(read_bin_result_matrix) export(set_boundary_ZQ) export(set_initial_conditions) @@ -30,6 +30,7 @@ export(sic_run_steady) export(sic_run_unsteady) export(sic_write_par) export(split_reach) +export(update_portion_abscissas) import(magrittr) import(utils) import(xml2) diff --git a/R/dem_to_reach.R b/R/dem_to_reach.R index 135059303245093201a3f035046510082c1e5d45..f640f39c1471b8f1cb68317707ae7cf81d49d73a 100644 --- a/R/dem_to_reach.R +++ b/R/dem_to_reach.R @@ -3,7 +3,7 @@ #' 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 node_coords a 2x2 [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 @@ -66,6 +66,8 @@ dem_to_reach_txt <- function(dem, node_coords, space_step, section_width, nb_poi #' @param section_centers See return value of [get_section_centers] #' @export dem_to_reach <- function(dem, node_coords, section_centers, section_width, nb_points = 50) { + if (nrow(node_coords) != 2 || ncol(node_coords) != 2) + stop("`node_coords` should be a matrix of dimensions 2x2") 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) @@ -95,8 +97,12 @@ dem_to_section <- function(dem, node_coords, section_center, section_width, nb_p 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), ]))}) + 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) + m <- m[!is.na(m[, 2]), ] colnames(m) <- c("x", "z") return(m) } diff --git a/R/loadConfig.R b/R/loadConfig.R index 98ea65c4c6360e6a55759caf733437f1da216871..fb1e2a82e0802a3d1d0400fffb52911c92957eb1 100644 --- a/R/loadConfig.R +++ b/R/loadConfig.R @@ -15,8 +15,9 @@ #' path: Path/To/Xml/Project/File.xml #' ``` #' -#' Moreover, the `sic_path` can be defined with the environment variable "SICPATH". This setting is a default which is overwritten -#' by the path defined in the user YAML file, which is also overwritten by the one define by the `sic_path` argument. +#' Moreover, the `sic_path` can be defined with the environment variable "SICPATH". +#' This setting is a default which is overwritten by the path defined in the user YAML file, +#' which is also overwritten by the one define by the `sic_path` argument. #' #' @return A configuration as it is returned by [config::get]. #' @@ -34,6 +35,8 @@ #' - `INTERF`: default `INTERF` parameter injected in command line arguments of TALWEG, FLUVIA, SIRENE #' - `project` #' - `path`: Path to the XML project file (should be defined by `xml_path` parameter or in the user config file) +#' - `stricklers`: A vector of 2 [numeric] values representing respectively +#' the default minor and medium bed Strickler coefficients to apply with [sic_import_reaches] #' #' @export #' diff --git a/R/merge_reaches.R b/R/merge.ReachTxt.R similarity index 82% rename from R/merge_reaches.R rename to R/merge.ReachTxt.R index 16c31059b00be40e21ea0f6273e8afe072ca6148..9654bc516175da053b71ffe6dbdfbcbb5b2339b5 100644 --- a/R/merge_reaches.R +++ b/R/merge.ReachTxt.R @@ -1,6 +1,8 @@ #' Merge several *ReachTxt* objects into one #' -#' @param ... *ReachTxt* objects +#' @param x First *ReachTxt* object to merge +#' @param y Second *ReachTxt* object to merge +#' @param ... Other *ReachTxt* objects to merge #' #' @return a *ReachTxt* object (See [create_uniform_reach_txt] and [dem_to_reach]) containing the merged reaches. #' @export @@ -28,9 +30,9 @@ #' 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(...) +#' reach <- merge(min_reach, maj_reach) +merge.ReachTxt <- function(x, y = NULL, ...) { + reaches <- c(list(x, y), list(...)) lapply(reaches, function(reach) { if (!inherits(reach, "ReachTxt")) stop("Parameters must be of class ReachTxt") }) diff --git a/R/sic_import_reaches.R b/R/sic_import_reaches.R index 71c1a2e073875d77132f979fec035dbe5aee5450..bc2e98476c138fa0ca6b58daaa2818f05aaf1c06 100644 --- a/R/sic_import_reaches.R +++ b/R/sic_import_reaches.R @@ -25,7 +25,7 @@ #' 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) +#' reach <- merge(min_reach, maj_reach) #' reaches <- split_reach(reach, seq(0, 10000, 5000)) #' #' \dontrun{ @@ -55,14 +55,70 @@ sic_import_reaches <- function(reaches, import_mode = "ImportXml_UPDATE", cfg = type = "cmd2" ) logger::log_debug(cmd_line) - shell( + ret <- shell( cmd_line, wait = T, translate = T ) + update_portion_abscissas(cfg) + return(ret) } write_reach_txt <- function(file, reach) { s <- unlist(reach) writeLines(s, file) } + +#' Update abscissas of Strickler portions in a SIC model +#' +#' @description +#' As [sic_import_reaches] redefines abscissas in the reaches, and as the Stricklers +#' are defined with abscissas, geometry importation can lead to an inconsistent model. +#' +#' So the aim of this function is to reset Strickler definitions of the whole model +#' to a default value with one Strickler portion by reach with correct abscissas. +#' +#' @template param_cfg +#' @param stricklers 2-length [numeric], Strickler coefficient to apply for minor bed and medium bed respectively +#' +#' @return Use for side effect on the XML project file. +#' +#' @family sic_import_reaches +#' +#' @export +#' +#' @examples +#' \dontrun{ +#' cfg <- cfg_tmp_project() +#' update_portion_abscissas(cfg, KMin = 50) +#' +#' # Display first Strickler definitions encountered in the XML file +#' x <- read_xml(cfg$project$path) +#' xml_find_first(x, "//Stricklers") +#' } +update_portion_abscissas <- function(cfg, + stricklers = cfg$project$stricklers) { + x <- read_xml(cfg$project$path) + x_biefs <- xml_find_all(x, "//Bief") + lapply(x_biefs, function(x) { + XD <- x %>% xml_find_first(".//Liste_Sections/SectionMin") %>% xml_attr("abscisse") + XF <- x %>% xml_find_last(".//Liste_Sections/SectionMin") %>% xml_attr("abscisse") + x_stricklers <- xml_find_all(x, ".//Stricklers") + lapply(x_stricklers, function(x) { + x_portion <- xml_find_first(x, "./Portion") + xml_attr(x_portion, "XD") <- XD + xml_attr(x_portion, "XF") <- XF + xml_set_text(xml_child(x_portion, "KMin"), as.character(stricklers[1])) + xml_set_text(xml_child(x_portion, "KMoy"), as.character(stricklers[2])) + # Delete other defined portions in the reach + other_portions <- xml_siblings(x_portion) + lapply(other_portions, xml_remove) + }) + }) + write_xml(x, cfg$project$path) +} + +xml_find_last <- function(x, path) { + x_all <- xml_find_all(x, path) + return(x_all[length(x_all)]) +} diff --git a/R/sic_inputs.R b/R/sic_inputs.R index f1c2d0bb724ce0de4434433374d45a704bdcf23b..dc6ee9e3731295f0f1e802e440bfe21590089836 100644 --- a/R/sic_inputs.R +++ b/R/sic_inputs.R @@ -49,7 +49,7 @@ sic_write_par <- function(cfg, scenario, sicInputs) { sParFileName = sic_get_par_filename(cfg, scenario) df = data.frame(C1 = "// PAR file for SIC", - C2 = "automatically generated by rsci2", + C2 = "automatically generated by rsic2", C3 = "", stringsAsFactors = FALSE) i = 0 diff --git a/R/sic_run_fortran.R b/R/sic_run_fortran.R index 759dd529b59ef59862cf2112666eeba6e090063d..8a658507d5d729a5eb75be1c84653c35cc9f762b 100644 --- a/R/sic_run_fortran.R +++ b/R/sic_run_fortran.R @@ -112,7 +112,10 @@ sic_run_steady <- function(cfg, scenario, variant = 0, sicInputs = NULL, params #' @export sic_run_unsteady <- function(cfg, scenario = iniParams[4], variant = iniParams[5], sicInputs = NULL, iniParams = NULL, params = list()) { if (!is.null(iniParams)) { - sic_run_steady(cfg, scenario = iniParams[1], variant = iniParams[1]) + sic_run_steady(cfg, + scenario = iniParams[1], + variant = iniParams[2], + sicInputs = sicInputs) set_initial_conditions(iniParams, cfg) } if (is.null(scenario)) stop("`scenario` should be defined") diff --git a/R/split_reach.R b/R/split_reach.R index ce02a351cc284107f447cbaa60d79c12757bd5b4..2f8aad137abb8d750c927926919ed7e75bd8b5ad 100644 --- a/R/split_reach.R +++ b/R/split_reach.R @@ -59,6 +59,14 @@ 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)] + # Change major bed distance of the last section (issue #13) + i_last <- length(sel_reach) + s_cfg <- strsplit(sel_reach[[i_last]][1], "$", fixed = TRUE)[[1]] + if (s_cfg[4] == " 1 ") { + s_cfg[3] = " 0 " + sel_reach[[i_last]][1] <- paste(s_cfg, collapse = "$") + } + class(sel_reach) <- c("ReachTxt", class(sel_reach)) return(sel_reach) } diff --git a/inst/config.yml b/inst/config.yml index ebdfcb2ef40a6684f722607bdfe12731c7a64598..c19fd43c67d466dce2f750895b426b212ed89783 100644 --- a/inst/config.yml +++ b/inst/config.yml @@ -11,3 +11,4 @@ default: INTERF: "0" project: path: "*** The path of your XML project file should be defined ***" + stricklers: [40, 40] diff --git a/man/dem_to_reach.Rd b/man/dem_to_reach.Rd index 828b099c70a7ecf01245bb6c856e02324ea55c20..562b123b4637d9e2c301a704abd47b7256c96382 100644 --- a/man/dem_to_reach.Rd +++ b/man/dem_to_reach.Rd @@ -20,7 +20,7 @@ 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{node_coords}{a 2x2 \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)} diff --git a/man/dem_to_section.Rd b/man/dem_to_section.Rd index 8346460fb318ba141c7abf1f7223d5a3972f2561..ab603a0d22082353a394f9bc7cd68ff7695a40b7 100644 --- a/man/dem_to_section.Rd +++ b/man/dem_to_section.Rd @@ -9,7 +9,7 @@ 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{node_coords}{a 2x2 \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} diff --git a/man/loadConfig.Rd b/man/loadConfig.Rd index 0e76006a0e91aa12876b8805587c121c87d6cc94..3307685d9a48ae6a87d2595cd974557080be14fe 100644 --- a/man/loadConfig.Rd +++ b/man/loadConfig.Rd @@ -44,6 +44,8 @@ Configuration of RSIC2 as the following structure: \item \code{project} \itemize{ \item \code{path}: Path to the XML project file (should be defined by \code{xml_path} parameter or in the user config file) +\item \code{stricklers}: A vector of 2 \link{numeric} values representing respectively +the default minor and medium bed Strickler coefficients to apply with \link{sic_import_reaches} } } } @@ -57,8 +59,9 @@ project: path: Path/To/Xml/Project/File.xml }\if{html}{\out{</div>}} -Moreover, the \code{sic_path} can be defined with the environment variable "SICPATH". This setting is a default which is overwritten -by the path defined in the user YAML file, which is also overwritten by the one define by the \code{sic_path} argument. +Moreover, the \code{sic_path} can be defined with the environment variable "SICPATH". +This setting is a default which is overwritten by the path defined in the user YAML file, +which is also overwritten by the one define by the \code{sic_path} argument. } \examples{ library(rsic2) diff --git a/man/merge_reaches.Rd b/man/merge.ReachTxt.Rd similarity index 75% rename from man/merge_reaches.Rd rename to man/merge.ReachTxt.Rd index 1ca02aa21a10fea5b1cb2d1ccae56a057a79cfae..371978ea89bd1943b40a31ebdb48eff62e660af9 100644 --- a/man/merge_reaches.Rd +++ b/man/merge.ReachTxt.Rd @@ -1,13 +1,17 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/merge_reaches.R -\name{merge_reaches} -\alias{merge_reaches} +% Please edit documentation in R/merge.ReachTxt.R +\name{merge.ReachTxt} +\alias{merge.ReachTxt} \title{Merge several \emph{ReachTxt} objects into one} \usage{ -merge_reaches(...) +\method{merge}{ReachTxt}(x, y = NULL, ...) } \arguments{ -\item{...}{\emph{ReachTxt} objects} +\item{x}{First \emph{ReachTxt} object to merge} + +\item{y}{Second \emph{ReachTxt} object to merge} + +\item{...}{Other \emph{ReachTxt} objects to merge} } \value{ a \emph{ReachTxt} object (See \link{create_uniform_reach_txt} and \link{dem_to_reach}) containing the merged reaches. @@ -38,5 +42,5 @@ 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) +reach <- merge(min_reach, maj_reach) } diff --git a/man/sic_import_reaches.Rd b/man/sic_import_reaches.Rd index 711ba97909dbc08e63b1572d6cb7c3cfa2938a9e..ff535fc8d366576ace596dbbe24da36b8ee0212a 100644 --- a/man/sic_import_reaches.Rd +++ b/man/sic_import_reaches.Rd @@ -41,7 +41,7 @@ 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) +reach <- merge(min_reach, maj_reach) reaches <- split_reach(reach, seq(0, 10000, 5000)) \dontrun{ diff --git a/man/update_portion_abscissas.Rd b/man/update_portion_abscissas.Rd new file mode 100644 index 0000000000000000000000000000000000000000..8ebc2a4fa92ce5bcfc2b8dde0d29a315d71aae78 --- /dev/null +++ b/man/update_portion_abscissas.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/sic_import_reaches.R +\name{update_portion_abscissas} +\alias{update_portion_abscissas} +\title{Update abscissas of Strickler portions in a SIC model} +\usage{ +update_portion_abscissas(cfg, stricklers = cfg$project$stricklers) +} +\arguments{ +\item{cfg}{a \link{config} object. Configuration to use. See \link{loadConfig} for details} + +\item{stricklers}{2-length \link{numeric}, Strickler coefficient to apply for minor bed and medium bed respectively} +} +\value{ +Use for side effect on the XML project file. +} +\description{ +As \link{sic_import_reaches} redefines abscissas in the reaches, and as the Stricklers +are defined with abscissas, geometry importation can lead to an inconsistent model. + +So the aim of this function is to reset Strickler definitions of the whole model +to a default value with one Strickler portion by reach with correct abscissas. +} +\examples{ +\dontrun{ +cfg <- cfg_tmp_project() +update_portion_abscissas(cfg, KMin = 50) + +# Display first Strickler definitions encountered in the XML file +x <- read_xml(cfg$project$path) +xml_find_first(x, "//Stricklers") +} +} +\concept{sic_import_reaches} diff --git a/tests/testthat/test-dem_to_reach.R b/tests/testthat/test-dem_to_reach.R index 8b137891791fe96927ad78e64b0aad7bded08bdc..cd15173ee79c5011a182cdad3d8706eaa877414c 100644 --- a/tests/testthat/test-dem_to_reach.R +++ b/tests/testthat/test-dem_to_reach.R @@ -1 +1,14 @@ +data("floodam_ead_dem") +dem <- terra::rast(floodam_ead_dem) +node_coords <- matrix(c(102550, 102550, 110000, 100000), ncol = 2) + +extent <- terra::ext(dem) +maj_section_width <- as.numeric(extent$xmax - extent$xmin) + +nb_sections <- ceiling(10000 / 25) + 1 +section_centers <- get_section_centers(node_coords, nb_sections) + +test_that("dem_to_reach works",{ + reach <- dem_to_reach(dem, node_coords, section_centers, section_width = maj_section_width) +}) diff --git a/tests/testthat/test-merge_reaches.R b/tests/testthat/test-merge.ReachTxt.R similarity index 86% rename from tests/testthat/test-merge_reaches.R rename to tests/testthat/test-merge.ReachTxt.R index ac7a2a1fd6efe93bee6e855702b69b91dd1d57e8..fa3b29f0ee8c1a556cf36f0b84222327abd83206 100644 --- a/tests/testthat/test-merge_reaches.R +++ b/tests/testthat/test-merge.ReachTxt.R @@ -5,7 +5,7 @@ profT <- list( ZF = 100, ZB = 100 + 2 ) -min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), +min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 500), upstream_bed_elevation = 10 + 2000 * 0.002, slope = 0.002, section_type = "T", @@ -15,12 +15,12 @@ min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), data("floodam_ead_dem") dem <- terra::rast(floodam_ead_dem) node_coords <- matrix(c(102550, 102550, 110000, 100000), ncol = 2) -space_step = 100 +space_step = 500 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) +test_that("merge.ReachTxt should work", { + reach <- merge(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 index d0ea19db98f1a3ed811f28322a3cfb0112d423e7..56d20c2be99bbd8c6368725ecf0045c93b15d1ee 100644 --- a/tests/testthat/test-sic_import_reaches.R +++ b/tests/testthat/test-sic_import_reaches.R @@ -19,7 +19,7 @@ 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) +reach <- merge(min_reach, maj_reach) reaches <- split_reach(reach, seq(0, 10000, 5000)) test_that("Geometry Import works", { diff --git a/tests/testthat/test-split_reach.R b/tests/testthat/test-split_reach.R index baea982190d17f9806c8ed725e7020797e2ef624..9947aa2427f377d418c3f16b612d2590b15a4359 100644 --- a/tests/testthat/test-split_reach.R +++ b/tests/testthat/test-split_reach.R @@ -5,7 +5,7 @@ profT <- list( ZF = 100, ZB = 100 + 2 ) -min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 100), +min_reach <- create_uniform_reach_txt(abscissas = seq(0, 10000, 500), upstream_bed_elevation = 10 + 2000 * 0.002, slope = 0.002, section_type = "T", @@ -15,10 +15,24 @@ 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) + expect_length(sel_reach, 5) }) test_that("split_reach works", { reaches <- split_reach(min_reach, seq(0, 10000, 2000)) expect_length(reaches, 5) }) + +test_that("Major bed section should have correct lengths", { + # 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 = 25 + section_width = 5000 + maj_reach <- dem_to_reach_txt(dem, node_coords, space_step, section_width, major_bed = TRUE) + reaches <- split_reach(maj_reach, seq(0, 10000, 5000)) + sn <- reaches[[1]][[length(reaches[[1]])]] + s_cfg <- strsplit(sn[1], "$", fixed = TRUE)[[1]] + expect_equal(s_cfg[3], " 0 ") +})