Commit d867f239 authored by Dorchies David's avatar Dorchies David
Browse files

Merge branch '13-crash-de-talweg-sur-les-distances-de-lit-majeur' into 'main'

Resolve "Crash de Talweg sur les distances de lit majeur"

Closes #14 and #13

See merge request !10
1 merge request!10Resolve "Crash de Talweg sur les distances de lit majeur"
Pipeline #34244 passed with stage
in 2 minutes and 36 seconds
Showing with 179 additions and 31 deletions
+179 -31
......@@ -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)
......
......@@ -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)
}
......
......@@ -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
#'
......
#' 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")
})
......
......@@ -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)])
}
......@@ -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
......
......@@ -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")
......
......@@ -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)
}
......@@ -11,3 +11,4 @@ default:
INTERF: "0"
project:
path: "*** The path of your XML project file should be defined ***"
stricklers: [40, 40]
......@@ -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)}
......
......@@ -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}
......
......@@ -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)
......
% 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)
}
......@@ -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{
......
% 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}
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)
})
......@@ -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))
})
......@@ -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", {
......
......@@ -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 ")
})
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment