From a855aaf57c82821564ef83f39183c32862c712db Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Tue, 8 Oct 2024 14:45:20 +0200 Subject: [PATCH 01/25] Removed citation of procrust in cutPlot function and replaced it a citation of bilinear interpolation. Added descriptions of last 3 pull-requests in the NEWS file. Changed the version from 2.1.13 to 2.1.14 in the description file, and added myself as an author --- DESCRIPTION | 3 ++- NEWS | 5 +++++ R/cutPlot.R | 2 +- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 3c0f039..effbb13 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: BIOMASS Type: Package Title: Estimating Aboveground Biomass and Its Uncertainty in Tropical Forests -Version: 2.1.13 +Version: 2.1.14 Date: 2024-03-19 Authors@R: c( person("Maxime", "Réjou-Méchain", email = "maxime.rejou@gmail.com", role = c("aut", "dtc")), @@ -14,6 +14,7 @@ Authors@R: c( person("Bruno", "Hérault", email = "bruno.herault@cirad.fr", role = c("aut")), person("Ted", "Feldpausch", email = "T.R.Feldpausch@exeter.ac.uk", role = c("dtc")), person("Philippe", "Verley", email = "philippe.verley@ird.fr ", role = c("ctb")) + person("Arthur", "Bailly", email = "arthur.bailly@ird.fr ", role = c("aut")) ) Description: Contains functions to estimate aboveground biomass/carbon and its uncertainty in tropical forests. These functions allow to (1) retrieve and to correct taxonomy, (2) estimate wood density and its uncertainty, diff --git a/NEWS b/NEWS index 119dd0f..a0a5a42 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,8 @@ +* version 2.1.14 +- Fix an issue about the inversion of subplot locations in cutPlot +- Change the automatic corner numbering from counter-clockwise to clockise direction in cutPlot +- Fix an issue with summaryByPlot when a plot do not contains any tree + * version 2.1.13 - Activate pkgdown diff --git a/R/cutPlot.R b/R/cutPlot.R index dbd1e3f..d43ef25 100644 --- a/R/cutPlot.R +++ b/R/cutPlot.R @@ -7,7 +7,7 @@ if (getRversion() >= "2.15.1") { #' Divides a plot in subplots #' #' This function divides a plot (or several plots) in subplots and returns the coordinates of the grid. -#' This function uses a procrust analysis to fit the rectangle you gave to the plot you have. +#' These coordinates are calculated by a bilinear interpolation with the projected corner coordinates as references. #' #' @param projCoord A data frame containing the projected coordinates of plot corners, with X and Y on the first and second column respectively #' @param plot A vector indicating the plot codes From dd39056a92a597450e125810c3c0b3faa0385349 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Tue, 8 Oct 2024 14:49:00 +0200 Subject: [PATCH 02/25] added a missing comma --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index effbb13..7bd5294 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,7 +13,7 @@ Authors@R: c( person("Jerome", "Chave", email = "jerome.chave@univ-tlse3.fr", role = c("dtc")), person("Bruno", "Hérault", email = "bruno.herault@cirad.fr", role = c("aut")), person("Ted", "Feldpausch", email = "T.R.Feldpausch@exeter.ac.uk", role = c("dtc")), - person("Philippe", "Verley", email = "philippe.verley@ird.fr ", role = c("ctb")) + person("Philippe", "Verley", email = "philippe.verley@ird.fr ", role = c("ctb")), person("Arthur", "Bailly", email = "arthur.bailly@ird.fr ", role = c("aut")) ) Description: Contains functions to estimate aboveground biomass/carbon and its uncertainty in tropical forests. From b63a24ebc8b7a24efe3f8a94b853b166de504f20 Mon Sep 17 00:00:00 2001 From: Guillaume Cornu <21357644+billy34@users.noreply.github.com> Date: Wed, 9 Oct 2024 09:57:48 +0200 Subject: [PATCH 03/25] Update DESCRIPTION missing comma in authors --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index effbb13..7bd5294 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -13,7 +13,7 @@ Authors@R: c( person("Jerome", "Chave", email = "jerome.chave@univ-tlse3.fr", role = c("dtc")), person("Bruno", "Hérault", email = "bruno.herault@cirad.fr", role = c("aut")), person("Ted", "Feldpausch", email = "T.R.Feldpausch@exeter.ac.uk", role = c("dtc")), - person("Philippe", "Verley", email = "philippe.verley@ird.fr ", role = c("ctb")) + person("Philippe", "Verley", email = "philippe.verley@ird.fr ", role = c("ctb")), person("Arthur", "Bailly", email = "arthur.bailly@ird.fr ", role = c("aut")) ) Description: Contains functions to estimate aboveground biomass/carbon and its uncertainty in tropical forests. From 9d214bd06a1fd0716b348929c109b7f4eabc7bf9 Mon Sep 17 00:00:00 2001 From: Guillaume Cornu Date: Wed, 9 Oct 2024 15:02:52 +0200 Subject: [PATCH 04/25] allow multiple spaces between genus and species "Terminalia superba" was ill detected because there was 2 spaces --- R/correctTaxo.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/correctTaxo.R b/R/correctTaxo.R index 228b6c5..5220a77 100644 --- a/R/correctTaxo.R +++ b/R/correctTaxo.R @@ -117,7 +117,7 @@ correctTaxo <- function(genus, species = NULL, score = 0.5, useCache = FALSE, ve # sub-function definition ------------------------------------------------- # split x always returning count columns (padding with NA) - tstrsplit_NA <- function(x, pattern = " ", count = 2) { + tstrsplit_NA <- function(x, pattern = "\\s+", count = 2) { # NOTE extraneous columns ignored maybe better paste them together split <- utils::head(tstrsplit(x, pattern), count) From 354ca077c6dc5357c6c4fd9d17e5c298f27682e3 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 9 Oct 2024 15:16:32 +0200 Subject: [PATCH 05/25] Added the use of bilinear interpolation in the NEWS file --- NEWS | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS b/NEWS index a0a5a42..df7e76f 100644 --- a/NEWS +++ b/NEWS @@ -2,6 +2,7 @@ - Fix an issue about the inversion of subplot locations in cutPlot - Change the automatic corner numbering from counter-clockwise to clockise direction in cutPlot - Fix an issue with summaryByPlot when a plot do not contains any tree +- Replace the procrust analyses used to calculate sub-plot corners with a bilinear interpolation that takes the coordinates of the projected plot corners as references. * version 2.1.13 - Activate pkgdown From 96756a0290763a5ed2c059457198353dc0fe659c Mon Sep 17 00:00:00 2001 From: Guillaume Cornu Date: Wed, 9 Oct 2024 15:28:23 +0200 Subject: [PATCH 06/25] sanity cleaning of queried species names --- R/correctTaxo.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/correctTaxo.R b/R/correctTaxo.R index 5220a77..32f05b6 100644 --- a/R/correctTaxo.R +++ b/R/correctTaxo.R @@ -130,7 +130,7 @@ correctTaxo <- function(genus, species = NULL, score = 0.5, useCache = FALSE, ve # Data preparation -------------------------------------------------------- - genus <- as.character(genus) + genus <- stringr::str_squish(as.character(genus)) if (is.null(species)) { @@ -143,7 +143,7 @@ correctTaxo <- function(genus, species = NULL, score = 0.5, useCache = FALSE, ve # split genus (query) userTaxo[, c("genus", "species") := tstrsplit_NA(query)] } else { - species <- as.character(species) + species <- stringr::str_squish(as.character(species)) # Create a dataframe with the original values userTaxo <- data.table( @@ -156,10 +156,10 @@ correctTaxo <- function(genus, species = NULL, score = 0.5, useCache = FALSE, ve } # If there is an empty genus - userTaxo[genus == "", ":="(genus = NA_character_, species = NA_character_, query = NA_character_)] + userTaxo[is.na(genus) | (genus == ""), ":="(genus = NA_character_, species = NA_character_, query = NA_character_)] # If there is empty species - userTaxo[species == "", ":="(species = NA_character_, query = gsub(" ", "", query))] + userTaxo[is.na(species) | (species == ""), ":="(species = NA_character_, query = gsub(" ", "", query))] # get unique values qryTaxo <- unique(userTaxo[!is.na(query)]) From 5fd83b40f0a3093a8c530ae4b69d2bc2edf7440a Mon Sep 17 00:00:00 2001 From: Guillaume Cornu Date: Wed, 9 Oct 2024 16:50:38 +0200 Subject: [PATCH 07/25] remove dependency on stringr --- R/correctTaxo.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/R/correctTaxo.R b/R/correctTaxo.R index 32f05b6..36069a0 100644 --- a/R/correctTaxo.R +++ b/R/correctTaxo.R @@ -130,7 +130,13 @@ correctTaxo <- function(genus, species = NULL, score = 0.5, useCache = FALSE, ve # Data preparation -------------------------------------------------------- - genus <- stringr::str_squish(as.character(genus)) + # remove spaces at beginning and end, and remove extra spacing + squish <- function(x) { + x <- gsub("(^\\s+)|(\\s+$)", "", x) + gsub("\\s+", " ", x) + } + + genus <- squish(as.character(genus)) if (is.null(species)) { @@ -143,7 +149,7 @@ correctTaxo <- function(genus, species = NULL, score = 0.5, useCache = FALSE, ve # split genus (query) userTaxo[, c("genus", "species") := tstrsplit_NA(query)] } else { - species <- stringr::str_squish(as.character(species)) + species <- squish(as.character(species)) # Create a dataframe with the original values userTaxo <- data.table( From 378a4ff56ff165f65968c3b5029e03182e384b33 Mon Sep 17 00:00:00 2001 From: lamonica-d <104020730+lamonica-d@users.noreply.github.com> Date: Tue, 15 Oct 2024 10:41:14 +0200 Subject: [PATCH 08/25] Create test-coverage.yaml --- .github/workflows/test-coverage.yaml | 61 ++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 .github/workflows/test-coverage.yaml diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml new file mode 100644 index 0000000..9882260 --- /dev/null +++ b/.github/workflows/test-coverage.yaml @@ -0,0 +1,61 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + +name: test-coverage.yaml + +permissions: read-all + +jobs: + test-coverage: + runs-on: ubuntu-latest + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::covr, any::xml2 + needs: coverage + + - name: Test coverage + run: | + cov <- covr::package_coverage( + quiet = FALSE, + clean = FALSE, + install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") + ) + covr::to_cobertura(cov) + shell: Rscript {0} + + - uses: codecov/codecov-action@v4 + with: + fail_ci_if_error: ${{ github.event_name != 'pull_request' && true || false }} + file: ./cobertura.xml + plugin: noop + disable_search: true + token: ${{ secrets.CODECOV_TOKEN }} + + - name: Show testthat output + if: always() + run: | + ## -------------------------------------------------------------------- + find '${{ runner.temp }}/package' -name 'testthat.Rout*' -exec cat '{}' \; || true + shell: bash + + - name: Upload test results + if: failure() + uses: actions/upload-artifact@v4 + with: + name: coverage-test-failures + path: ${{ runner.temp }}/package From 6fed132da174e90b8d71f5f224cfa288672efc98 Mon Sep 17 00:00:00 2001 From: lamonica-d <104020730+lamonica-d@users.noreply.github.com> Date: Tue, 15 Oct 2024 15:43:14 +0200 Subject: [PATCH 09/25] Update test-coverage.yaml --- .github/workflows/test-coverage.yaml | 9 --------- 1 file changed, 9 deletions(-) diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 9882260..e5049ca 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -35,17 +35,8 @@ jobs: clean = FALSE, install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") ) - covr::to_cobertura(cov) shell: Rscript {0} - - uses: codecov/codecov-action@v4 - with: - fail_ci_if_error: ${{ github.event_name != 'pull_request' && true || false }} - file: ./cobertura.xml - plugin: noop - disable_search: true - token: ${{ secrets.CODECOV_TOKEN }} - - name: Show testthat output if: always() run: | From 557fa0e29709266b5170a5043fd1e13f880f1a7b Mon Sep 17 00:00:00 2001 From: lamonica-d <104020730+lamonica-d@users.noreply.github.com> Date: Tue, 15 Oct 2024 15:54:49 +0200 Subject: [PATCH 10/25] Update test-coverage.yaml --- .github/workflows/test-coverage.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index e5049ca..4d4a04a 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -30,7 +30,7 @@ jobs: - name: Test coverage run: | - cov <- covr::package_coverage( + covr::package_coverage( quiet = FALSE, clean = FALSE, install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") From e9cf1479d5c15fec1a75e58a1aae073d6d39300e Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Fri, 18 Oct 2024 19:05:16 +0200 Subject: [PATCH 11/25] Transformed plot to ggplot (and return as an element of output list) for the HDmodel function. The function now display every model plot fit when drawGraph is TRUE, even if there are multiple plots. --- R/modelHD.R | 158 ++++++++++++++++++------------------------ vignettes/BIOMASS.Rmd | 8 ++- 2 files changed, 74 insertions(+), 92 deletions(-) diff --git a/R/modelHD.R b/R/modelHD.R index 1ecb0ac..58c365c 100644 --- a/R/modelHD.R +++ b/R/modelHD.R @@ -2,10 +2,8 @@ #' #' This function fits and compares (optional) height-diameter models. #' -#' @param D Vector with diameter measurements (in cm). NA values are accepted but a -#' minimum of 10 valid entries (i.e. having a corresponding height in H) is required. -#' @param H Vector with total height measurements (in m). NA values are accepted but a minimum of 10 valid entries -#' (i.e. having a corresponding diameter in D) is required. +#' @param D Vector with diameter measurements (in cm). NA values are accepted but a minimum of 10 valid entries (i.e. having a corresponding height in H) is required. +#' @param H Vector with total height measurements (in m). NA values are accepted but a minimum of 10 valid entries (i.e. having a corresponding diameter in D) is required. #' @param method Method used to fit the relationship. #' To be chosen between: #' - log1, log2 @@ -26,9 +24,10 @@ #' #' #' @return -#' If plot is NULL or has a single value, a single list is returned. If there is more than one plot, +#' If `plot` is NULL or has a single value, a single list is returned. If there is more than one plot, #' multiple embedded lists are returned with plots as the list names. -#' Returns a list if the parameter model is not null: +#' +#' If `model` is not null (model comparison), returns a list : #' - `input`: list of the data used to construct the model (list(H, D)) #' - `model`: outputs of the model (same outputs as given by [stats::lm()], [stats::nls()]) #' - `RSE`: Residual Standard Error of the model @@ -39,19 +38,20 @@ #' - `formula`: Formula of the model #' - `method`: Name of the method used to construct the model #' - `predicted`: Predicted height values +#' - `fitPlot`: a ggplot object containing the model fitting plot #' #' -#' If the parameter model is null, the function return a graph with all the methods for +#' If the parameter model is null, the function return a plot with all the methods for #' comparison, the function also returns a data.frame with: -#' - `method`: The method that had been used to construct the graph -#' - `color`: The color of the curve in the graph +#' - `method`: The method that had been used to construct the plot +#' - `color`: The color of the curve in the plot #' - `RSE`: Residual Standard Error of the model #' - `RSElog`: Residual Standard Error of the log model (`NULL` if other model) #' - `Average_bias`: The average bias for the model #' #' #' -#' @author Maxime REJOU-MECHAIN, Arthur PERE, Ariane TANGUY +#' @author Maxime REJOU-MECHAIN, Arthur PERE, Ariane TANGUY, Arthur Bailly #' @seealso [retrieveH()] #' #' @export @@ -79,13 +79,10 @@ #' @importFrom stats SSmicmen lm median na.omit quantile rnorm sd predict coef #' @importFrom utils data #' @importFrom data.table data.table +#' @importFrom ggplot2 ggplot aes geom_point geom_smooth labs theme_minimal theme scale_x_continuous scale_y_continuous modelHD <- function(D, H, method = NULL, useWeight = FALSE, drawGraph = FALSE, plot = NULL) { - # # To maintain user's original options - oldpar <- par(no.readonly = TRUE) - on.exit(par(oldpar)) - # parameters verification ------------------------------------------------- # Check if there is enough data to compute an accurate model @@ -93,62 +90,55 @@ modelHD <- function(D, H, method = NULL, useWeight = FALSE, drawGraph = FALSE, p if (nbNonNA < 15) { stop("The data has not enough height data (less than 15 non NA)") } - if (length(H) != length(D)) { stop("Your vector D and H do not have the same length") } - if (!is.null(method)) { method <- tolower(method) } - methods <- c("log1", "log2", "weibull", "michaelis") if (!is.null(method) && !(method %in% methods)) { stop("Chose your method among those ones : ", paste(methods, collapse = ", ")) } - if (!is.logical(useWeight)) { stop("UseWeight argument must be a boolean") } - if (!is.logical(drawGraph)) { stop("drawGraph argument must be a boolean") } - if (!is.null(plot) && !length(plot) %in% c(1, length(D))) { stop("The length of the 'plot' vector must be either 1 or the length of D") } - - - if (!is.null(plot)) { - drawGraph <- FALSE - } - - - - # If there is a plot ID --------------------------------------------------- - + # if (!is.null(plot)) { + # drawGraph <- FALSE + # } + + + # Multiple plots managment --------------------------------------------------- + # if there is multiple plots in the plot vector if (!is.null(plot) && length(unique(plot)) != 1) { Hdata <- data.table(H = H, D = D, plot = plot) - + output <- lapply(split(Hdata, by = "plot", keep.by = TRUE), function(subData) { suppressMessages(modelHD( subData$D, subData$H, method, useWeight, drawGraph, unique(subData$plot) )) }) - + if (is.null(method)) { message("To build a HD model you must use the parameter 'method' in this function") } - + return(output) } + + # functions ---------------------------------------------------------------- - # fonction to choose the fonction + # function to select the method modSelect <- function(Hdata, method, useGraph = FALSE) { output <- list() @@ -158,7 +148,6 @@ modelHD <- function(D, H, method = NULL, useWeight = FALSE, drawGraph = FALSE, p mod <- loglogFunction(Hdata, method) output$RSElog <- summary(mod)$sigma - # Baskerville correction 1972 output$Hpredict <- exp(predict(mod) + 0.5 * output$RSElog^2) @@ -189,65 +178,45 @@ modelHD <- function(D, H, method = NULL, useWeight = FALSE, drawGraph = FALSE, p output$residuals <- res output$mod <- mod - return(output) } - # function to draw the beining of the graph drawPlotBegin <- function(givenMethod = FALSE, plotId) { main_title <- ifelse(givenMethod == FALSE, "Model comparison", paste("Selected model : ", givenMethod)) - - par(mar = c(5, 5, 3, 3)) - plot(Hdata$D, Hdata$H, - pch = 20, cex = 0.5, col = "grey50", log = "xy", las = 1, - xlab = "D (cm)", ylab = "H (m)", cex.lab = 1.8, cex.axis = 1.5, - main = main_title, - cex.main = 2, axes = FALSE, frame.plot = FALSE - ) - grid(col = "grey80", lty = 1, equilogs = FALSE) - axis(side = 1, lty = "blank", las = 1) - axis(side = 2, lty = "blank", las = 1) + + starting_plot <- ggplot(data = Hdata, mapping = aes(x=D, y=H)) + + geom_point(col="grey50") + + labs(title = main_title, x="D (cm)", y="H (m)") + + #coord_trans(x = "log", y = "log") + theme_minimal() + scale_x_continuous(transform = "log10", n.breaks = 8, )+ scale_y_continuous(transform = "log10", n.breaks = 8) + + theme_minimal() + theme(plot.title = element_text(hjust = 0.5)) + return(starting_plot) } # Data processing --------------------------------------------------------- - Hdata <- data.table(H = H, D = D) Hdata <- na.omit(Hdata) # Remove NA values weight <- NULL # Vector of diameter used only for visualisation purpose - D_Plot <- data.frame(D = Hdata[, seq(floor(min(D)), ceiling(max(D)), 0.5)]) - + #D_Plot <- data.frame(D = Hdata[, seq(floor(min(D)), ceiling(max(D)), 0.5)]) + D_Plot <- data.frame(D = Hdata[, 10^seq(log10(floor(min(D))), log10(ceiling(max(D))), l=100 )]) + # If the measures need to be weighted if (useWeight == TRUE) { weight <- Hdata[, D^2 * H] } # weight is proportional to tree volume - - # If we gave the function a method ---------------------------------------- - + # If one method is supplied ---------------------------------------------- if (!is.null(method)) { output <- modSelect(Hdata, method, drawGraph) - ####### if drawGraph is true - if (drawGraph) { - drawPlotBegin(method, plot) - - lines(D_Plot$D, output$Hpredict_plot, lwd = 2, col = "blue") - legend("bottomright", c("Data", "Model selected"), - lty = c(3, 1), lwd = c(3, 3), - col = c("grey", "blue"), cex = 1.5 - ) - } - ################## Return the model chosen - # Results (RSE, model coefficient, residuals, R?) - out <- list( input = list(H = Hdata$H, D = Hdata$D), model = output$mod, @@ -264,45 +233,56 @@ modelHD <- function(D, H, method = NULL, useWeight = FALSE, drawGraph = FALSE, p out$RSElog <- output$RSElog } + ####### if drawGraph is true + if (drawGraph) { + fitPlot <- drawPlotBegin(method, plot) + + geom_smooth( + data = data.frame(x = D_Plot$D, y = output$Hpredict_plot), mapping = aes(x,y), + formula = y~x, method = "loess" + ) + print(fitPlot) + out$fitPlot <- fitPlot + } + return(out) - } else { + } else { # If no method was supplied # Compare Models ---------------------------------------------------------- - if (is.null(plot)) { - drawPlotBegin(plotId = plot) - } - color <- c("blue", "green", "orange", "purple") - result <- rbindlist(lapply(1:length(methods), function(i) { + output_list <- list() + df_plot_list <- list() + for(i in 1:length(methods)) { method <- methods[i] - out <- modSelect(Hdata, method, useGraph = is.null(plot)) + out <- modSelect(Hdata, method, useGraph = drawGraph) - output <- list( - method = method, color = color[i], + output_list[[i]] <- list( + method = method, RSE = out$RSE, # Residual standard error RSElog = out$RSElog, Average_bias = out$Average_bias ) - if (is.null(plot)) { - lines(D_Plot$D, out$Hpredict_plot, lwd = 2, col = color[i], lty = i) - } - if (!is.null(plot)) { - output[["color"]] <- NULL + if (drawGraph) { + df_plot_list[[i]] <- data.frame(x = D_Plot$D, y = out$Hpredict_plot, method=methods[i]) } + } + result <- rbindlist(output_list,fill=T) - return(output) - }), fill = TRUE) - - if (is.null(plot)) { - legend("bottomright", methods, - lty = 1:5, lwd = 2, cex = 1, - col = color - ) + if (drawGraph) { + fitPlot <- drawPlotBegin(plotId = plot) + df_plot <- rbindlist(df_plot_list) + fitPlot <- fitPlot + + geom_smooth( + data = df_plot, mapping = aes(x,y,colour = method,linetype = method), + formula = y~x, method = "loess") + + theme(legend.position = "bottom") + print(fitPlot) } message("To build a HD model you must use the parameter 'method' in this function") + + return(data.frame(result)) } } diff --git a/vignettes/BIOMASS.Rmd b/vignettes/BIOMASS.Rmd index 817e887..356356a 100644 --- a/vignettes/BIOMASS.Rmd +++ b/vignettes/BIOMASS.Rmd @@ -53,7 +53,7 @@ KarnatakaForestsub <- droplevels(KarnatakaForest[selecPlot, ]) **First, check for any typo in the taxonomy** ```{r eval=test, cache=CACHE} -Taxo <- correctTaxo(genus = KarnatakaForestsub$genus, species = KarnatakaForestsub$species, useCache = FALSE, verbose = FALSE) +Taxo <- correctTaxo(genus = KarnatakaForestsub$genus, species = KarnatakaForestsub$species, useCache = TRUE, verbose = FALSE) KarnatakaForestsub$genusCorr <- Taxo$genusCorrected KarnatakaForestsub$speciesCorr <- Taxo$speciesCorrected ``` @@ -128,14 +128,16 @@ HDmodel <- modelHD( ) ``` -**Compute models specific to given stands** +**Compute models specific to given stands (given plots)** ```{r, cache=CACHE} HDmodelPerPlot <- modelHD( D = NouraguesHD$D, H = NouraguesHD$H, method = "weibull", - useWeight = TRUE, plot = NouraguesHD$plotId + useWeight = TRUE, plot = NouraguesHD$plotId, drawGraph = T ) +# To get model parameters by stand : ResHD <- t(sapply(HDmodelPerPlot, function(x) c(coef(x$model), RSE = x$RSE))) kable(ResHD, row.names = TRUE, digits = 3) + ``` From 6a471075cd4a53ff3c58db42df4a5871fb0037a9 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Mon, 21 Oct 2024 11:56:14 +0200 Subject: [PATCH 12/25] - Added element_blank as an import from ggplot2 - Improved modelHD documentation - Update NAMESPACE --- NAMESPACE | 7 +++++ R/checkPlotCoord.R | 4 +-- R/cutPlot.R | 2 +- R/modelHD.R | 26 +++++++++------ man/cutPlot.Rd | 2 +- man/modelHD.Rd | 40 ++++++++++++++---------- tests/testthat/test_00_modelHD.R | 14 ++++----- tests/testthat/test_01_correctCoordGPS.R | 12 +++---- 8 files changed, 63 insertions(+), 44 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index bd6aa13..dd36c22 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -42,16 +42,23 @@ importFrom(data.table,tstrsplit) importFrom(ggplot2,aes) importFrom(ggplot2,arrow) importFrom(ggplot2,coord_equal) +importFrom(ggplot2,element_blank) +importFrom(ggplot2,element_text) importFrom(ggplot2,geom_point) importFrom(ggplot2,geom_polygon) importFrom(ggplot2,geom_segment) +importFrom(ggplot2,geom_smooth) importFrom(ggplot2,geom_text) importFrom(ggplot2,ggplot) importFrom(ggplot2,ggtitle) +importFrom(ggplot2,labs) importFrom(ggplot2,scale_color_manual) importFrom(ggplot2,scale_shape_manual) +importFrom(ggplot2,scale_x_continuous) +importFrom(ggplot2,scale_y_continuous) importFrom(ggplot2,theme) importFrom(ggplot2,theme_minimal) +importFrom(ggplot2,unit) importFrom(grDevices,terrain.colors) importFrom(graphics,axis) importFrom(graphics,grid) diff --git a/R/checkPlotCoord.R b/R/checkPlotCoord.R index 2a02190..7d1394d 100644 --- a/R/checkPlotCoord.R +++ b/R/checkPlotCoord.R @@ -34,7 +34,7 @@ #' #' @importFrom data.table data.table := setnames #' @importFrom sf st_multipoint st_polygon st_sfc -#' @importFrom ggplot2 ggplot aes geom_point geom_segment geom_polygon geom_text scale_shape_manual scale_color_manual ggtitle theme_minimal theme coord_equal arrow unit +#' @importFrom ggplot2 ggplot aes geom_point geom_segment geom_polygon geom_text scale_shape_manual scale_color_manual ggtitle theme_minimal theme coord_equal arrow unit element_blank #' #' @author Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY #' @@ -78,7 +78,7 @@ #' checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPScorners, cornerID=NULL, maxDist = 15, rmOutliers = TRUE, drawPlot = TRUE, treeCoord = NULL) { - # parameters verification ------------------------------------------------- + # Checking arguments ------------------------------------------------- if (is.null(longlat) && is.null(projCoord)) { stop("Give at least one set of coordinates: longlat or projCoord") diff --git a/R/cutPlot.R b/R/cutPlot.R index c082f5d..73476d2 100644 --- a/R/cutPlot.R +++ b/R/cutPlot.R @@ -4,7 +4,7 @@ if (getRversion() >= "2.15.1") { )) } -#' Divides a plot in subplots +#' Divides one ore more plots into subplots #' #' This function divides a plot (or several plots) in subplots and returns the coordinates of the grid. #' This function uses a procrust analysis to fit the rectangle you gave to the plot you have. diff --git a/R/modelHD.R b/R/modelHD.R index 58c365c..21a00e4 100644 --- a/R/modelHD.R +++ b/R/modelHD.R @@ -44,7 +44,6 @@ #' If the parameter model is null, the function return a plot with all the methods for #' comparison, the function also returns a data.frame with: #' - `method`: The method that had been used to construct the plot -#' - `color`: The color of the curve in the plot #' - `RSE`: Residual Standard Error of the model #' - `RSElog`: Residual Standard Error of the log model (`NULL` if other model) #' - `Average_bias`: The average bias for the model @@ -61,29 +60,36 @@ #' # Load a data set #' data(NouraguesHD) #' -#' # To model the height from a dataset +#' # Fit H-D models for the Nouragues dataset #' \donttest{ #' HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, drawGraph = TRUE) #' } #' -#' # If the method needed is known -#' HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, method = "weibull", drawGraph = TRUE) -#' HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, method = "log1", drawGraph = TRUE) +#' # For a chosen model +#' HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, +#' method = "log2", drawGraph = TRUE) #' #' # Using weights #' HDmodel <- modelHD( -#' D = NouraguesHD$D, H = NouraguesHD$H, method = "weibull", useWeight = TRUE, -#' drawGraph = TRUE -#' ) +#' D = NouraguesHD$D, H = NouraguesHD$H, +#' method = "log2", useWeight = TRUE, +#' drawGraph = TRUE) +#' +#' # With multiple stands (plots) +#' HDmodel <- modelHD( +#' D = NouraguesHD$D, H = NouraguesHD$H, +#' method = "log2", useWeight = TRUE, +#' plot = NouraguesHD$plotId, drawGraph = TRUE) +#' #' @importFrom graphics legend lines par plot grid axis #' @importFrom stats SSmicmen lm median na.omit quantile rnorm sd predict coef #' @importFrom utils data #' @importFrom data.table data.table -#' @importFrom ggplot2 ggplot aes geom_point geom_smooth labs theme_minimal theme scale_x_continuous scale_y_continuous +#' @importFrom ggplot2 ggplot aes geom_point geom_smooth labs theme_minimal theme scale_x_continuous scale_y_continuous element_text modelHD <- function(D, H, method = NULL, useWeight = FALSE, drawGraph = FALSE, plot = NULL) { - # parameters verification ------------------------------------------------- + # Checking arguments ------------------------------------------------- # Check if there is enough data to compute an accurate model nbNonNA <- sum(!is.na(H)) diff --git a/man/cutPlot.Rd b/man/cutPlot.Rd index c6e6f8b..e514478 100644 --- a/man/cutPlot.Rd +++ b/man/cutPlot.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/cutPlot.R \name{cutPlot} \alias{cutPlot} -\title{Divides one ore more plots in subplots} +\title{Divides one ore more plots into subplots} \usage{ cutPlot(projCoord, plot, cornerNum, gridsize = 100, dimX = 200, dimY = 200) } diff --git a/man/modelHD.Rd b/man/modelHD.Rd index 46a59d7..d90fbac 100644 --- a/man/modelHD.Rd +++ b/man/modelHD.Rd @@ -7,11 +7,9 @@ modelHD(D, H, method = NULL, useWeight = FALSE, drawGraph = FALSE, plot = NULL) } \arguments{ -\item{D}{Vector with diameter measurements (in cm). NA values are accepted but a -minimum of 10 valid entries (i.e. having a corresponding height in H) is required.} +\item{D}{Vector with diameter measurements (in cm). NA values are accepted but a minimum of 10 valid entries (i.e. having a corresponding height in H) is required.} -\item{H}{Vector with total height measurements (in m). NA values are accepted but a minimum of 10 valid entries -(i.e. having a corresponding diameter in D) is required.} +\item{H}{Vector with total height measurements (in m). NA values are accepted but a minimum of 10 valid entries (i.e. having a corresponding diameter in D) is required.} \item{method}{Method used to fit the relationship. To be chosen between: @@ -36,9 +34,10 @@ volume, so that larger trees have a stronger influence during the construction o stand-specific HD models.} } \value{ -If plot is NULL or has a single value, a single list is returned. If there is more than one plot, +If \code{plot} is NULL or has a single value, a single list is returned. If there is more than one plot, multiple embedded lists are returned with plots as the list names. -Returns a list if the parameter model is not null: + +If \code{model} is not null (model comparison), returns a list : \itemize{ \item \code{input}: list of the data used to construct the model (list(H, D)) \item \code{model}: outputs of the model (same outputs as given by \code{\link[stats:lm]{stats::lm()}}, \code{\link[stats:nls]{stats::nls()}}) @@ -50,13 +49,13 @@ Returns a list if the parameter model is not null: \item \code{formula}: Formula of the model \item \code{method}: Name of the method used to construct the model \item \code{predicted}: Predicted height values +\item \code{fitPlot}: a ggplot object containing the model fitting plot } -If the parameter model is null, the function return a graph with all the methods for +If the parameter model is null, the function return a plot with all the methods for comparison, the function also returns a data.frame with: \itemize{ -\item \code{method}: The method that had been used to construct the graph -\item \code{color}: The color of the curve in the graph +\item \code{method}: The method that had been used to construct the plot \item \code{RSE}: Residual Standard Error of the model \item \code{RSElog}: Residual Standard Error of the log model (\code{NULL} if other model) \item \code{Average_bias}: The average bias for the model @@ -74,24 +73,31 @@ where RSE is the Residual Standard Error). # Load a data set data(NouraguesHD) -# To model the height from a dataset +# Fit H-D models for the Nouragues dataset \donttest{ HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, drawGraph = TRUE) } -# If the method needed is known -HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, method = "weibull", drawGraph = TRUE) -HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, method = "log1", drawGraph = TRUE) +# For a chosen model +HDmodel <- modelHD(D = NouraguesHD$D, H = NouraguesHD$H, +method = "log2", drawGraph = TRUE) # Using weights HDmodel <- modelHD( - D = NouraguesHD$D, H = NouraguesHD$H, method = "weibull", useWeight = TRUE, - drawGraph = TRUE -) + D = NouraguesHD$D, H = NouraguesHD$H, + method = "log2", useWeight = TRUE, + drawGraph = TRUE) + +# With multiple stands (plots) +HDmodel <- modelHD( + D = NouraguesHD$D, H = NouraguesHD$H, + method = "log2", useWeight = TRUE, + plot = NouraguesHD$plotId, drawGraph = TRUE) + } \seealso{ \code{\link[=retrieveH]{retrieveH()}} } \author{ -Maxime REJOU-MECHAIN, Arthur PERE, Ariane TANGUY +Maxime REJOU-MECHAIN, Arthur PERE, Ariane TANGUY, Arthur Bailly } diff --git a/tests/testthat/test_00_modelHD.R b/tests/testthat/test_00_modelHD.R index 9ac0872..dd2d5e1 100644 --- a/tests/testthat/test_00_modelHD.R +++ b/tests/testthat/test_00_modelHD.R @@ -1,4 +1,4 @@ -require(data.table) +#require(data.table) D <- NouraguesHD$D H <- NouraguesHD$H @@ -58,13 +58,13 @@ test_that("Without parameters", { Res <- expect_message(modelHD(D, H, useWeight = TRUE), "build a HD model") expect_is(Res, "data.frame") - expect_equal(ncol(Res), 5) + expect_equal(ncol(Res),4) - res <- "method color RSE RSElog Average_bias -log1 blue 4.305060 0.2231136 0.004227454 -log2 green 4.222718 0.2215495 0.003121671 -weibull orange 4.307951 NA 0.002823978 -michaelis purple 4.294488 NA 0.014564152 + res <- "method RSE RSElog Average_bias +log1 4.305060 0.2231136 0.004227454 +log2 4.222718 0.2215495 0.003121671 +weibull 4.307951 NA 0.002823978 +michaelis 4.294488 NA 0.014564152 " diff --git a/tests/testthat/test_01_correctCoordGPS.R b/tests/testthat/test_01_correctCoordGPS.R index 3971462..e10cb97 100644 --- a/tests/testthat/test_01_correctCoordGPS.R +++ b/tests/testthat/test_01_correctCoordGPS.R @@ -71,26 +71,26 @@ test_that("correct coord GPS in UTM", { }) test_that("correct coord GPS error", { - expect_error(correctCoordGPS(), "Give at least one set of coordinates") + expect_error(suppressWarnings(correctCoordGPS()), "Give at least one set of coordinates") expect_error( - correctCoordGPS(projCoord = projCoord, coordRel = coordRel, rangeX = 52, rangeY = 53,rmOutliers = FALSE), + suppressWarnings(correctCoordGPS(projCoord = projCoord, coordRel = coordRel, rangeX = 52, rangeY = 53,rmOutliers = FALSE)), "must be of length 2" ) expect_error( - correctCoordGPS(projCoord = projCoord, coordRel = coordRel, rangeX = c(0, 100), rangeY = c(0, 100), maxDist = c(15, 0),rmOutliers = FALSE), + suppressWarnings(correctCoordGPS(projCoord = projCoord, coordRel = coordRel, rangeX = c(0, 100), rangeY = c(0, 100), maxDist = c(15, 0),rmOutliers = FALSE)), "Your argument maxDist must be of length 1" ) expect_error( - correctCoordGPS(projCoord = projCoord, coordRel = coordRel, rangeX = c(0, 40), rangeY = c(0, 40),rmOutliers = FALSE), + suppressWarnings(correctCoordGPS(projCoord = projCoord, coordRel = coordRel, rangeX = c(0, 40), rangeY = c(0, 40),rmOutliers = FALSE)), "coordRel must be inside the X and Y ranges" ) expect_error( - correctCoordGPS(projCoord = projCoord, coordRel = rbind(coordRel, c(40, 40)), rangeX = c(0, 100), rangeY = c(0, 100),rmOutliers = FALSE), + suppressWarnings(correctCoordGPS(projCoord = projCoord, coordRel = rbind(coordRel, c(40, 40)), rangeX = c(0, 100), rangeY = c(0, 100),rmOutliers = FALSE)), "same dimension" ) skip_if_not_installed("proj4") expect_error( - correctCoordGPS(longlat = c(15, 12), projCoord = projCoord,rmOutliers = FALSE), + suppressWarnings(correctCoordGPS(longlat = c(15, 12), projCoord = projCoord,rmOutliers = FALSE)), "Give only one set of coordinates" ) }) From 5b92f521bdd6317f0908e8957e009b473b56cbf1 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Mon, 21 Oct 2024 14:07:05 +0200 Subject: [PATCH 13/25] Replaced setDT(copy()) by data.table() --- R/getTaxonomy.R | 4 ++-- R/getWoodDensity.R | 10 +++++----- data-raw/apgfamily.R | 2 +- data-raw/genusFamily.R | 2 +- data-raw/sd_10.R | 2 +- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/R/getTaxonomy.R b/R/getTaxonomy.R index d467f28..d9429f0 100644 --- a/R/getTaxonomy.R +++ b/R/getTaxonomy.R @@ -32,7 +32,7 @@ getTaxonomy <- function(genus, findOrder = FALSE) { ################## 1. Retrieve the Family # Load taxonomical data (sourced from Angiosperm Phylogeny Website, http://www.mobot.org/MOBOT/research/APweb/) - genusFamily <- setDT(copy(BIOMASS::genusFamily)) + genusFamily <- data.table(BIOMASS::genusFamily) setkey(genusFamily, genus) # Create ids @@ -48,7 +48,7 @@ getTaxonomy <- function(genus, findOrder = FALSE) { ################## 2. Retrieve the Order if (findOrder == TRUE) { - apgFamilies <- setDT(copy(BIOMASS::apgFamilies)) + apgFamilies <- data.table(BIOMASS::apgFamilies) genusFam <- merge(genusFam, apgFamilies, by.x = "family", by.y = "famAPG", all.x = TRUE) genusFam <- genusFam[, .(id, inputGenus, family, order)] diff --git a/R/getWoodDensity.R b/R/getWoodDensity.R index efd72db..a8e11d6 100644 --- a/R/getWoodDensity.R +++ b/R/getWoodDensity.R @@ -111,7 +111,7 @@ if (getRversion() >= "2.15.1") { #' } #' @seealso [wdData], [sd_10] #' @keywords Wood density -#' @importFrom data.table data.table := setDF setDT setkey copy chmatch %chin% +#' @importFrom data.table data.table := setDF setDT setkey chmatch %chin% #' getWoodDensity <- function(genus, species, stand = NULL, family = NULL, region = "World", addWoodDensityData = NULL, verbose = TRUE) { @@ -148,11 +148,11 @@ getWoodDensity <- function(genus, species, stand = NULL, family = NULL, region = # Data processing --------------------------------------------------------- # Load global wood density database downloaded from http://datadryad.org/handle/10255/dryad.235 - wdData <- setDT(copy(BIOMASS::wdData)) + wdData <- data.table(BIOMASS::wdData) # Load the mean standard deviation observed at the species, Genus or Family level # in the Dryad dataset when at least 10 individuals are considered - sd_10 <- setDT(copy(BIOMASS::sd_10)) + sd_10 <- data.table(BIOMASS::sd_10) sd_tot <- sd(wdData$wd) Region <- tolower(region) @@ -175,7 +175,7 @@ getWoodDensity <- function(genus, species, stand = NULL, family = NULL, region = if (!is.null(addWoodDensityData)) { setDT(addWoodDensityData) if (!("family" %in% names(addWoodDensityData))) { - genusFamily <- setDT(copy(BIOMASS::genusFamily)) + genusFamily <- data.table(BIOMASS::genusFamily) addWoodDensityData[genusFamily, on = "genus", family := i.family] } addWoodDensityData <- addWoodDensityData[!is.na(wd), ] @@ -195,7 +195,7 @@ getWoodDensity <- function(genus, species, stand = NULL, family = NULL, region = inputData[, family := as.character(family)] } else { if (!exists("genusFamily", inherits = FALSE)) { - genusFamily <- setDT(copy(BIOMASS::genusFamily)) + genusFamily <- data.table(BIOMASS::genusFamily) } inputData[genusFamily, on = "genus", family := i.family] } diff --git a/data-raw/apgfamily.R b/data-raw/apgfamily.R index 9015930..f5eb17d 100644 --- a/data-raw/apgfamily.R +++ b/data-raw/apgfamily.R @@ -22,7 +22,7 @@ setDF(apgFamilies) usethis::use_data(apgFamilies, compress = "xz", overwrite = FALSE) -apgFamiliesOrg <- setDT(copy(BIOMASS::apgFamilies)) +apgFamiliesOrg <- data.table(BIOMASS::apgFamilies) apgFamiliesOrg <- apgFamiliesOrg[base::order(order)] apgFamilies <- apgFamilies[base::order(order)] diff --git a/data-raw/genusFamily.R b/data-raw/genusFamily.R index 6fa32c7..3e0f9d5 100644 --- a/data-raw/genusFamily.R +++ b/data-raw/genusFamily.R @@ -261,7 +261,7 @@ usethis::use_data(genusFamily, compress = "xz") genusFamily <- genusFamily[base::order(family)] -genusFamilyOrg <- setDT(copy(BIOMASS::genusFamily)) +genusFamilyOrg <- data.table(BIOMASS::genusFamily) genusFamilyOrg <- genusFamilyOrg[base::order(family)] result <- merge(genusFamilyOrg, genusFamily, all = T, by = "genus") diff --git a/data-raw/sd_10.R b/data-raw/sd_10.R index 78783c4..565767a 100644 --- a/data-raw/sd_10.R +++ b/data-raw/sd_10.R @@ -7,7 +7,7 @@ require(data.table) # ) # data_wd[, c("genus", "species") := tstrsplit(Binomial, " ", keep = 1:2)] -data_wd <- setDT(copy(BIOMASS::wdData)) +data_wd <- data.table(BIOMASS::wdData) species <- data_wd[, .(family, genus, wd, .N), by = species][N >= 10][, .(sd = sd(wd)), by = species][, mean(sd)] From e03d96c97ad5c902408b465848ae68eac67d3193 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 23 Oct 2024 18:00:53 +0200 Subject: [PATCH 14/25] - Modified bilinear interpolation to generalized bilinear interpolation that is able to convert relative coordinates to projected coordinates even if the plot is a rectangle, and even if the relative origin is not (0;0) - bilinearInterpolation doesn't take into account cornerNumber anymore, making it more general - Modified corner number assignment to corner row arrangement in checkPlotCoord() - Added bilinearInterpolation tests --- NAMESPACE | 3 +- R/bilinearInterpolation.R | 96 +++++++++++++------ R/checkPlotCoord.R | 60 ++++++------ man/BIOMASS-package.Rd | 1 + man/bilinearInterpolation.Rd | 33 +++++-- man/cutPlot.Rd | 2 +- tests/testthat/test-bilinearInterpolation.R | 71 ++++++++++++++ tests/testthat/test_01_small_function.R | 18 ---- ...SpatialRepresentationAndManagementPlot.Rmd | 2 +- 9 files changed, 195 insertions(+), 91 deletions(-) create mode 100644 tests/testthat/test-bilinearInterpolation.R diff --git a/NAMESPACE b/NAMESPACE index dd36c22..9e41af3 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ export(AGBmonteCarlo) export(attributeTree) export(attributeTreeCoord) +export(bilinearInterpolation) export(cacheManager) export(cachePath) export(checkPlotCoord) @@ -28,11 +29,11 @@ importFrom(data.table,":=") importFrom(data.table,as.data.table) importFrom(data.table,between) importFrom(data.table,chmatch) -importFrom(data.table,copy) importFrom(data.table,data.table) importFrom(data.table,first) importFrom(data.table,fread) importFrom(data.table,fwrite) +importFrom(data.table,is.data.table) importFrom(data.table,rbindlist) importFrom(data.table,setDF) importFrom(data.table,setDT) diff --git a/R/bilinearInterpolation.R b/R/bilinearInterpolation.R index 4eba3fa..df3ce87 100644 --- a/R/bilinearInterpolation.R +++ b/R/bilinearInterpolation.R @@ -1,42 +1,82 @@ -#' Bilinear interpolation of X and Y absolute coordinates +#' Generalized bilinear interpolation of coordinates #' -#' Apply a bilinear interpolation to convert relative coordinates into absolute coordinates, using the 4 corners absolute coordinates of the plot. +#' @description Apply a generalized bilinear interpolation to convert any coordinates from one original coordinate system to another, using the plot's 4 corner coordinates of both system. #' -#' @param relCoord a matrix : the relative coordinates to transform -#' @param cornerCoord a data.table : the absolute corners coordinates and its numbering following a clockwise direction (names resp. "X","Y","cornerNum") -#' @param dimX a vector indicating the size of the plot on the X axis, in meters and in the relative coordinates system -#' @param dimY a vector indicating the size of the plot on the Y axis, in meters and in the relative coordinates system +#' @details +#' The plot represented by the 4 coordinates in fromCornerCoord must have 4 right angles, i.e. a rectangular (or square) plot. +#' +#' @references C. -C. Wei and C. -H. Chen, "Generalized Bilinear Interpolation of Motion Vectors for Quad-Tree Mesh," 2008 International Conference on Intelligent Information Hiding and Multimedia Signal Processing, Harbin, China, 2008, pp. 635-638, doi: 10.1109/IIH-MSP.2008.283. +#' +#' +#' @param coord a matrix or data.frame : coordinates to be transformed, with X and Y corresponding to the first two columns +#' @param fromCornerCoord a matrix or data.frame : corner coordinates of the plot in the original coordinate system, with X and Y corresponding to the first two columns +#' @param toCornerCoord a matrix or data.frame : corner coordinates of the plot in the coordinate system to be projected, with the same line order as fromCornerCoord and , with X and Y corresponding to the first two columns +#' @param orderedCorner a logical, indicating if fromCornerCoord and toCornerCoord rows are sorted in correct order (clockwise or counter-clockwise) #' -#' @return a data.table containing the absolute coordinates of desired points +#' @return a data.frame containing the converted coordinates #' -#' @importFrom data.table data.table := +#' @importFrom data.table is.data.table #' -#' @keywords Internal bilinear interpolation +#' @export +#' +#' @keywords generalized bilinear interpolation #' #' @author Arthur Bailly #' -bilinearInterpolation = function(relCoord, cornerCoord, dimX, dimY) { - - # See https://en.wikipedia.org/wiki/Bilinear_interpolation#Alternative_matrix_form for details +#' @examples +#' fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) +#' rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) +#' toCornerCoord <- as.matrix(fromCornerCoord) %*% rot_mat +#' toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") +#' coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) +#' projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) +#' # plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(fromCornerCoord) ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") + + +bilinearInterpolation = function(coord, fromCornerCoord, toCornerCoord, orderedCorner = F) { - setnames(cornerCoord,c("X","Y","cornerNum")) + # Parameters verification + if(nrow(fromCornerCoord)!=4 | nrow(toCornerCoord)!=4 | nrow(fromCornerCoord)!=nrow(fromCornerCoord)) { + stop("fromCornerCoord and toCornerCoord must have 4 rows representing the 4 corners of the plot") + } + if(!(is.data.frame(coord) | is.matrix(coord) | is.data.table(coord))){ + stop("tree coordinates must be a data.frame, a matrix or a data.table") + } + if(is.data.table(coord)) coord <- data.frame(coord) + if(is.data.table(fromCornerCoord) | is.data.table(toCornerCoord)) { + fromCornerCoord <- data.frame(fromCornerCoord) + toCornerCoord <- data.frame(toCornerCoord) + } - denominator <- (dimX-0)*(dimY-0) + # Sorting rows if necessary + centroid <- colMeans(fromCornerCoord) + if(!orderedCorner) { + # Sort fromCornerCoord and toCornerCoord rows in a counter-clockwise direction + angles <- atan2(fromCornerCoord[, 2] - centroid[2], fromCornerCoord[, 1] - centroid[1]) + fromCornerCoord <- fromCornerCoord[order(angles), ] + toCornerCoord <- toCornerCoord[order(angles), ] + } + + # Verification of a rectangular plot for fromCornerCoord + if(!all(abs(dist(rbind(fromCornerCoord,centroid))[c(4,7,9,10)] - mean(dist(rbind(fromCornerCoord,centroid))[c(4,7,9,10)]))<0.1)) { + stop("The plot in the relative coordinate system is not a rectangle (or a square). You may consider using trustGPScorners = F") + } - multMat <- matrix(c(dimX*dimY , -dimX*0 , -0*dimY , 0*0, - -dimY,0,dimY,-0, - -dimX,dimX,0,-0, - 1,-1,-1,1) , ncol = 4) + x_A <- fromCornerCoord[1,1] ; x_B <- fromCornerCoord[2,1] ; x_C <- fromCornerCoord[3,1] ; x_D <- fromCornerCoord[4,1] + y_A <- fromCornerCoord[1,2] ; y_B <- fromCornerCoord[2,2] ; y_C <- fromCornerCoord[3,2] ; y_D <- fromCornerCoord[4,2] + u_A <- toCornerCoord[1,1] ; u_B <- toCornerCoord[2,1]; u_C <- toCornerCoord[3,1]; u_D <- toCornerCoord[4,1] + v_A <- toCornerCoord[1,2] ; v_B <- toCornerCoord[2,2] ; v_C <- toCornerCoord[3,2] ; v_D <- toCornerCoord[4,2] - # Apply the bilinear interpolation formula to a single coordinate point (x;y) - bilinearInterpolationFunction = function(x,y) { - apply(cornerCoord[order(cornerNum),.(X,Y)][c(1,2,4,3),] , 2 , function(i) 1/denominator * i %*% multMat %*% matrix(c(1,x,y,x*y)) ) + apply_bilinear_interpolation <- function(x,y) { + rate_A <- (1-(x-x_A)/(x_C-x_A)) * (1-(y-y_A)/(y_C-y_A)) + rate_B <- (1-(x-x_B)/(x_D-x_B)) * (1-(y-y_B)/(y_D-y_B)) + rate_C <- (1-(x-x_C)/(x_A-x_C)) * (1-(y-y_C)/(y_A-y_C)) + rate_D <- (1-(x-x_D)/(x_B-x_D)) * (1-(y-y_D)/(y_B-y_D)) + return(data.frame( + X = rate_A*u_A + rate_B*u_B + rate_C*u_C + rate_D*u_D, + Y = rate_A*v_A + rate_B*v_B + rate_C*v_C + rate_D*v_D) + ) } - # Apply the bilinearInterpolationFunction to all points of relCoord - absCoord = apply(relCoord , 1 , function(dat) { - bilinearInterpolationFunction(x=dat[1] , y=dat[2]) - }) - - return(data.table(XAbs=absCoord[1,] , YAbs=absCoord[2,])) -} \ No newline at end of file + return(apply_bilinear_interpolation(x=coord[,1],y=coord[,2])) +} diff --git a/R/checkPlotCoord.R b/R/checkPlotCoord.R index 7d1394d..76a5227 100644 --- a/R/checkPlotCoord.R +++ b/R/checkPlotCoord.R @@ -32,7 +32,7 @@ #' #' @export #' -#' @importFrom data.table data.table := setnames +#' @importFrom data.table data.table := setnames %between% #' @importFrom sf st_multipoint st_polygon st_sfc #' @importFrom ggplot2 ggplot aes geom_point geom_segment geom_polygon geom_text scale_shape_manual scale_color_manual ggtitle theme_minimal theme coord_equal arrow unit element_blank #' @@ -144,7 +144,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc if(trustGPScorners == TRUE) { - if(nrow(projCoord)!= 4) { # if multiple measures of each corner, then do the mean of coordinates and look for outliers + if(nrow(projCoord)!= 4) { # if multiple measures of each corner, then do the mean of coordinates and search for outliers cornerCoord <- data.table(cbind(projCoord[,1:2], relCoord[,1:2],cornerID=cornerID)) setnames(cornerCoord, c("X","Y","Xrel","Yrel","cornerID")) @@ -173,14 +173,13 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc } } - cornerCoord <- cornerCoord[ , c("Xmean","Ymean","Xrel","Yrel","cornerID")] + cornerCoord <- data.frame(cornerCoord[ , c("Xmean","Ymean","Xrel","Yrel","cornerID")]) cornerCoord <- unique(cornerCoord) - setnames(cornerCoord, colnames(cornerCoord), c("X", "Y","Xrel","Yrel","cornerID")) + colnames(cornerCoord) <- c("X", "Y","Xrel","Yrel","cornerID") } else { # if exactly 1 measures for each corner - - cornerCoord <- data.table(cbind(projCoord[,1:2], relCoord[,1:2])) - setnames(cornerCoord, c("X","Y","Xrel","Yrel")) + cornerCoord <- data.frame(cbind(projCoord[,1:2], relCoord[,1:2])) + colnames(cornerCoord) <- c("X","Y","Xrel","Yrel") outliers <- data.frame() if(!is.null(cornerID)) cornerCoord$cornerID = cornerID } @@ -233,28 +232,16 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc cornerProjCoord <- as.matrix(cornerRelCoord[,1:2]) %*% procrustRes$rotation cornerProjCoord <- sweep(cornerProjCoord, 2, procrustRes$translation, FUN = "+") - cornerCoord <- data.table(cbind(cornerProjCoord[,1:2], cornerRelCoord[,1:2])) - setnames(cornerCoord, c("X","Y","Xrel","Yrel")) + cornerCoord <- data.frame(cbind(cornerProjCoord[,1:2], cornerRelCoord[,1:2])) + colnames(cornerCoord) <- c("X","Y","Xrel","Yrel") if(!is.null(cornerID)) cornerCoord$cornerID <- unique(cornerID) } # End trustGPScorners = "FALSE" - # Assign corner numbers in a clockwise direction ----------------------------- - m1 <- cornerCoord[ rank(Xrel) <= 2, ] - tmp1 <- m1[rank(Yrel),] - m2 <- cornerCoord[ rank(Xrel) > 2, ] - tmp2 <- m2[rank(Yrel),] - cornerCoord <- rbind(tmp1,tmp2) - cornerCoord$cornerNum <- c(1,2,4,3) - - # if(! is.null(originalCorner)) { - # # Shift the corner numbers to have corner 1 on the origin - # shift <- 5 - which(cornerCoord$cornerID == originalCorner) - # cornerCoord$cornerNum <- (cornerCoord$cornerNum + shift) %% 4 - # cornerCoord$cornerNum[cornerCoord$cornerNum == 0] <- 4 - # } - - cornerCoord <- cornerCoord[ order(cornerNum), ] + # Sort cornerCoord rows in a counter-clockwise direction --------------------- + centroid <- colMeans(cornerCoord[,c("Xrel","Yrel")]) + angles <- base::atan2(cornerCoord[["Yrel"]] - centroid[2], cornerCoord[["Xrel"]] - centroid[1]) + cornerCoord <- cornerCoord[order(angles), ] # Create a polygon ----------------------------------------------------------- cornerPolygon <- st_multipoint(as.matrix(rbind(cornerCoord[,c("X","Y")], cornerCoord[1,c("X","Y")]))) @@ -263,16 +250,23 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc # Calculate projected tree coordinates from relative tree coordinates -------- if(!is.null(treeCoord)) { + if(any(! treeCoord[,1] %between% range(cornerCoord$Xrel)) | any(! treeCoord[,2] %between% range(cornerCoord$Yrel))) { + warning("Be careful, one or more trees are not inside the plot defined by relCoord") + } if(trustGPScorners) { - treeProjCoord <- bilinearInterpolation(relCoord = treeCoord[,1:2], - cornerCoord = cornerCoord[,c("X","Y","cornerNum")], - dimX = diff(range(cornerCoord$Xrel)), - dimY = diff(range(cornerCoord$Yrel)) - ) + treeProjCoord <- bilinearInterpolation(coord = treeCoord[,1:2], + fromCornerCoord = cornerCoord[,c("Xrel","Yrel")], + toCornerCoord = cornerCoord[,c("X","Y")], + orderedCorner = T) + # treeProjCoord <- bilinearInterpolation(relCoord = treeCoord[,1:2], + # cornerCoord = cornerCoord[,c("X","Y","cornerNum")], + # dimX = diff(range(cornerCoord$Xrel)), + # dimY = diff(range(cornerCoord$Yrel)) + # ) colnames(treeProjCoord) <- c("X","Y") } else { treeProjCoord <- as.matrix(treeCoord[,1:2]) %*% procrustRes$rotation - treeProjCoord <- data.table(sweep(treeProjCoord, 2, procrustRes$translation, FUN = "+")) + treeProjCoord <- data.frame(sweep(treeProjCoord, 2, procrustRes$translation, FUN = "+")) colnames(treeProjCoord) <- c("X","Y") } } @@ -322,7 +316,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc } output <- list( - cornerCoord = data.frame(cornerCoord), + cornerCoord = cornerCoord, polygon = cornerPolygon, outliers = outliers ) if (drawPlot) { @@ -332,7 +326,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc output$codeUTM <- codeUTM } if (!is.null(treeCoord)) { - output$treeProjCoord <- data.frame(treeProjCoord[,c("X","Y")]) + output$treeProjCoord <- treeProjCoord[,c("X","Y")] } return(output) } diff --git a/man/BIOMASS-package.Rd b/man/BIOMASS-package.Rd index 7160919..a55af88 100644 --- a/man/BIOMASS-package.Rd +++ b/man/BIOMASS-package.Rd @@ -169,6 +169,7 @@ Authors: \item Ariane Tanguy \item Camille Piponiot \item Bruno Hérault \email{bruno.herault@cirad.fr} + \item Arthur Bailly \email{arthur.bailly@ird.fr} } Other contributors: diff --git a/man/bilinearInterpolation.Rd b/man/bilinearInterpolation.Rd index 975c70a..2115a17 100644 --- a/man/bilinearInterpolation.Rd +++ b/man/bilinearInterpolation.Rd @@ -2,28 +2,43 @@ % Please edit documentation in R/bilinearInterpolation.R \name{bilinearInterpolation} \alias{bilinearInterpolation} -\title{Bilinear interpolation of X and Y absolute coordinates} +\title{Generalized bilinear interpolation of coordinates} \usage{ -bilinearInterpolation(relCoord, cornerCoord, dimX, dimY) +bilinearInterpolation(coord, fromCornerCoord, toCornerCoord, orderedCorner = F) } \arguments{ -\item{relCoord}{a matrix : the relative coordinates to transform} +\item{coord}{a matrix or data.frame : coordinates to be transformed, with X and Y corresponding to the first two columns} -\item{cornerCoord}{a data.table : the absolute corners coordinates and its numbering following a clockwise direction (names resp. "X","Y","cornerNum")} +\item{fromCornerCoord}{a matrix or data.frame : corner coordinates of the plot in the original coordinate system, with X and Y corresponding to the first two columns} -\item{dimX}{a vector indicating the size of the plot on the X axis, in meters and in the relative coordinates system} +\item{toCornerCoord}{a matrix or data.frame : corner coordinates of the plot in the coordinate system to be projected, with the same line order as fromCornerCoord and , with X and Y corresponding to the first two columns} -\item{dimY}{a vector indicating the size of the plot on the Y axis, in meters and in the relative coordinates system} +\item{orderedCorner}{a logical, indicating if fromCornerCoord and toCornerCoord rows are sorted in correct order (clockwise or counter-clockwise)} } \value{ -a data.table containing the absolute coordinates of desired points +a data.frame containing the converted coordinates } \description{ -Apply a bilinear interpolation to convert relative coordinates into absolute coordinates, using the 4 corners absolute coordinates of the plot. +Apply a generalized bilinear interpolation to convert any coordinates from one original coordinate system to another, using the plot's 4 corner coordinates of both system. +} +\details{ +The plot represented by the 4 coordinates in fromCornerCoord must have 4 right angles, i.e. a rectangular (or square) plot. +} +\examples{ +fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) +rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) +toCornerCoord <- as.matrix(fromCornerCoord) \%*\% rot_mat +toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") +coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) +projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) +# plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(fromCornerCoord) ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") +} +\references{ +C. -C. Wei and C. -H. Chen, "Generalized Bilinear Interpolation of Motion Vectors for Quad-Tree Mesh," 2008 International Conference on Intelligent Information Hiding and Multimedia Signal Processing, Harbin, China, 2008, pp. 635-638, doi: 10.1109/IIH-MSP.2008.283. } \author{ Arthur Bailly } -\keyword{Internal} \keyword{bilinear} +\keyword{generalized} \keyword{interpolation} diff --git a/man/cutPlot.Rd b/man/cutPlot.Rd index e514478..5eb8548 100644 --- a/man/cutPlot.Rd +++ b/man/cutPlot.Rd @@ -32,7 +32,7 @@ Returns a data-frame containing as many rows as there are corners corresponding } \description{ This function divides a plot (or several plots) in subplots and returns the coordinates of the grid. -This function uses a procrust analysis to fit the rectangle you gave to the plot you have. +These coordinates are calculated by a bilinear interpolation with the projected corner coordinates as references. } \examples{ diff --git a/tests/testthat/test-bilinearInterpolation.R b/tests/testthat/test-bilinearInterpolation.R new file mode 100644 index 0000000..37bf6da --- /dev/null +++ b/tests/testthat/test-bilinearInterpolation.R @@ -0,0 +1,71 @@ +context("Bilinear interpolation") + + +test_that("bilinearInterpolation error", { + fromCornerCoord <- data.frame(x=1:5,y=1:5) + toCornerCoord <- data.frame(x=1:4,y=1:4) + expect_error(bilinearInterpolation(coord = NULL, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) , "fromCornerCoord and toCornerCoord must have 4 rows representing the 4 corners of the plot") + expect_error(bilinearInterpolation(coord = NULL, fromCornerCoord = toCornerCoord, toCornerCoord = fromCornerCoord) , "fromCornerCoord and toCornerCoord must have 4 rows representing the 4 corners of the plot") + fromCornerCoord <- data.frame(x=c(0,0,10,10), y=c(0,10,10,0)) + toCornerCoord <- fromCornerCoord+50 + coord <- c(5,5) + expect_error(bilinearInterpolation(coord = coord, fromCornerCoord = toCornerCoord, toCornerCoord = fromCornerCoord) , "tree coordinates must be a data.frame, a matrix or a data.table") + coord <- data.frame(x=5,y=5) + fromCornerCoord[3,1] <- 11 + expect_error(bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) , "You may consider using trustGPScorners = F") +}) + +test_that("bilinearInterpolation function", { + + fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) + rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) + toCornerCoord <- as.matrix(fromCornerCoord) %*% rot_mat + toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") + #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) + coord <- data.frame(x=c(25,60),y=c(25,40)) + projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) + #plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(fromCornerCoord) ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") + + expect_is(projCoord,"data.frame") + expect_equal(nrow(projCoord), nrow(coord)) + expect_equal(names(projCoord), c("X","Y")) + + expect_equal(projCoord, data.frame(X=c(59.15064,81.96152),Y=c(134.1506,164.6410)) , tolerance = 1e-6) + + + # Corner deformation + set.seed(52) + toCornerCoord <- toCornerCoord + runif(1,-10,10) + # toCornerCoord <- fromCornerCoord + # toCornerCoord[2,] = c(103,-5) + # toCornerCoord[3,] = c(10,52) + # toCornerCoord[4,] = c(95,48) + # toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") + projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) + expect_equal(projCoord, data.frame(X=c(52.39103,75.20192),Y=c(127.3910,157.8814)) , tolerance = 1e-6) + #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) + #plot(fromCornerCoord, xlim=c(-10,150),ylim=c(-5,200)) ; points(coord , col="blue") ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") + + + # If the origin is not (0;0) + fromCornerCoord <- expand.grid(X = c(50, 150), Y = c(50, 100)) + toCornerCoord <- sweep(fromCornerCoord, 2, c(60,60), FUN = "+") + coord <- data.frame(x=c(25,60)+50,y=c(25,40)+50) + #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) + 50 + projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) + #plot(fromCornerCoord , xlim=c(45,210),ylim=c(45,160) , asp=1) ; points(coord , col="blue") ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") + expect_equivalent(projCoord, coord+60) + + # To negative coordinates : + fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) + rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) + toCornerCoord <- as.matrix(fromCornerCoord) %*% rot_mat + toCornerCoord <- sweep(toCornerCoord, 2, c(150,100), FUN = "-") + coord <- data.frame(x=c(25,60),y=c(25,40)) + #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) + projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) + #plot(fromCornerCoord , xlim=c(-200,110), ylim=c(-110,60) , asp=1) ; points(coord , col="blue") ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") + expect_equal(projCoord, data.frame(X=c(-140.8494 ,-118.0385 ),Y=c(-65.84936,-35.35898)) , tolerance = 1e-6) + +}) + diff --git a/tests/testthat/test_01_small_function.R b/tests/testthat/test_01_small_function.R index 3998d80..c1e7489 100644 --- a/tests/testthat/test_01_small_function.R +++ b/tests/testthat/test_01_small_function.R @@ -213,24 +213,6 @@ test_that("Analysis of Procrust", { }) -context("Bilinear interpolation") -test_that("Bilinear interpolation", { - relCoord <- matrix(c(0,0 , 5,5 , 10,10) , ncol = 2 , byrow = T) - cornerCoord <- data.table(expand.grid(X = c(0, 100), Y = c(0, 50)) , corner = c(1,4,2,3)) - res <- bilinearInterpolation(relCoord, cornerCoord, dimX = 20, dimY=10) - expect_is(res, "data.table") - expect_equal(nrow(res), nrow(relCoord)) - expect_equal(names(res), c("XAbs","YAbs")) - expect_equivalent(res, data.table(c(0,25,50), c(0,25,50))) - - # Test if corners are numbered in counter-clockwise direction (only when relative plot is not a square, otherwise it always works) - cornerCoord <- data.table(expand.grid(X = c(0, 100), Y = c(0, 50)) , corner = c(1,2,4,3)) - resCounterClockwise <- bilinearInterpolation(relCoord, cornerCoord = cornerCoord, dimX = 20, dimY=10) - expect_error(expect_equal(res,resCounterClockwise)) -}) - - - context("lat long to UTM") test_that("lat long to UTM", { long <- c(-52.68, -51.12, -53.11) diff --git a/vignettes/VignetteSpatialRepresentationAndManagementPlot.Rmd b/vignettes/VignetteSpatialRepresentationAndManagementPlot.Rmd index c2d1336..877cc46 100644 --- a/vignettes/VignetteSpatialRepresentationAndManagementPlot.Rmd +++ b/vignettes/VignetteSpatialRepresentationAndManagementPlot.Rmd @@ -98,7 +98,7 @@ If only **4 GPS measurements** have been taken **with a high degree of accuracy* ```{r dpi=200} check_plot <- checkPlotCoord( longlat = cornerCoord[c("Long", "Lat")], - # OR projCoord = cornerCoord[c("xProj", "yProj")], + # OR, if exists : projCoord = cornerCoord[c("xProj", "yProj")], relCoord = cornerCoord[c("xRel", "yRel")], trustGPScorners = T, cornerID = cornerCoord$Corners, drawPlot = TRUE, From 70b922c29005cf489f7a7496c1424a5f46386fed Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 24 Oct 2024 09:46:03 +0200 Subject: [PATCH 15/25] Adapted cutPlot function to bilinearInterpolation changes --- R/cutPlot.R | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/R/cutPlot.R b/R/cutPlot.R index 2b1bc2b..1d15feb 100644 --- a/R/cutPlot.R +++ b/R/cutPlot.R @@ -98,9 +98,10 @@ cutPlot <- function(projCoord, plot, cornerNum, gridsize = 100, dimX = 200, dimY )) # Transformation of relative grid coordinates into absolute coordinates - absCoord <- bilinearInterpolation(relCoord = gridMat , cornerCoord = data[,.(X,Y,cornerNum)] ,dimX = plotDimX, dimY = plotDimY ) + #absCoord <- bilinearInterpolation(relCoord = gridMat , cornerCoord = data[,.(X,Y,cornerNum)] ,dimX = plotDimX, dimY = plotDimY ) + absCoord <- bilinearInterpolation(coord = gridMat , fromCornerCoord = relCoordMat , toCornerCoord = absCoordMat ) - return(data.table(XRel = gridMat[, 1], YRel = gridMat[, 2], absCoord[, 1], absCoord[, 2])) + return(data.table(XRel = gridMat[, 1], YRel = gridMat[, 2], XAbs=absCoord[, 1], YAbs=absCoord[, 2])) } # Apply gridFunction to all plots From 161e9209b7c8858ed15dd75454998764793ede1b Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 24 Oct 2024 10:04:15 +0200 Subject: [PATCH 16/25] Modified checkPlotCoord tests to adapt them to bilinearInterpolation modifications --- tests/testthat/test-checkCoordPlot.R | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/testthat/test-checkCoordPlot.R b/tests/testthat/test-checkCoordPlot.R index 1cff25b..659f42c 100644 --- a/tests/testthat/test-checkCoordPlot.R +++ b/tests/testthat/test-checkCoordPlot.R @@ -51,7 +51,7 @@ test_that("checkPlotCoord outputs and outliers", { expect_is(outputs$outliers, "data.frame") expect_equal(dim(outputs$outliers), c(10,3)) - expect_equal(dim(outputs$cornerCoord), c(4, 6)) + expect_equal(dim(outputs$cornerCoord), c(4, 5)) # with max dist equal 25 there isn't outliers anymore expect_failure(expect_warning( @@ -102,11 +102,11 @@ test_that("checkPlotCoord, trustGPScorners", { checkPlotCoord(projCoord = projCoord[-(1:3),], relCoord = relCoord[-(1:3),], trustGPScorners = T, rmOutliers = T, cornerID = cornerID[-(1:3)], drawPlot = F),"At least one corner has less than 5 measurements. We suggest using the argument trustGPScorners = FALSE" ) res_trustGPScorners_T <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) - expect_equivalent(res_trustGPScorners_T$cornerCoord[,1:2] , unique(projCoord0)) + expect_equivalent(res_trustGPScorners_T$cornerCoord[,1:2] , unique(projCoord0)[c(1,4,3,2),]) res_trustGPScorners_F <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F, rmOutliers = T, cornerID = cornerID, drawPlot = F) - expect_equal(res_trustGPScorners_T$cornerCoord[,3:6] , res_trustGPScorners_F$cornerCoord[,3:6]) - expect_equal(res_trustGPScorners_T$cornerCoord[,1:2], round(res_trustGPScorners_F$cornerCoord[,1:2],digits = -1)) + expect_equivalent(res_trustGPScorners_T$cornerCoord[,3:5] , res_trustGPScorners_F$cornerCoord[,3:5]) + expect_equivalent(res_trustGPScorners_T$cornerCoord[,1:2], round(res_trustGPScorners_F$cornerCoord[,1:2],digits = -1)) # Test when there's only 4 measures and trustGPScorners = F expect_equal(res_trustGPScorners_F , checkPlotCoord(projCoord = projCoord0[c(1,8,15,22),], relCoord = relCoord[c(1,8,15,22),], trustGPScorners = F, rmOutliers = T, drawPlot = F, cornerID = cornerID[c(1,8,15,22)])) @@ -131,13 +131,13 @@ test_that("checkPlotCoord, origin corner", { res_NE <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord_NE, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) res_SW <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord_SW, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) - expect_equal(res_SW$cornerCoord[c("Xrel","Yrel","cornerNum")] , res_NE$cornerCoord[c("Xrel","Yrel","cornerNum")]) + expect_equivalent(res_SW$cornerCoord[c("Xrel","Yrel")] , res_NE$cornerCoord[,c("Xrel","Yrel")]) expect_equivalent(res_SW$cornerCoord[c("X","Y","cornerID")] , res_NE$cornerCoord[c(3,4,1,2),c("X","Y","cornerID")]) # Test when the origin of the relative coordinates is not(0;0) relCoord_SW <- relCoord_SW + 200 res_SW_200 <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord_SW, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) - expect_equal(res_SW$cornerCoord[c("X","Y","cornerID","cornerNum")] , res_SW_200$cornerCoord[c("X","Y","cornerID","cornerNum")]) + expect_equal(res_SW$cornerCoord[c("X","Y","cornerID")] , res_SW_200$cornerCoord[c("X","Y","cornerID")]) expect_equal(res_SW$cornerCoord[c("Xrel","Yrel")]+200 , res_SW_200$cornerCoord[c("Xrel","Yrel")]) }) From e2c86fab29a42451357b8ad8654dc15d761188c2 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Tue, 5 Nov 2024 15:18:38 +0100 Subject: [PATCH 17/25] - Changed bilinearInterpolation name function to bilinear_interpolation (and also the arguments of the fonction, eg, fromCornerCoord -> from_corner_coord,etc.) - Added tests for bilinear_interpolation - colnames of the bilinear_interpolation return are now the same as to_corner_coord inputs, or (x_interp , y_interp) is not supplied --- NAMESPACE | 4 +- R/bilinearInterpolation.R | 82 --------------- R/bilinear_interpolation.R | 90 +++++++++++++++++ R/checkPlotCoord.R | 15 +-- R/cutPlot.R | 3 +- man/bilinearInterpolation.Rd | 44 -------- man/bilinear_interpolation.Rd | 49 +++++++++ man/checkPlotCoord.Rd | 2 +- tests/testthat/test-bilinearInterpolation.R | 106 ++++++++++++-------- tests/testthat/test-checkCoordPlot.R | 2 +- 10 files changed, 213 insertions(+), 184 deletions(-) delete mode 100644 R/bilinearInterpolation.R create mode 100644 R/bilinear_interpolation.R delete mode 100644 man/bilinearInterpolation.Rd create mode 100644 man/bilinear_interpolation.Rd diff --git a/NAMESPACE b/NAMESPACE index 9e41af3..c585726 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,7 +3,7 @@ export(AGBmonteCarlo) export(attributeTree) export(attributeTreeCoord) -export(bilinearInterpolation) +export(bilinear_interpolation) export(cacheManager) export(cachePath) export(checkPlotCoord) @@ -14,6 +14,7 @@ export(correctCoordGPS) export(correctTaxo) export(createCache) export(cutPlot) +export(divide_plot) export(getBioclimParam) export(getTaxonomy) export(getWoodDensity) @@ -37,6 +38,7 @@ importFrom(data.table,is.data.table) importFrom(data.table,rbindlist) importFrom(data.table,setDF) importFrom(data.table,setDT) +importFrom(data.table,setcolorder) importFrom(data.table,setkey) importFrom(data.table,setnames) importFrom(data.table,tstrsplit) diff --git a/R/bilinearInterpolation.R b/R/bilinearInterpolation.R deleted file mode 100644 index df3ce87..0000000 --- a/R/bilinearInterpolation.R +++ /dev/null @@ -1,82 +0,0 @@ -#' Generalized bilinear interpolation of coordinates -#' -#' @description Apply a generalized bilinear interpolation to convert any coordinates from one original coordinate system to another, using the plot's 4 corner coordinates of both system. -#' -#' @details -#' The plot represented by the 4 coordinates in fromCornerCoord must have 4 right angles, i.e. a rectangular (or square) plot. -#' -#' @references C. -C. Wei and C. -H. Chen, "Generalized Bilinear Interpolation of Motion Vectors for Quad-Tree Mesh," 2008 International Conference on Intelligent Information Hiding and Multimedia Signal Processing, Harbin, China, 2008, pp. 635-638, doi: 10.1109/IIH-MSP.2008.283. -#' -#' -#' @param coord a matrix or data.frame : coordinates to be transformed, with X and Y corresponding to the first two columns -#' @param fromCornerCoord a matrix or data.frame : corner coordinates of the plot in the original coordinate system, with X and Y corresponding to the first two columns -#' @param toCornerCoord a matrix or data.frame : corner coordinates of the plot in the coordinate system to be projected, with the same line order as fromCornerCoord and , with X and Y corresponding to the first two columns -#' @param orderedCorner a logical, indicating if fromCornerCoord and toCornerCoord rows are sorted in correct order (clockwise or counter-clockwise) -#' -#' @return a data.frame containing the converted coordinates -#' -#' @importFrom data.table is.data.table -#' -#' @export -#' -#' @keywords generalized bilinear interpolation -#' -#' @author Arthur Bailly -#' -#' @examples -#' fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) -#' rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) -#' toCornerCoord <- as.matrix(fromCornerCoord) %*% rot_mat -#' toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") -#' coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) -#' projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) -#' # plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(fromCornerCoord) ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") - - -bilinearInterpolation = function(coord, fromCornerCoord, toCornerCoord, orderedCorner = F) { - - # Parameters verification - if(nrow(fromCornerCoord)!=4 | nrow(toCornerCoord)!=4 | nrow(fromCornerCoord)!=nrow(fromCornerCoord)) { - stop("fromCornerCoord and toCornerCoord must have 4 rows representing the 4 corners of the plot") - } - if(!(is.data.frame(coord) | is.matrix(coord) | is.data.table(coord))){ - stop("tree coordinates must be a data.frame, a matrix or a data.table") - } - if(is.data.table(coord)) coord <- data.frame(coord) - if(is.data.table(fromCornerCoord) | is.data.table(toCornerCoord)) { - fromCornerCoord <- data.frame(fromCornerCoord) - toCornerCoord <- data.frame(toCornerCoord) - } - - # Sorting rows if necessary - centroid <- colMeans(fromCornerCoord) - if(!orderedCorner) { - # Sort fromCornerCoord and toCornerCoord rows in a counter-clockwise direction - angles <- atan2(fromCornerCoord[, 2] - centroid[2], fromCornerCoord[, 1] - centroid[1]) - fromCornerCoord <- fromCornerCoord[order(angles), ] - toCornerCoord <- toCornerCoord[order(angles), ] - } - - # Verification of a rectangular plot for fromCornerCoord - if(!all(abs(dist(rbind(fromCornerCoord,centroid))[c(4,7,9,10)] - mean(dist(rbind(fromCornerCoord,centroid))[c(4,7,9,10)]))<0.1)) { - stop("The plot in the relative coordinate system is not a rectangle (or a square). You may consider using trustGPScorners = F") - } - - x_A <- fromCornerCoord[1,1] ; x_B <- fromCornerCoord[2,1] ; x_C <- fromCornerCoord[3,1] ; x_D <- fromCornerCoord[4,1] - y_A <- fromCornerCoord[1,2] ; y_B <- fromCornerCoord[2,2] ; y_C <- fromCornerCoord[3,2] ; y_D <- fromCornerCoord[4,2] - u_A <- toCornerCoord[1,1] ; u_B <- toCornerCoord[2,1]; u_C <- toCornerCoord[3,1]; u_D <- toCornerCoord[4,1] - v_A <- toCornerCoord[1,2] ; v_B <- toCornerCoord[2,2] ; v_C <- toCornerCoord[3,2] ; v_D <- toCornerCoord[4,2] - - apply_bilinear_interpolation <- function(x,y) { - rate_A <- (1-(x-x_A)/(x_C-x_A)) * (1-(y-y_A)/(y_C-y_A)) - rate_B <- (1-(x-x_B)/(x_D-x_B)) * (1-(y-y_B)/(y_D-y_B)) - rate_C <- (1-(x-x_C)/(x_A-x_C)) * (1-(y-y_C)/(y_A-y_C)) - rate_D <- (1-(x-x_D)/(x_B-x_D)) * (1-(y-y_D)/(y_B-y_D)) - return(data.frame( - X = rate_A*u_A + rate_B*u_B + rate_C*u_C + rate_D*u_D, - Y = rate_A*v_A + rate_B*v_B + rate_C*v_C + rate_D*v_D) - ) - } - - return(apply_bilinear_interpolation(x=coord[,1],y=coord[,2])) -} diff --git a/R/bilinear_interpolation.R b/R/bilinear_interpolation.R new file mode 100644 index 0000000..a12058f --- /dev/null +++ b/R/bilinear_interpolation.R @@ -0,0 +1,90 @@ +#' Generalized bilinear interpolation of coordinates +#' +#' @description Apply a generalized bilinear interpolation to convert any coordinates from one original coordinate system to another, using the plot's 4 corner coordinates of both system. +#' +#' @details +#' The plot represented by the 4 coordinates in from_corner_coord must have 4 right angles, i.e. a rectangular (or square) plot. +#' +#' @references C. -C. Wei and C. -H. Chen, "Generalized Bilinear Interpolation of Motion Vectors for Quad-Tree Mesh," 2008 International Conference on Intelligent Information Hiding and Multimedia Signal Processing, Harbin, China, 2008, pp. 635-638, doi: 10.1109/IIH-MSP.2008.283. +#' +#' +#' @param coord a matrix or data.frame : coordinates to be transformed, with X and Y corresponding to the first two columns +#' @param from_corner_coord a matrix or data.frame : corner coordinates of the plot in the original coordinate system, with X and Y corresponding to the first two columns +#' @param to_corner_coord a matrix or data.frame : corner coordinates of the plot in the coordinate system to be projected, with the same line order as from_corner_coord and , with X and Y corresponding to the first two columns +#' @param ordered_corner a logical, indicating if from_corner_coord and to_corner_coord rows are sorted in correct order (clockwise or counter-clockwise) +#' +#' @return a data.frame containing the converted coordinates +#' +#' @importFrom data.table is.data.table +#' +#' @export +#' +#' @keywords generalized bilinear interpolation +#' +#' @author Arthur Bailly +#' +#' @examples +#' from_corner_coord <- expand.grid(X = c(0, 100), Y = c(0, 50)) +#' rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) +#' to_corner_coord <- as.matrix(from_corner_coord) %*% rot_mat +#' to_corner_coord <- sweep(to_corner_coord, 2, c(50,100), FUN = "+") +#' coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) +#' projCoord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) +#' # plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(from_corner_coord) ; points(projCoord , col="purple") ; points(to_corner_coord, col="red") + + +bilinear_interpolation = function(coord, from_corner_coord, to_corner_coord, ordered_corner = F) { + + # Parameters verification + if(nrow(from_corner_coord)!=4 | nrow(to_corner_coord)!=4 | nrow(from_corner_coord)!=nrow(from_corner_coord)) { + stop("from_corner_coord and to_corner_coord must have 4 rows representing the 4 corners of the plot") + } + if(!(is.data.frame(coord) | is.matrix(coord) | is.data.table(coord))){ + stop("tree coordinates must be a data.frame, a matrix or a data.table") + } + if(is.data.table(coord)) coord <- data.frame(coord) + if(is.data.table(from_corner_coord) | is.data.table(to_corner_coord)) { + from_corner_coord <- data.frame(from_corner_coord) + to_corner_coord <- data.frame(to_corner_coord) + } + + # to_corner_coord colnames attribution + if(is.null(colnames(to_corner_coord))) { + to_corner_coord <- to_corner_coord[,1:2] + colnames(to_corner_coord) <- c("x_interp","y_interp") + } + + # Sorting rows if necessary + centroid <- colMeans(from_corner_coord[,1:2]) + if(!ordered_corner) { + # Sort from_corner_coord and to_corner_coord rows in a counter-clockwise direction + angles <- atan2(from_corner_coord[, 2] - centroid[2], from_corner_coord[, 1] - centroid[1]) + from_corner_coord <- from_corner_coord[order(angles), ] + to_corner_coord <- to_corner_coord[order(angles), ] + } + + # Verification of a rectangular plot for from_corner_coord + if(!all(abs(dist(rbind(from_corner_coord[,1:2],centroid))[c(4,7,9,10)] - mean(dist(rbind(from_corner_coord[,1:2],centroid))[c(4,7,9,10)]))<0.1)) { + stop("The plot in the relative coordinate system is not a rectangle (or a square). You may consider using trustGPScorners = F") + } + + x_A <- from_corner_coord[1,1] ; x_B <- from_corner_coord[2,1] ; x_C <- from_corner_coord[3,1] ; x_D <- from_corner_coord[4,1] + y_A <- from_corner_coord[1,2] ; y_B <- from_corner_coord[2,2] ; y_C <- from_corner_coord[3,2] ; y_D <- from_corner_coord[4,2] + u_A <- to_corner_coord[1,1] ; u_B <- to_corner_coord[2,1]; u_C <- to_corner_coord[3,1]; u_D <- to_corner_coord[4,1] + v_A <- to_corner_coord[1,2] ; v_B <- to_corner_coord[2,2] ; v_C <- to_corner_coord[3,2] ; v_D <- to_corner_coord[4,2] + + apply_bilinear_interpolation <- function(x,y,to_corner_coord_colnames) { + rate_A <- (1-(x-x_A)/(x_C-x_A)) * (1-(y-y_A)/(y_C-y_A)) + rate_B <- (1-(x-x_B)/(x_D-x_B)) * (1-(y-y_B)/(y_D-y_B)) + rate_C <- (1-(x-x_C)/(x_A-x_C)) * (1-(y-y_C)/(y_A-y_C)) + rate_D <- (1-(x-x_D)/(x_B-x_D)) * (1-(y-y_D)/(y_B-y_D)) + interp_df <- data.frame( + rate_A*u_A + rate_B*u_B + rate_C*u_C + rate_D*u_D, + rate_A*v_A + rate_B*v_B + rate_C*v_C + rate_D*v_D + ) + setnames(interp_df, new = to_corner_coord_colnames) + interp_df + } + + return(apply_bilinear_interpolation(x=coord[,1],y=coord[,2],to_corner_coord_colnames=colnames(to_corner_coord)[1:2])) +} diff --git a/R/checkPlotCoord.R b/R/checkPlotCoord.R index 76a5227..0774350 100644 --- a/R/checkPlotCoord.R +++ b/R/checkPlotCoord.R @@ -18,7 +18,7 @@ #' @param maxDist a numeric giving the maximum distance (in meters) above which GPS measurements should be considered outliers (default 15 m) #' @param rmOutliers a logical indicating if detected outliers are removed from the coordinate calculation #' @param drawPlot a logical indicating if the plot design should be displayed and returned -#' @param treeCoord a data frame containing at least the tree coordinates in the relative coordinates system (that of the field), with X and Y corresponding to the first and second columns respectively +#' @param treeCoord (optional) a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively #' #' @author Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY #' @@ -254,15 +254,10 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc warning("Be careful, one or more trees are not inside the plot defined by relCoord") } if(trustGPScorners) { - treeProjCoord <- bilinearInterpolation(coord = treeCoord[,1:2], - fromCornerCoord = cornerCoord[,c("Xrel","Yrel")], - toCornerCoord = cornerCoord[,c("X","Y")], - orderedCorner = T) - # treeProjCoord <- bilinearInterpolation(relCoord = treeCoord[,1:2], - # cornerCoord = cornerCoord[,c("X","Y","cornerNum")], - # dimX = diff(range(cornerCoord$Xrel)), - # dimY = diff(range(cornerCoord$Yrel)) - # ) + treeProjCoord <- bilinear_interpolation(coord = treeCoord[,1:2], + from_corner_coord = cornerCoord[,c("Xrel","Yrel")], + to_corner_coord = cornerCoord[,c("X","Y")], + ordered_corner = T) colnames(treeProjCoord) <- c("X","Y") } else { treeProjCoord <- as.matrix(treeCoord[,1:2]) %*% procrustRes$rotation diff --git a/R/cutPlot.R b/R/cutPlot.R index 1d15feb..389316e 100644 --- a/R/cutPlot.R +++ b/R/cutPlot.R @@ -98,8 +98,7 @@ cutPlot <- function(projCoord, plot, cornerNum, gridsize = 100, dimX = 200, dimY )) # Transformation of relative grid coordinates into absolute coordinates - #absCoord <- bilinearInterpolation(relCoord = gridMat , cornerCoord = data[,.(X,Y,cornerNum)] ,dimX = plotDimX, dimY = plotDimY ) - absCoord <- bilinearInterpolation(coord = gridMat , fromCornerCoord = relCoordMat , toCornerCoord = absCoordMat ) + absCoord <- bilinear_interpolation(coord = gridMat , from_corner_coord = relCoordMat , to_corner_coord = absCoordMat ) return(data.table(XRel = gridMat[, 1], YRel = gridMat[, 2], XAbs=absCoord[, 1], YAbs=absCoord[, 2])) } diff --git a/man/bilinearInterpolation.Rd b/man/bilinearInterpolation.Rd deleted file mode 100644 index 2115a17..0000000 --- a/man/bilinearInterpolation.Rd +++ /dev/null @@ -1,44 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/bilinearInterpolation.R -\name{bilinearInterpolation} -\alias{bilinearInterpolation} -\title{Generalized bilinear interpolation of coordinates} -\usage{ -bilinearInterpolation(coord, fromCornerCoord, toCornerCoord, orderedCorner = F) -} -\arguments{ -\item{coord}{a matrix or data.frame : coordinates to be transformed, with X and Y corresponding to the first two columns} - -\item{fromCornerCoord}{a matrix or data.frame : corner coordinates of the plot in the original coordinate system, with X and Y corresponding to the first two columns} - -\item{toCornerCoord}{a matrix or data.frame : corner coordinates of the plot in the coordinate system to be projected, with the same line order as fromCornerCoord and , with X and Y corresponding to the first two columns} - -\item{orderedCorner}{a logical, indicating if fromCornerCoord and toCornerCoord rows are sorted in correct order (clockwise or counter-clockwise)} -} -\value{ -a data.frame containing the converted coordinates -} -\description{ -Apply a generalized bilinear interpolation to convert any coordinates from one original coordinate system to another, using the plot's 4 corner coordinates of both system. -} -\details{ -The plot represented by the 4 coordinates in fromCornerCoord must have 4 right angles, i.e. a rectangular (or square) plot. -} -\examples{ -fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) -rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) -toCornerCoord <- as.matrix(fromCornerCoord) \%*\% rot_mat -toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") -coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) -projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) -# plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(fromCornerCoord) ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") -} -\references{ -C. -C. Wei and C. -H. Chen, "Generalized Bilinear Interpolation of Motion Vectors for Quad-Tree Mesh," 2008 International Conference on Intelligent Information Hiding and Multimedia Signal Processing, Harbin, China, 2008, pp. 635-638, doi: 10.1109/IIH-MSP.2008.283. -} -\author{ -Arthur Bailly -} -\keyword{bilinear} -\keyword{generalized} -\keyword{interpolation} diff --git a/man/bilinear_interpolation.Rd b/man/bilinear_interpolation.Rd new file mode 100644 index 0000000..a64210e --- /dev/null +++ b/man/bilinear_interpolation.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/bilinear_interpolation.R +\name{bilinear_interpolation} +\alias{bilinear_interpolation} +\title{Generalized bilinear interpolation of coordinates} +\usage{ +bilinear_interpolation( + coord, + from_corner_coord, + to_corner_coord, + ordered_corner = F +) +} +\arguments{ +\item{coord}{a matrix or data.frame : coordinates to be transformed, with X and Y corresponding to the first two columns} + +\item{from_corner_coord}{a matrix or data.frame : corner coordinates of the plot in the original coordinate system, with X and Y corresponding to the first two columns} + +\item{to_corner_coord}{a matrix or data.frame : corner coordinates of the plot in the coordinate system to be projected, with the same line order as from_corner_coord and , with X and Y corresponding to the first two columns} + +\item{ordered_corner}{a logical, indicating if from_corner_coord and to_corner_coord rows are sorted in correct order (clockwise or counter-clockwise)} +} +\value{ +a data.frame containing the converted coordinates +} +\description{ +Apply a generalized bilinear interpolation to convert any coordinates from one original coordinate system to another, using the plot's 4 corner coordinates of both system. +} +\details{ +The plot represented by the 4 coordinates in from_corner_coord must have 4 right angles, i.e. a rectangular (or square) plot. +} +\examples{ +from_corner_coord <- expand.grid(X = c(0, 100), Y = c(0, 50)) +rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) +to_corner_coord <- as.matrix(from_corner_coord) \%*\% rot_mat +to_corner_coord <- sweep(to_corner_coord, 2, c(50,100), FUN = "+") +coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) +projCoord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) +# plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(from_corner_coord) ; points(projCoord , col="purple") ; points(to_corner_coord, col="red") +} +\references{ +C. -C. Wei and C. -H. Chen, "Generalized Bilinear Interpolation of Motion Vectors for Quad-Tree Mesh," 2008 International Conference on Intelligent Information Hiding and Multimedia Signal Processing, Harbin, China, 2008, pp. 635-638, doi: 10.1109/IIH-MSP.2008.283. +} +\author{ +Arthur Bailly +} +\keyword{bilinear} +\keyword{generalized} +\keyword{interpolation} diff --git a/man/checkPlotCoord.Rd b/man/checkPlotCoord.Rd index 344238c..aee5517 100644 --- a/man/checkPlotCoord.Rd +++ b/man/checkPlotCoord.Rd @@ -33,7 +33,7 @@ checkPlotCoord( \item{drawPlot}{a logical indicating if the plot design should be displayed and returned} -\item{treeCoord}{a data frame containing at least the tree coordinates in the relative coordinates system (that of the field), with X and Y corresponding to the first and second columns respectively} +\item{treeCoord}{(optional) a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively} } \value{ Returns a list including : diff --git a/tests/testthat/test-bilinearInterpolation.R b/tests/testthat/test-bilinearInterpolation.R index 37bf6da..1517407 100644 --- a/tests/testthat/test-bilinearInterpolation.R +++ b/tests/testthat/test-bilinearInterpolation.R @@ -1,71 +1,91 @@ context("Bilinear interpolation") -test_that("bilinearInterpolation error", { - fromCornerCoord <- data.frame(x=1:5,y=1:5) - toCornerCoord <- data.frame(x=1:4,y=1:4) - expect_error(bilinearInterpolation(coord = NULL, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) , "fromCornerCoord and toCornerCoord must have 4 rows representing the 4 corners of the plot") - expect_error(bilinearInterpolation(coord = NULL, fromCornerCoord = toCornerCoord, toCornerCoord = fromCornerCoord) , "fromCornerCoord and toCornerCoord must have 4 rows representing the 4 corners of the plot") - fromCornerCoord <- data.frame(x=c(0,0,10,10), y=c(0,10,10,0)) - toCornerCoord <- fromCornerCoord+50 +test_that("bilinear_interpolation error", { + from_corner_coord <- data.frame(x=1:5,y=1:5) + to_corner_coord <- data.frame(Xproj=1:4,Yproj=1:4) + expect_error(bilinear_interpolation(coord = NULL, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) , "from_corner_coord and to_corner_coord must have 4 rows representing the 4 corners of the plot") + expect_error(bilinear_interpolation(coord = NULL, from_corner_coord = to_corner_coord, to_corner_coord = from_corner_coord) , "from_corner_coord and to_corner_coord must have 4 rows representing the 4 corners of the plot") + from_corner_coord <- data.frame(x=c(0,0,10,10), y=c(0,10,10,0)) + to_corner_coord <- from_corner_coord+50 coord <- c(5,5) - expect_error(bilinearInterpolation(coord = coord, fromCornerCoord = toCornerCoord, toCornerCoord = fromCornerCoord) , "tree coordinates must be a data.frame, a matrix or a data.table") + expect_error(bilinear_interpolation(coord = coord, from_corner_coord = to_corner_coord, to_corner_coord = from_corner_coord) , "tree coordinates must be a data.frame, a matrix or a data.table") coord <- data.frame(x=5,y=5) - fromCornerCoord[3,1] <- 11 - expect_error(bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) , "You may consider using trustGPScorners = F") + # Error when from_corner_coord is not a rectangle + from_corner_coord[3,1] <- 11 + expect_error(bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) , "You may consider using trustGPScorners = F") }) -test_that("bilinearInterpolation function", { +test_that("bilinear_interpolation function", { - fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) + # Test when matrices are supplied and to_corner_coord colnames are not + from_corner_coord <- as.matrix(expand.grid(X = c(0, 100), Y = c(0, 50))) rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) - toCornerCoord <- as.matrix(fromCornerCoord) %*% rot_mat - toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") + to_corner_coord <- as.matrix(from_corner_coord) %*% rot_mat + to_corner_coord <- sweep(to_corner_coord, 2, c(50,100), FUN = "+") + coord <- as.matrix(data.frame(x=c(25,60),y=c(25,40))) + proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) - coord <- data.frame(x=c(25,60),y=c(25,40)) - projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) - #plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(fromCornerCoord) ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") - - expect_is(projCoord,"data.frame") - expect_equal(nrow(projCoord), nrow(coord)) - expect_equal(names(projCoord), c("X","Y")) + #plot(coord, xlim=c(-10,150),ylim=c(-5,200), col="blue") ; points(from_corner_coord) ; points(proj_coord , col="purple") ; points(to_corner_coord, col="red") - expect_equal(projCoord, data.frame(X=c(59.15064,81.96152),Y=c(134.1506,164.6410)) , tolerance = 1e-6) + expect_is(proj_coord,"data.frame") + expect_equal(nrow(proj_coord), nrow(coord)) + expect_equal(names(proj_coord), c("Xinterp","Yinterp")) + expect_equal(proj_coord, data.frame(Xinterp=c(59.15064,81.96152),Yinterp=c(134.1506,164.6410)) , tolerance = 1e-6) + # Test when to_corner_coord colnames are supplied + colnames(to_corner_coord) <- c("Xproj","Yproj") + proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) + expect_equal(names(proj_coord), c("Xproj","Yproj")) # Corner deformation set.seed(52) - toCornerCoord <- toCornerCoord + runif(1,-10,10) - # toCornerCoord <- fromCornerCoord - # toCornerCoord[2,] = c(103,-5) - # toCornerCoord[3,] = c(10,52) - # toCornerCoord[4,] = c(95,48) - # toCornerCoord <- sweep(toCornerCoord, 2, c(50,100), FUN = "+") - projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) - expect_equal(projCoord, data.frame(X=c(52.39103,75.20192),Y=c(127.3910,157.8814)) , tolerance = 1e-6) + to_corner_coord <- to_corner_coord + runif(1,-10,10) + # to_corner_coord <- from_corner_coord + # to_corner_coord[2,] = c(103,-5) + # to_corner_coord[3,] = c(10,52) + # to_corner_coord[4,] = c(95,48) + # to_corner_coord <- sweep(to_corner_coord, 2, c(50,100), FUN = "+") + proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) + expect_equal(proj_coord, data.frame(Xproj=c(52.39103,75.20192),Yproj=c(127.3910,157.8814)) , tolerance = 1e-6) #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) - #plot(fromCornerCoord, xlim=c(-10,150),ylim=c(-5,200)) ; points(coord , col="blue") ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") - + #plot(from_corner_coord, xlim=c(-10,150),ylim=c(-5,200)) ; points(coord , col="blue") ; points(proj_coord , col="purple") ; points(to_corner_coord, col="red") # If the origin is not (0;0) - fromCornerCoord <- expand.grid(X = c(50, 150), Y = c(50, 100)) - toCornerCoord <- sweep(fromCornerCoord, 2, c(60,60), FUN = "+") + from_corner_coord <- expand.grid(X = c(50, 150), Y = c(50, 100)) + to_corner_coord <- sweep(from_corner_coord, 2, c(60,60), FUN = "+") coord <- data.frame(x=c(25,60)+50,y=c(25,40)+50) #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) + 50 - projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) - #plot(fromCornerCoord , xlim=c(45,210),ylim=c(45,160) , asp=1) ; points(coord , col="blue") ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") - expect_equivalent(projCoord, coord+60) + proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) + #plot(from_corner_coord , xlim=c(45,210),ylim=c(45,160) , asp=1) ; points(coord , col="blue") ; points(proj_coord , col="purple") ; points(to_corner_coord, col="red") + expect_equivalent(proj_coord, coord+60) # To negative coordinates : - fromCornerCoord <- expand.grid(X = c(0, 100), Y = c(0, 50)) + from_corner_coord <- expand.grid(X = c(0, 100), Y = c(0, 50)) rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) - toCornerCoord <- as.matrix(fromCornerCoord) %*% rot_mat - toCornerCoord <- sweep(toCornerCoord, 2, c(150,100), FUN = "-") + to_corner_coord <- as.matrix(from_corner_coord) %*% rot_mat + to_corner_coord <- sweep(to_corner_coord, 2, c(150,100), FUN = "-") coord <- data.frame(x=c(25,60),y=c(25,40)) #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) - projCoord = bilinearInterpolation(coord = coord, fromCornerCoord = fromCornerCoord, toCornerCoord = toCornerCoord) - #plot(fromCornerCoord , xlim=c(-200,110), ylim=c(-110,60) , asp=1) ; points(coord , col="blue") ; points(projCoord , col="purple") ; points(toCornerCoord, col="red") - expect_equal(projCoord, data.frame(X=c(-140.8494 ,-118.0385 ),Y=c(-65.84936,-35.35898)) , tolerance = 1e-6) + proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) + #plot(from_corner_coord , xlim=c(-200,110), ylim=c(-110,60) , asp=1) ; points(coord , col="blue") ; points(proj_coord , col="purple") ; points(to_corner_coord, col="red") + expect_equal(proj_coord, data.frame(Xinterp=c(-140.8494 ,-118.0385 ),Yinterp=c(-65.84936,-35.35898)) , tolerance = 1e-6) + + # If the origin is the NE corner : + to_corner_coord <- sweep(from_corner_coord, 2, c(60,60), FUN = "+") + from_corner_coord <- from_corner_coord[c(4,3,2,1),] + coord <- data.frame(x=c(25,60),y=c(25,40)) + proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) + expect_equal(proj_coord, data.frame(X=c(135,100 ),Y=c(85,70))) + + # Test when there are multiple columns + coord <- data.frame(coord) + coord$col <- rep("",2) + from_corner_coord <- data.frame(from_corner_coord) + from_corner_coord$col <- rep("",4) + to_corner_coord <- data.frame(to_corner_coord) + to_corner_coord$col <- rep("",4) + proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) }) diff --git a/tests/testthat/test-checkCoordPlot.R b/tests/testthat/test-checkCoordPlot.R index 659f42c..8579e9b 100644 --- a/tests/testthat/test-checkCoordPlot.R +++ b/tests/testthat/test-checkCoordPlot.R @@ -119,7 +119,7 @@ test_that("checkPlotCoord, origin corner", { ) cornerID <- rep(c("SW","NW","SE","NE"),e=7) - # Test when the origin is not the South-East corner + # Test when the origin is not the South-West corner relCoord_NE <- data.frame( X = rep(c(100,100,0,0),e=7), Y = rep(c(100,0,0,100),e=7) From 23c3ccf7425a0a1f84a2c71392e0419a3be47692 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Tue, 5 Nov 2024 18:50:17 +0100 Subject: [PATCH 18/25] Transformed cutPlot into divide_plot. The function now : - divide a plot in the relative coordinates system by taking the coordinates of the corners - if corner projected coordinates are supplied, a bilinear interpolation is made on subcorner to calcul their coordinates in the projected system - if dealing with multiple plot, the user has to manage it with the plot_ID_corner argument - no corner numbering is needed anymore - no dimX and dimY are needed neither - the function also replace the attributeTree function by supplying the tree_coord agrument (and the plot_ID_tree argument if dealing with multiple plots) --- R/divide_plot.R | 191 +++++++++++++++++++++++++++++++ man/divide_plot.Rd | 61 ++++++++++ tests/testthat/test-dividePlot.R | 177 ++++++++++++++++++++++++++++ 3 files changed, 429 insertions(+) create mode 100644 R/divide_plot.R create mode 100644 man/divide_plot.Rd create mode 100644 tests/testthat/test-dividePlot.R diff --git a/R/divide_plot.R b/R/divide_plot.R new file mode 100644 index 0000000..fd78d65 --- /dev/null +++ b/R/divide_plot.R @@ -0,0 +1,191 @@ +#' Divides one ore more plots into subplots +#' +#' This function divides a plot (or several plots) into subplots in the relative coordinates system, and returns the coordinates of subplot corners +#' If corner coordinates in the projected coordinate system is supplied (proj_coord), subplot projected coordinates are calculated by a bilinear interpolation in relation with plot corner relative coordinates. +#' +#' @param rel_coord a data frame containing the relative (local) coordinates of plot corners, with X and Y on the first and second column respectively +#' @param proj_coord (optional) a data frame containing the projected coordinates of plot corners, with X and Y on the first and second column respectively, and with the same row order than rel_coord +#' @param grid_size a vector indicating the dimensions of grid cells (resp. X and Y dimensions). If only one value is given, the grid cells will be considered as squares. +#' @param plot_ID_corner if dealing with multiple plots : a vector indicating plot IDs for corners. +#' @param tree_coord (otpional) a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively +#' @param plot_ID_tree if dealing with multiple plots : a vector indicating tree plot IDs. +#' +#' @return Returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : +#' - `plot_ID_corner`: If dealing with multiple plots : the plot code +#' - `subplot`: The automatically generated subplot code +#' - `Xrel`: The relative X-axis coordinates of subplots corners +#' - `Yrel`: The relative Y-axis coordinates of subplots corners +#' - `X`: If proj_coord is supplied : the projected X-axis coordinates of subplots corners +#' - `Y`: If proj_coord is supplied : the projected Y-axis coordinates of subplots corners +#' +#' @export +#' @author Arthur PERE, Arthur BAILLY +#' @importFrom data.table data.table := setnames setcolorder +#' @examples +#' +#' coord <- data.frame(X = c(0, 200, 0, 200), Y = c(0, 0, 200, 200)) + 5000 +#' cornerNum <- c(1, 2, 4, 3) +#' plot <- rep("plot1", 4) +#' +#' cut <- cutPlot(coord, plot, cornerNum, grid_size = 100, dimX = 200, dimY = 200) +#' +#' # plot the result +#' plot(coord, main = "example", xlim = c(4900, 5300), ylim = c(4900, 5300), asp = 1) +#' text(coord, labels = cornerNum, pos = 1) +#' points(cut$XAbs, cut$YAbs, pch = "+") +#' legend("bottomright", legend = c("orignal", "cut"), pch = c("o", "+")) +#' +divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner = NULL, tree_coord = NULL, plot_ID_tree = NULL) { + + # Parameters verification ---------------------------------------------------- + if (is.matrix(rel_coord)) { + rel_coord <- data.frame(x=rel_coord[,1],y=rel_coord[,2]) + } + if (!is.null(proj_coord) && !is.data.frame(proj_coord)) { + proj_coord <- data.frame(x_proj=proj_coord[,1], y_proj=proj_coord[,2]) + } + if (!is.null(plot_ID_corner) && nrow(rel_coord) != length(plot_ID_corner)) { + stop("The length of plot_ID_corner and the number of rows of rel_coord are different") + } + if(!length(grid_size) %in% c(1,2)) { + stop("The length of grid_size must be equal to 1 or 2\nIf you want to divide several plots with different grid sizes, you must apply yourself the function for each plot") + } + if(nrow(rel_coord)%% 4 !=0){ + stop("rel_coord does'nt contain exactly 4 corners by plot") + } + if(nrow(rel_coord)!=4 & is.null(plot_ID_corner)){ + stop("You must supply plot_ID_corner if you have more than one plot in your data") + } + if(!is.null(plot_ID_corner) && nrow(rel_coord)/length(unique(plot_ID_corner)) !=4 ) { + stop("plot_ID_corner must contain the same number of plots as rel_coord") + } + if(!is.null(proj_coord) && any(dim(proj_coord) != dim(rel_coord))) { + stop("rel_coord and proj_coord are not of the same dimension") + } + if(!is.null(tree_coord) && (!is.matrix(tree_coord) & !is.data.frame(tree_coord))){ + stop("tree_coord must be a matrix or a data frame") + } + if(nrow(rel_coord)!=4 & !is.null(tree_coord) & is.null(plot_ID_tree)){ + stop("You must supply plot_ID_tree if you have more than one plot in your data") + } + if(!is.null(tree_coord) && !is.null(plot_ID_tree) && nrow(tree_coord)!=length(plot_ID_tree)) { + stop("The length of plot_ID_tree and the number of rows of tree_coord are different") + } + if(nrow(rel_coord)!=4 && !is.null(plot_ID_corner) && !is.null(plot_ID_tree) && any(! unique(plot_ID_tree) %in% unique(plot_ID_corner))) { + stop("Some plot_ID_tree are not found in plot_ID_corner") + } + + + # Formatting data ----------------------------------------------------------- + + if(length(grid_size)!=2) grid_size = rep(grid_size,2) + + rel_coord <- as.data.table(rel_coord) + + if(!is.null(plot_ID_corner)){ + rel_coord[,plot_ID_corner:=plot_ID_corner] + } else{ + rel_coord[,plot_ID_corner:=""] + } + setcolorder(rel_coord , "plot_ID_corner" , after = 2) # if rel_coord has more than 2 col + + # Sorting rows in a counter-clockwise direction and check for non-rectangular plot + sort_rows <- function(dat) { + centroid <- colMeans(dat[,1:2]) + angles <- atan2(dat[[2]] - centroid[2], dat[[1]] - centroid[1]) + dat <- dat[order(angles), ] + # check for non-rectangular plot : distances between centroid and corners must be equals + if(!all(abs(dist(rbind(dat[,1:2],as.data.frame.list(centroid)))[c(4,7,9,10)] - mean(dist(rbind(dat[,1:2],as.data.frame.list(centroid)))[c(4,7,9,10)]))<0.1)) { + stop("The plot in the relative coordinate system is not a rectangle (or a square). BIOMASS package can't deal with non-rectangular plot") + } + return(dat) + } + + if(!is.null(proj_coord)) { + combine_coord <- cbind(rel_coord[,1:3] , proj_coord) + combine_coord <- combine_coord[, sort_rows(.SD), by = plot_ID_corner] + } else { + combine_coord <- rel_coord[, sort_rows(.SD), by = plot_ID_corner] + } + + # Dividing plots ----------------------------------------------------------- + + # Grids the plot from the relative coordinates and calculates the projected coordinates of the grid points. + divide_plot_fct <- function(dat, grid_size) { + + # Controlling for a regular grid + if((max(dat[[2]])-min(dat[[2]])) %% grid_size[1] != 0) { + stop("The x-dimension of the plot is not a multiple of the x-dimension of the grid size") + } + if((max(dat[[3]])-min(dat[[3]])) %% grid_size[2] != 0) { + stop("The y-dimension of the plot is not a multiple of the y-dimension of the grid size") + } + + plot_grid <- data.table(as.matrix(expand.grid( + Xrel = seq(min(dat[[2]]), max(dat[[2]]), by = grid_size[1]), + Yrel = seq(min(dat[[3]]), max(dat[[3]]), by = grid_size[2]) + ))) + plot_grid[,plot_ID_corner:=unique(dat$plot_ID_corner)] + + # Attributing subplots names to each corner and adding shared subplot corners + plot_grid <- rbindlist(apply(plot_grid[Xrel < max(Xrel) & Yrel < max(Yrel),], 1, function(grid_dat) { + X <- as.numeric(grid_dat[["Xrel"]]) + Y <- as.numeric(grid_dat[["Yrel"]]) + plot_grid[ + (Xrel == X & Yrel == Y) | (Xrel == X + grid_size[1] & Yrel == Y) | (Xrel == X + grid_size[1] & Yrel == Y + grid_size[2]) | (Xrel == X & Yrel == Y + grid_size[2]), + .(subplot_id = paste(plot_ID_corner, (X-min(plot_grid$Xrel)) / grid_size[1], (Y-min(plot_grid$Yrel)) / grid_size[2], sep = "_"),Xrel, Yrel)] + })) + setnames(plot_grid,2:3,names(combine_coord)[2:3]) + + # Transformation of relative grid coordinates into projected coordinates if supplied + if(!is.null(proj_coord)) { + plot_grid <- cbind(plot_grid,bilinear_interpolation(coord = plot_grid[,2:3] , from_corner_coord = dat[,2:3] , to_corner_coord = dat[,4:5], ordered_corner = T)) + } + return(plot_grid) + } + + # Apply divide_plot_fct to all plots + sub_corner_coord <- combine_coord[, divide_plot_fct(.SD, grid_size), by = plot_ID_corner, .SDcols = colnames(combine_coord)] + + # Assigning trees to subplots ------------------------------------------------ + if(!is.null(tree_coord)) { + tree_coord <- as.data.table(tree_coord) + if(!is.null(plot_ID_tree)){ + tree_coord[,plot_ID_tree:=plot_ID_tree] + } else{ + tree_coord[,plot_ID_tree:=""] + } + + invisible(lapply(split(sub_corner_coord, by = "subplot_id", keep.by = TRUE), function(dat) { + tree_coord[ plot_ID_tree == dat$plot_ID_corner[1] & + get(names(tree_coord)[1]) %between% range(dat[[3]]) & + get(names(tree_coord)[2]) %between% range(dat[[4]]), + subplot_id := dat$subplot_id[1] ] + })) + + if (anyNA(tree_coord[, subplot_id])) { + warning("One or more trees could not be assigned to a subplot (not in a subplot area)") + } + if(is.null(plot_ID_tree)) { + tree_coord[ , c("subplot_id","plot_ID_tree") := list(paste0("subplot",subplot_id),NULL)] + } + } + + + # returns -------------------------------------------------------------------- + + if(is.null(plot_ID_corner)) { + sub_corner_coord[ , c("subplot_id","plot_ID_corner") := list(paste0("subplot",subplot_id),NULL)] + } + + + if(is.null(tree_coord)) { + output <- data.frame(sub_corner_coord) + } else { + output <- list(sub_corner_coord = data.frame(sub_corner_coord), + tree_coord = tree_coord) + } + + return(output) +} + diff --git a/man/divide_plot.Rd b/man/divide_plot.Rd new file mode 100644 index 0000000..b3300ff --- /dev/null +++ b/man/divide_plot.Rd @@ -0,0 +1,61 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/divide_plot.R +\name{divide_plot} +\alias{divide_plot} +\title{Divides one ore more plots into subplots} +\usage{ +divide_plot( + rel_coord, + proj_coord = NULL, + grid_size, + plot_ID_corner = NULL, + tree_coord = NULL, + plot_ID_tree = NULL +) +} +\arguments{ +\item{rel_coord}{a data frame containing the relative (local) coordinates of plot corners, with X and Y on the first and second column respectively} + +\item{proj_coord}{(optional) a data frame containing the projected coordinates of plot corners, with X and Y on the first and second column respectively, and with the same row order than rel_coord} + +\item{grid_size}{a vector indicating the dimensions of grid cells (resp. X and Y dimensions). If only one value is given, the grid cells will be considered as squares.} + +\item{plot_ID_corner}{if dealing with multiple plots : a vector indicating plot IDs for corners.} + +\item{tree_coord}{(otpional) a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively} + +\item{plot_ID_tree}{if dealing with multiple plots : a vector indicating tree plot IDs.} +} +\value{ +Returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : +\itemize{ +\item \code{plot_ID_corner}: If dealing with multiple plots : the plot code +\item \code{subplot}: The automatically generated subplot code +\item \code{Xrel}: The relative X-axis coordinates of subplots corners +\item \code{Yrel}: The relative Y-axis coordinates of subplots corners +\item \code{X}: If proj_coord is supplied : the projected X-axis coordinates of subplots corners +\item \code{Y}: If proj_coord is supplied : the projected Y-axis coordinates of subplots corners +} +} +\description{ +This function divides a plot (or several plots) into subplots in the relative coordinates system, and returns the coordinates of subplot corners +If corner coordinates in the projected coordinate system is supplied (proj_coord), subplot projected coordinates are calculated by a bilinear interpolation in relation with plot corner relative coordinates. +} +\examples{ + +coord <- data.frame(X = c(0, 200, 0, 200), Y = c(0, 0, 200, 200)) + 5000 +cornerNum <- c(1, 2, 4, 3) +plot <- rep("plot1", 4) + +cut <- cutPlot(coord, plot, cornerNum, grid_size = 100, dimX = 200, dimY = 200) + +# plot the result +plot(coord, main = "example", xlim = c(4900, 5300), ylim = c(4900, 5300), asp = 1) +text(coord, labels = cornerNum, pos = 1) +points(cut$XAbs, cut$YAbs, pch = "+") +legend("bottomright", legend = c("orignal", "cut"), pch = c("o", "+")) + +} +\author{ +Arthur PERE, Arthur BAILLY +} diff --git a/tests/testthat/test-dividePlot.R b/tests/testthat/test-dividePlot.R new file mode 100644 index 0000000..a61a691 --- /dev/null +++ b/tests/testthat/test-dividePlot.R @@ -0,0 +1,177 @@ +context("divide Plot") + +test_that("divide_plot error", { + + rel_coord <- expand.grid(X = c(0, 100), Y = c(0, 100)) + + expect_error(divide_plot(rel_coord, grid_size = 25, plot_ID_corner = "plot"), "The length of plot_ID_corner and the number of rows of rel_coord are different") + expect_error(divide_plot(rel_coord, grid_size = c(25,25,25)), "you must apply yourself the function for each plot") + expect_error(divide_plot(rbind(rel_coord,c(0,0)),grid_size = 25), "rel_coord does'nt contain exactly 4 corners by plot") + expect_error(divide_plot(rbind(rel_coord,rel_coord),grid_size = 25), "You must supply plot_ID_corner if you have more than one plot in your data") + expect_error(divide_plot(rbind(rel_coord,rel_coord), grid_size = 25, plot_ID_corner = rep(1:4,2)), "plot_ID_corner must contain the same number of plots as rel_coord") + + corner_proj_coord <- sweep(rel_coord[1:3,], 2, c(150,150), "+") + expect_error(divide_plot(rel_coord, proj_coord = corner_proj_coord ,grid_size = 25), "rel_coord and proj_coord are not of the same dimension") + + expect_error(divide_plot(rel_coord, grid_size = 25, tree_coord = c(10,10)), "tree_coord must be a matrix or a data frame") + + expect_error(divide_plot(rel_coord, grid_size = c(30,25)) , "The x-dimension of the plot is not a multiple of the x-dimension of the grid size") + expect_error(divide_plot(rel_coord, grid_size = c(25,30)) , "The y-dimension of the plot is not a multiple of the y-dimension of the grid size") + + tree_coord <- data.frame(x =runif(100,0,100), y=runif(100,0,100)) + expect_error(divide_plot(rbind(rel_coord,rel_coord), grid_size = 25, plot_ID_corner = rep(1:2,e=4) , tree_coord = tree_coord), "You must supply plot_ID_tree if you have more than one plot in your data") + + # when the plot is not a rectangle + rel_coord <- data.frame(X=c(0,100,0,110),y=c(0,0,100,100)) + expect_error(divide_plot(rel_coord, grid_size = 25) , "BIOMASS package can't deal with non-rectangular plot") + +}) + + +test_that("divide_plot on relative coordinates only", { + + # Test when rel_coord is a matrix (colnames aren't supplied) + rel_coord <- matrix(c(0,100,0,100,0,0,100,100), ncol=2) + subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x=c(0,50,0,50), y=c(0,0,50,50))) + + # Test when rel_coord is a data.table + rel_coord <- data.table(expand.grid(x_field = c(0, 100), y_field = c(0, 100))) + subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x_field=c(0,50,0,50), y_field=c(0,0,50,50))) + + # Test when rel_coord is a data.frame + rel_coord <- expand.grid(x_field = c(0, 100), y_field = c(0, 100)) + subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x_field=c(0,50,0,50), y_field=c(0,0,50,50))) + + # Test rectangular division + rect_subplots <- divide_plot(rel_coord, grid_size = c(25,50)) + expect_equivalent(rect_subplots[29:32,] , data.frame(subplot_id=rep("subplot_3_1",4),x_field=c(75,100,75,100), y_field=c(50,50,100,100))) + + # Test when the origin is not (0;0) + rel_coord_2 <- expand.grid(x_field = c(10, 110), y_field = c(10, 110)) + subplots <- divide_plot(rel_coord = rel_coord_2, grid_size = 50) + expect_equivalent(subplots[13:16,] , data.frame(subplot_id=rep("subplot_1_1",4),x=c(60,110,60,110), y=c(60,60,110,110))) + + # Test multiple plots + multiple_rel_coord <- rbind(rel_coord , rel_coord_2) + multiple_rel_coord$plotID <- rep(c("plot1","plot2"),e=4) + multiple_subplots <- divide_plot(rel_coord = multiple_rel_coord[,c("x_field","y_field")], grid_size = 50, plot_ID_corner = multiple_rel_coord$plotID) + expect_equivalent(multiple_subplots[29:32,] , data.frame(plotID=rep("plot2",4),subplot_id=rep("plot2_1_1",4),x_field=c(60,110,60,110), y_field=c(60,60,110,110))) + + # Test rectangular plot + rel_coord <- expand.grid(x = c(0, 100), y = c(0, 50)) + subplots <- divide_plot(rel_coord = rel_coord, grid_size = c(50,25)) + #ggplot(subplots,aes(x=x,y=y,label=subplot_id)) + geom_point() + geom_text(position = position_jitter()) + coord_equal() + expect_equivalent(subplots[5:8,] , data.frame(subplot_id=rep("subplot_1_0",4),x=c(50,100,50,100), y=c(0,0,25,25))) +}) + +test_that("divide_plot with projected coordinates", { + + # Test when proj_coord is a matrix (colnames aren't supplied) + rel_coord <- expand.grid(x = c(0, 100), y = c(0, 100)) + rot_mat <- matrix(c(cos(-pi/6),sin(-pi/6),-sin(-pi/6),cos(-pi/6)),nrow=2) + proj_coord <- as.matrix(rel_coord) %*% rot_mat + proj_coord <- sweep(proj_coord, 2, c(110,110), FUN = "+") + + subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50, proj_coord = proj_coord) + #ggplot(subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x,y=y)) + geom_text(position = position_jitter()) + coord_equal() + expect_equivalent(subplots[13:16,c("subplot_id","x_proj","y_proj")] , + data.frame(subplot_id="subplot_1_1", + x_proj=c(128.3013,171.6025,103.3013,146.6025), + y_proj=c(178.3013,203.3013,221.6025,246.6025)), + tol = 1e-5) + + # Test when proj_coord is a data.table + proj_coord <- as.data.table(proj_coord) + colnames(proj_coord) <- c("Xproj","Yproj") + subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50, proj_coord = proj_coord) + expect_equivalent(subplots[13:16,c("subplot_id","Xproj","Yproj")] , + data.frame(subplot_id="subplot_1_1", + Xproj=c(128.3013,171.6025,103.3013,146.6025), + Yproj=c(178.3013,203.3013,221.6025,246.6025)), + tol = 1e-5) + + # Test when proj_coord is a data.frame + proj_coord <- as.data.frame(proj_coord) + colnames(proj_coord) <- c("X_proj","Y_proj") + subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50, proj_coord = proj_coord) + expect_equal(names(subplots) , c("subplot_id", "x", "y", "X_proj", "Y_proj")) + + # Test rectangular plot with rectangular division + rel_coord <- expand.grid(x = c(0, 100), y = c(0, 75)) + proj_coord <- as.matrix(rel_coord) %*% rot_mat + proj_coord <- sweep(proj_coord, 2, c(110,110), FUN = "+") + rect_subplots <- divide_plot(rel_coord, grid_size = c(50,25), proj_coord = proj_coord) + #ggplot(rect_subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x,y=y)) + geom_text(position = position_jitter()) + coord_equal() + expect_equivalent(rect_subplots[21:24,c(1,4,5)], + data.frame(subplot_id=rep("subplot_1_2",4), x_proj=c(128.3013,171.6025,115.8013,159.1025), y_proj=c(178.3013,203.3013,199.9519,224.9519)), + tol=1e-5) + + # Test when the origin is the NE corner : + inv_proj_coord <- proj_coord[c(4,3,2,1),] + rect_subplots <- divide_plot(rel_coord, grid_size = c(50,25), proj_coord = inv_proj_coord) + #ggplot(rect_subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x,y=y)) + geom_text(position = position_jitter()) + coord_equal() + expect_equivalent(rect_subplots[1:4,c(1,4,5)], + data.frame(subplot_id=rep("subplot_0_0",4), + x_proj=c(128.3013,171.6025,115.8013,159.1025)[c(4,3,2,1)], + y_proj=c(178.3013,203.3013,199.9519,224.9519)[c(4,3,2,1)]), + tol=1e-5) + + # Test multiple plots + rel_coord_1 <- expand.grid(x = c(0, 150), y = c(0, 100)) + proj_coord_1 <- as.matrix(rel_coord_1) %*% rot_mat + proj_coord_1 <- sweep(proj_coord_1, 2, c(250,110), FUN = "+") + rel_coord_2 <- expand.grid(x = c(0, 100), y = c(0, 100)) + proj_coord_2 <- as.matrix(rel_coord_2) %*% rot_mat + proj_coord_2 <- sweep(proj_coord_2, 2, c(60,110), FUN = "+") + proj_coord_2 <- proj_coord_2[c(4,3,2,1),] + + multiple_rel_coord <- rbind(rel_coord_1 , rel_coord_2) + multiple_proj_coord <- rbind(proj_coord_1,proj_coord_2) + multiple_rel_coord$plotID <- rep(c("plot1","plot2"),e=4) + multiple_subplots <- divide_plot(rel_coord = multiple_rel_coord[,c("x","y")], grid_size = 50, proj_coord = multiple_proj_coord, plot_ID_corner = multiple_rel_coord$plotID) + ggplot(multiple_subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x,y=y)) + geom_text(position = position_jitter()) + coord_equal() + expect_equivalent(multiple_subplots[c(24,40),] , + data.frame(plotID=c("plot1","plot2"),subplot_id=c("plot1_2_1","plot2_1_1"), + x=c(150,100), y=c(100,100), x_proj=c(329.9038,60), y_proj=c(271.6025,110)), + tol=1e-5) + +}) + +test_that("divide_plot with tree coordinates", { + + rel_coord_1 <- expand.grid(x = c(0, 100), y = c(0, 100)) + + # Test warning when a tree is not in any subplot + expect_warning(divide_plot(rel_coord = rel_coord_1, grid_size = 50, tree_coord = data.frame(x=110,y=110))) + + tree_coord_1 <- expand.grid(X = seq(0,100,10), Y = seq(0,100,10)) + subplots <- divide_plot(rel_coord = rel_coord_1, grid_size = 50, tree_coord = tree_coord_1) + expect_equal(subplots$tree_coord$subplot_id[c(1,6,56,61)],c("subplot_0_0","subplot_1_0","subplot_0_1","subplot_1_1")) + + # Test with multiple plots + rel_coord_2 <- expand.grid(x = c(0, 100), y = c(0, 100)) + multiple_rel_coord <- rbind(rel_coord_1 , rel_coord_2) + plot_ID_corner <- rep(c("plot1","plot2"),e=4) + tree_coord_2 <- expand.grid(X = seq(5,95,20), Y = seq(5,95,20)) + multiple_tree_coord <- rbind(tree_coord_1,tree_coord_2) + plot_ID_tree <- c(rep("plot1",nrow(tree_coord_1)),rep("plot3",nrow(tree_coord_2))) + + # Test error + expect_error( + divide_plot(rel_coord = multiple_rel_coord[,c("x","y")], grid_size = 50, plot_ID_corner = plot_ID_corner, + tree_coord = multiple_tree_coord, plot_ID_tree = plot_ID_tree) , + "Some plot_ID_tree are not found in plot_ID_corner") + + plot_ID_tree <- c(rep("plot1",nrow(tree_coord_1)),rep("plot2",nrow(tree_coord_2))) + + multiple_subplots <- divide_plot(rel_coord = multiple_rel_coord[,c("x","y")], grid_size = 50, + plot_ID_corner = plot_ID_corner, + tree_coord = multiple_tree_coord, plot_ID_tree = plot_ID_tree) + + expect_equal(multiple_subplots$tree_coord$subplot_id[c(122,125,137,140)], + c("plot2_0_0","plot2_1_0","plot2_0_1","plot2_1_1")) + +}) \ No newline at end of file From 9be34f1b7ed3ef9e826a6f56b748f9887d283b34 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 6 Nov 2024 16:14:36 +0100 Subject: [PATCH 19/25] - Improved documentations - Updated examples - Added divide_plot in the vignette - tree_coord now accept as many columns as wanted - Added tests for multiple tree_coord columns - Renamed vignette --- R/checkPlotCoord.R | 2 +- R/divide_plot.R | 72 ++++++++++++------- man/checkPlotCoord.Rd | 2 +- man/divide_plot.Rd | 57 +++++++++------ tests/testthat/test-dividePlot.R | 10 ++- ...lized_trees_and_forest_stands_metrics.Rmd} | 48 +++++++++++-- 6 files changed, 135 insertions(+), 56 deletions(-) rename vignettes/{VignetteSpatialRepresentationAndManagementPlot.Rmd => Vignette_spatialized_trees_and_forest_stands_metrics.Rmd} (78%) diff --git a/R/checkPlotCoord.R b/R/checkPlotCoord.R index 0774350..9dd5895 100644 --- a/R/checkPlotCoord.R +++ b/R/checkPlotCoord.R @@ -4,7 +4,7 @@ #' Quality check of plot corner and tree coordinates. #' #' @details -#' If trustGPScorners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the maxDist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. +#' If trustGPScorners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the maxDist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. Be aware that this projection only works if the plot, in the relative coordinates system, is rectangular (ie, has 4 right angles). #' #' If trustGPScorners is FALSE, corner coordinates in the projected coordinate system are calculated by a procrust analysis that preserves the shape and dimensions of the plot in the local coordinate system. Outlier corners are also identified sequentially and projected coordinates of the trees are calculated by applying the resulting procrust analysis. #' diff --git a/R/divide_plot.R b/R/divide_plot.R index fd78d65..df11fc7 100644 --- a/R/divide_plot.R +++ b/R/divide_plot.R @@ -1,39 +1,53 @@ #' Divides one ore more plots into subplots #' -#' This function divides a plot (or several plots) into subplots in the relative coordinates system, and returns the coordinates of subplot corners -#' If corner coordinates in the projected coordinate system is supplied (proj_coord), subplot projected coordinates are calculated by a bilinear interpolation in relation with plot corner relative coordinates. +#' @description +#' This function divides a plot (or several plots) into subplots in the relative coordinates system, and returns the coordinates of subplot corners. #' -#' @param rel_coord a data frame containing the relative (local) coordinates of plot corners, with X and Y on the first and second column respectively -#' @param proj_coord (optional) a data frame containing the projected coordinates of plot corners, with X and Y on the first and second column respectively, and with the same row order than rel_coord -#' @param grid_size a vector indicating the dimensions of grid cells (resp. X and Y dimensions). If only one value is given, the grid cells will be considered as squares. +#' @details +#' If corner coordinates in the projected coordinate system are supplied (proj_coord), projected coordinates of subplot corners are calculated by a bilinear interpolation in relation with relative coordinates of plot corners. Be aware that this bilinear interpolation only works if the plot in the relative coordinates system is rectangular (ie, has 4 right angles). +#' +#' @param rel_coord a data frame containing the relative (local) coordinates of plot corners, with X and Y corresponding to the first and second column respectively +#' @param proj_coord a data frame containing the projected coordinates of plot corners, with X and Y corresponding to the first and second column respectively, and with the same row order than rel_coord +#' @param grid_size a vector indicating the dimensions of grid cells (resp. X and Y dimensions). If only one value is given, grid cells will be considered as squares. #' @param plot_ID_corner if dealing with multiple plots : a vector indicating plot IDs for corners. -#' @param tree_coord (otpional) a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively +#' @param tree_coord a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively #' @param plot_ID_tree if dealing with multiple plots : a vector indicating tree plot IDs. #' -#' @return Returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : +#' @return If tree_coord isn't supplied, returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : #' - `plot_ID_corner`: If dealing with multiple plots : the plot code -#' - `subplot`: The automatically generated subplot code -#' - `Xrel`: The relative X-axis coordinates of subplots corners -#' - `Yrel`: The relative Y-axis coordinates of subplots corners -#' - `X`: If proj_coord is supplied : the projected X-axis coordinates of subplots corners -#' - `Y`: If proj_coord is supplied : the projected Y-axis coordinates of subplots corners +#' - `subplot_id`: The automatically generated subplot code, using the following rule : subplot_X_Y +#' - `x` and `y` (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. +#' - `x_proj` and `y_proj` (or the column names provided by proj_coord) : if proj_coord is supplied, the projected X-axis and Y-axis coordinates of subplots corners +#' +#' If tree_coord is supplied, returns a list containing the previous data-frame and a data-frame containing for each tree : +#' - the relative coordinates +#' - the subplot_id associated +#' - the rest of tree_coord columns supplied #' #' @export #' @author Arthur PERE, Arthur BAILLY #' @importFrom data.table data.table := setnames setcolorder #' @examples #' -#' coord <- data.frame(X = c(0, 200, 0, 200), Y = c(0, 0, 200, 200)) + 5000 -#' cornerNum <- c(1, 2, 4, 3) -#' plot <- rep("plot1", 4) +#' # Rectangular plot and grid cells +#' rel_coord <- data.frame(x_rel = c(0, 200, 0, 200), y_rel = c(0, 0, 100, 100)) +#' subplots <- divide_plot(rel_coord, grid_size = c(100,50)) +#' +#' # Squared plot +#' rel_coord <- data.frame(x_rel = c(0, 200, 0, 200), y_rel = c(0, 0, 200, 200)) +#' proj_coord <- data.frame(x_proj = c(110, 190, 60, 145), y_proj = c(110, 160, 196, 245)) +#' subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100) #' -#' cut <- cutPlot(coord, plot, cornerNum, grid_size = 100, dimX = 200, dimY = 200) +#' tree_coord <- data.frame(x_proj = runif(50,0,200), y_proj = runif(50,0,200)) +#' subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100, tree_coord = tree_coord) #' -#' # plot the result -#' plot(coord, main = "example", xlim = c(4900, 5300), ylim = c(4900, 5300), asp = 1) -#' text(coord, labels = cornerNum, pos = 1) -#' points(cut$XAbs, cut$YAbs, pch = "+") -#' legend("bottomright", legend = c("orignal", "cut"), pch = c("o", "+")) +#' # Dealing with multiple plots +#' rel_coord <- rbind(rel_coord, rel_coord) +#' proj_coord <- rbind(proj_coord, proj_coord + 200) +#' tree_coord <- rbind(tree_coord, data.frame(x_proj = runif(50,0,200), y_proj = runif(50,0,200))) +#' plot_ID_corner <- rep(c("plot1","plot2"), e=4) +#' plot_ID_tree <- rep(c("plot1","plot2"), e=50) +#' subplots <- divide_plot(rel_coord, proj_coord, 100, plot_ID_corner, tree_coord, plot_ID_tree) #' divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner = NULL, tree_coord = NULL, plot_ID_tree = NULL) { @@ -151,13 +165,15 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner if(!is.null(tree_coord)) { tree_coord <- as.data.table(tree_coord) if(!is.null(plot_ID_tree)){ - tree_coord[,plot_ID_tree:=plot_ID_tree] + tree_coord[,plot_id:=plot_ID_tree] + setcolorder(tree_coord , "plot_id" , after = 2) # if tree_coord has more than 2 col } else{ - tree_coord[,plot_ID_tree:=""] + tree_coord[,plot_id:=""] + setcolorder(tree_coord , "plot_id" , after = 2) # if tree_coord has more than 2 col } invisible(lapply(split(sub_corner_coord, by = "subplot_id", keep.by = TRUE), function(dat) { - tree_coord[ plot_ID_tree == dat$plot_ID_corner[1] & + tree_coord[ plot_id == dat$plot_ID_corner[1] & get(names(tree_coord)[1]) %between% range(dat[[3]]) & get(names(tree_coord)[2]) %between% range(dat[[4]]), subplot_id := dat$subplot_id[1] ] @@ -167,7 +183,10 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner warning("One or more trees could not be assigned to a subplot (not in a subplot area)") } if(is.null(plot_ID_tree)) { - tree_coord[ , c("subplot_id","plot_ID_tree") := list(paste0("subplot",subplot_id),NULL)] + tree_coord[ , c("subplot_id","plot_id") := list(paste0("subplot",subplot_id),NULL)] + setcolorder(tree_coord , "subplot_id" , after = 2) + } else { + setcolorder(tree_coord , "subplot_id" , after = 3) } } @@ -178,12 +197,11 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner sub_corner_coord[ , c("subplot_id","plot_ID_corner") := list(paste0("subplot",subplot_id),NULL)] } - if(is.null(tree_coord)) { output <- data.frame(sub_corner_coord) } else { output <- list(sub_corner_coord = data.frame(sub_corner_coord), - tree_coord = tree_coord) + tree_coord = data.frame(tree_coord)) } return(output) diff --git a/man/checkPlotCoord.Rd b/man/checkPlotCoord.Rd index aee5517..0b5a520 100644 --- a/man/checkPlotCoord.Rd +++ b/man/checkPlotCoord.Rd @@ -50,7 +50,7 @@ Returns a list including : Quality check of plot corner and tree coordinates. } \details{ -If trustGPScorners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the maxDist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. +If trustGPScorners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the maxDist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. Be aware that this projection only works if the plot, in the relative coordinates system, is rectangular (ie, has 4 right angles). If trustGPScorners is FALSE, corner coordinates in the projected coordinate system are calculated by a procrust analysis that preserves the shape and dimensions of the plot in the local coordinate system. Outlier corners are also identified sequentially and projected coordinates of the trees are calculated by applying the resulting procrust analysis. diff --git a/man/divide_plot.Rd b/man/divide_plot.Rd index b3300ff..56c072a 100644 --- a/man/divide_plot.Rd +++ b/man/divide_plot.Rd @@ -14,46 +14,61 @@ divide_plot( ) } \arguments{ -\item{rel_coord}{a data frame containing the relative (local) coordinates of plot corners, with X and Y on the first and second column respectively} +\item{rel_coord}{a data frame containing the relative (local) coordinates of plot corners, with X and Y corresponding to the first and second column respectively} -\item{proj_coord}{(optional) a data frame containing the projected coordinates of plot corners, with X and Y on the first and second column respectively, and with the same row order than rel_coord} +\item{proj_coord}{a data frame containing the projected coordinates of plot corners, with X and Y corresponding to the first and second column respectively, and with the same row order than rel_coord} -\item{grid_size}{a vector indicating the dimensions of grid cells (resp. X and Y dimensions). If only one value is given, the grid cells will be considered as squares.} +\item{grid_size}{a vector indicating the dimensions of grid cells (resp. X and Y dimensions). If only one value is given, grid cells will be considered as squares.} \item{plot_ID_corner}{if dealing with multiple plots : a vector indicating plot IDs for corners.} -\item{tree_coord}{(otpional) a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively} +\item{tree_coord}{a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively} \item{plot_ID_tree}{if dealing with multiple plots : a vector indicating tree plot IDs.} } \value{ -Returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : +If tree_coord isn't supplied, returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : \itemize{ \item \code{plot_ID_corner}: If dealing with multiple plots : the plot code -\item \code{subplot}: The automatically generated subplot code -\item \code{Xrel}: The relative X-axis coordinates of subplots corners -\item \code{Yrel}: The relative Y-axis coordinates of subplots corners -\item \code{X}: If proj_coord is supplied : the projected X-axis coordinates of subplots corners -\item \code{Y}: If proj_coord is supplied : the projected Y-axis coordinates of subplots corners +\item \code{subplot_id}: The automatically generated subplot code, using the following rule : subplot_X_Y +\item \code{x} and \code{y} (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. +\item \code{x_proj} and \code{y_proj} (or the column names provided by proj_coord) : if proj_coord is supplied, the projected X-axis and Y-axis coordinates of subplots corners +} + +If tree_coord is supplied, returns a list containing the previous data-frame and a data-frame containing for each tree : +\itemize{ +\item the relative coordinates +\item the subplot_id associated +\item the rest of tree_coord columns supplied } } \description{ -This function divides a plot (or several plots) into subplots in the relative coordinates system, and returns the coordinates of subplot corners -If corner coordinates in the projected coordinate system is supplied (proj_coord), subplot projected coordinates are calculated by a bilinear interpolation in relation with plot corner relative coordinates. +This function divides a plot (or several plots) into subplots in the relative coordinates system, and returns the coordinates of subplot corners. +} +\details{ +If corner coordinates in the projected coordinate system are supplied (proj_coord), projected coordinates of subplot corners are calculated by a bilinear interpolation in relation with relative coordinates of plot corners. Be aware that this bilinear interpolation only works if the plot in the relative coordinates system is rectangular (ie, has 4 right angles). } \examples{ -coord <- data.frame(X = c(0, 200, 0, 200), Y = c(0, 0, 200, 200)) + 5000 -cornerNum <- c(1, 2, 4, 3) -plot <- rep("plot1", 4) +# Rectangular plot and grid cells +rel_coord <- data.frame(x_rel = c(0, 200, 0, 200), y_rel = c(0, 0, 100, 100)) +subplots <- divide_plot(rel_coord, grid_size = c(100,50)) + +# Squared plot +rel_coord <- data.frame(x_rel = c(0, 200, 0, 200), y_rel = c(0, 0, 200, 200)) +proj_coord <- data.frame(x_proj = c(110, 190, 60, 145), y_proj = c(110, 160, 196, 245)) +subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100) -cut <- cutPlot(coord, plot, cornerNum, grid_size = 100, dimX = 200, dimY = 200) +tree_coord <- data.frame(x_proj = runif(50,0,200), y_proj = runif(50,0,200)) +subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100, tree_coord = tree_coord) -# plot the result -plot(coord, main = "example", xlim = c(4900, 5300), ylim = c(4900, 5300), asp = 1) -text(coord, labels = cornerNum, pos = 1) -points(cut$XAbs, cut$YAbs, pch = "+") -legend("bottomright", legend = c("orignal", "cut"), pch = c("o", "+")) +# Dealing with multiple plots +rel_coord <- rbind(rel_coord, rel_coord) +proj_coord <- rbind(proj_coord, proj_coord + 200) +tree_coord <- rbind(tree_coord, data.frame(x_proj = runif(50,0,200), y_proj = runif(50,0,200))) +plot_ID_corner <- rep(c("plot1","plot2"), e=4) +plot_ID_tree <- rep(c("plot1","plot2"), e=50) +subplots <- divide_plot(rel_coord, proj_coord, 100, plot_ID_corner, tree_coord, plot_ID_tree) } \author{ diff --git a/tests/testthat/test-dividePlot.R b/tests/testthat/test-dividePlot.R index a61a691..538aa0d 100644 --- a/tests/testthat/test-dividePlot.R +++ b/tests/testthat/test-dividePlot.R @@ -132,7 +132,7 @@ test_that("divide_plot with projected coordinates", { multiple_proj_coord <- rbind(proj_coord_1,proj_coord_2) multiple_rel_coord$plotID <- rep(c("plot1","plot2"),e=4) multiple_subplots <- divide_plot(rel_coord = multiple_rel_coord[,c("x","y")], grid_size = 50, proj_coord = multiple_proj_coord, plot_ID_corner = multiple_rel_coord$plotID) - ggplot(multiple_subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x,y=y)) + geom_text(position = position_jitter()) + coord_equal() + #ggplot(multiple_subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x,y=y)) + geom_text(position = position_jitter()) + coord_equal() expect_equivalent(multiple_subplots[c(24,40),] , data.frame(plotID=c("plot1","plot2"),subplot_id=c("plot1_2_1","plot2_1_1"), x=c(150,100), y=c(100,100), x_proj=c(329.9038,60), y_proj=c(271.6025,110)), @@ -151,11 +151,17 @@ test_that("divide_plot with tree coordinates", { subplots <- divide_plot(rel_coord = rel_coord_1, grid_size = 50, tree_coord = tree_coord_1) expect_equal(subplots$tree_coord$subplot_id[c(1,6,56,61)],c("subplot_0_0","subplot_1_0","subplot_0_1","subplot_1_1")) + # Test with multiple columns for tree_coord + tree_coord_1[c("H","AGB")] <- data.frame(H = rnorm(nrow(tree_coord_1),10,3), AGB = runif(nrow(tree_coord_1),0,2)) + subplots <- divide_plot(rel_coord = rel_coord_1, grid_size = 50, tree_coord = tree_coord_1) + expect_equal(names(subplots$tree_coord) , c("X","Y","subplot_id","H","AGB")) + # Test with multiple plots rel_coord_2 <- expand.grid(x = c(0, 100), y = c(0, 100)) multiple_rel_coord <- rbind(rel_coord_1 , rel_coord_2) plot_ID_corner <- rep(c("plot1","plot2"),e=4) tree_coord_2 <- expand.grid(X = seq(5,95,20), Y = seq(5,95,20)) + tree_coord_2[c("H","AGB")] <- data.frame(H = rnorm(nrow(tree_coord_2),10,3), AGB = runif(nrow(tree_coord_2),0,2)) multiple_tree_coord <- rbind(tree_coord_1,tree_coord_2) plot_ID_tree <- c(rep("plot1",nrow(tree_coord_1)),rep("plot3",nrow(tree_coord_2))) @@ -173,5 +179,7 @@ test_that("divide_plot with tree coordinates", { expect_equal(multiple_subplots$tree_coord$subplot_id[c(122,125,137,140)], c("plot2_0_0","plot2_1_0","plot2_0_1","plot2_1_1")) + expect_equal(names(multiple_subplots$tree_coord) , c("X","Y","plot_id","subplot_id","H","AGB")) + }) \ No newline at end of file diff --git a/vignettes/VignetteSpatialRepresentationAndManagementPlot.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd similarity index 78% rename from vignettes/VignetteSpatialRepresentationAndManagementPlot.Rmd rename to vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index 877cc46..2fe0178 100644 --- a/vignettes/VignetteSpatialRepresentationAndManagementPlot.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -20,7 +20,8 @@ knitr::opts_chunk$set( collapse = TRUE, echo = TRUE, comment = "#>", fig.align = "center" ) -require(BIOMASS) +#require(BIOMASS) +devtools::load_all() require(knitr) require(ggplot2) require(data.table) @@ -96,7 +97,7 @@ If **enough coordinates have been recorded for each corner** (for more informati If only **4 GPS measurements** have been taken **with a high degree of accuracy** (by a geometer, for example), or if you have averaged your measurements by yourself, you can supply these 4 GPS coordinates to the function. ```{r dpi=200} -check_plot <- checkPlotCoord( +check_plot_trust_GPS <- checkPlotCoord( longlat = cornerCoord[c("Long", "Lat")], # OR, if exists : projCoord = cornerCoord[c("xProj", "yProj")], relCoord = cornerCoord[c("xRel", "yRel")], @@ -115,7 +116,7 @@ Let's degrade the data to simulate the fact that we only have 8 GPS coordinates ```{r dpi=200} degradedCornerCoord <- cornerCoord[c(1:2,11:12,21:22,31:32),] -check_plot <- checkPlotCoord( +check_plot_trust_field <- checkPlotCoord( longlat = degradedCornerCoord[, c("Long", "Lat")], # OR projCoord = cornerCoord[c("xProj", "yProj")], relCoord = degradedCornerCoord[, c("xRel", "yRel")], @@ -155,7 +156,7 @@ The tree coordinates can be obtained via the `$treeProjCoord` output. You can also grep and modify the plot via the `$plotDesign` output which is a ggplot object. For exemple, to change the plot title : ```{r} -plot_to_change <- check_plot$plotDesign +plot_to_change <- check_plot_trust_GPS$plotDesign plot_to_change$labels$title <- "A nice title !" plot_to_change ``` @@ -173,6 +174,43 @@ WORK IN PROGRESS # Dividing plots -WORK IN PROGRESS +Dividing plots into several rectangular sub-plots is performed using the `divide_plot()` function. This function takes the relative coordinates of the plot corners for division (be aware that the plot must be rectangular, ie, have 4 right angles) : +```{r divide_plot} +subplots <- divide_plot( + rel_coord = check_plot_trust_field$cornerCoord[c("Xrel","Yrel")], + proj_coord = check_plot_trust_field$cornerCoord[c("X","Y")], #optional + grid_size = 50 # or c(50,50) + ) +``` +For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners (keeping column names) but can also assign to each tree its subplot with the `tree_coord` : + +```{r divide_plot_trees} +trees <- trees %>% relocate(xRel, yRel, .before = everything()) # to get the coordinates at the first two columns + +subplots <- divide_plot( + rel_coord = check_plot_trust_field$cornerCoord[c("Xrel","Yrel")], + proj_coord = check_plot_trust_field$cornerCoord[c("X","Y")], #optional + grid_size = 50, # or c(50,50) + tree_coord = trees + ) +``` + +Note that `tree_coord` can contain any number of columns, but the first two must be the coordinates of the trees in the relative system. + +The function can handle as many plots as you wish, using the `plot_ID_corner` and `plot_ID_tree` arguments : + +```{r divide_multiple_plots} +multiple_plots <- rbind(check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")], check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")]) +multiple_plots$plot_ID = rep(c("NB1","NB2"),e=4) +multiple_trees <- rbind(trees,trees %>% mutate(plot="NB2")) + +mulitple_subplots <- divide_plot( + rel_coord = multiple_plots[c("Xrel","Yrel")], + grid_size = 50, + plot_ID_corner = multiple_plots$plot_ID, + tree_coord = multiple_trees, + plot_ID_tree = multiple_trees$plot +) +``` From 8c7d06c8b9c61188b4c049264325cfbf812e67e7 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 6 Nov 2024 16:27:06 +0100 Subject: [PATCH 20/25] Fixed forgotten tests and decomment devtools::load_all() --- tests/testthat/test-bilinearInterpolation.R | 6 +++--- ...ignette_spatialized_trees_and_forest_stands_metrics.Rmd | 7 +++---- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/tests/testthat/test-bilinearInterpolation.R b/tests/testthat/test-bilinearInterpolation.R index 1517407..bbfdda7 100644 --- a/tests/testthat/test-bilinearInterpolation.R +++ b/tests/testthat/test-bilinearInterpolation.R @@ -30,8 +30,8 @@ test_that("bilinear_interpolation function", { expect_is(proj_coord,"data.frame") expect_equal(nrow(proj_coord), nrow(coord)) - expect_equal(names(proj_coord), c("Xinterp","Yinterp")) - expect_equal(proj_coord, data.frame(Xinterp=c(59.15064,81.96152),Yinterp=c(134.1506,164.6410)) , tolerance = 1e-6) + expect_equal(names(proj_coord), c("x_interp","y_interp")) + expect_equal(proj_coord, data.frame(x_interp=c(59.15064,81.96152),y_interp=c(134.1506,164.6410)) , tolerance = 1e-6) # Test when to_corner_coord colnames are supplied colnames(to_corner_coord) <- c("Xproj","Yproj") @@ -69,7 +69,7 @@ test_that("bilinear_interpolation function", { #coord <- expand.grid(X = seq(0,100,10), Y = seq(0,50,5)) proj_coord = bilinear_interpolation(coord = coord, from_corner_coord = from_corner_coord, to_corner_coord = to_corner_coord) #plot(from_corner_coord , xlim=c(-200,110), ylim=c(-110,60) , asp=1) ; points(coord , col="blue") ; points(proj_coord , col="purple") ; points(to_corner_coord, col="red") - expect_equal(proj_coord, data.frame(Xinterp=c(-140.8494 ,-118.0385 ),Yinterp=c(-65.84936,-35.35898)) , tolerance = 1e-6) + expect_equal(proj_coord, data.frame(x_interp=c(-140.8494 ,-118.0385 ),y_interp=c(-65.84936,-35.35898)) , tolerance = 1e-6) # If the origin is the NE corner : to_corner_coord <- sweep(from_corner_coord, 2, c(60,60), FUN = "+") diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index 2fe0178..d2e1cdf 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -1,5 +1,5 @@ --- -title: "Manage plot and tree coordinates with BIOMASS" +title: "Spatialized trees and forest stands metrics with BIOMASS" author: "Arthur Bailly" date: "`r Sys.Date()`" output: @@ -10,7 +10,7 @@ output: self_contained: yes theme: cayman vignette: > - %\VignetteIndexEntry{Manage plot and tree coordinates with BIOMASS} + %\VignetteIndexEntry{Spatialized trees and forest stands metrics with BIOMASS} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -20,8 +20,7 @@ knitr::opts_chunk$set( collapse = TRUE, echo = TRUE, comment = "#>", fig.align = "center" ) -#require(BIOMASS) -devtools::load_all() +require(BIOMASS) require(knitr) require(ggplot2) require(data.table) From 89eeb9b2b2c25ab4a804542349266a6000d57ba9 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 6 Nov 2024 16:33:35 +0100 Subject: [PATCH 21/25] fixed vignette error --- .../Vignette_spatialized_trees_and_forest_stands_metrics.Rmd | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index d2e1cdf..fd51579 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -25,6 +25,7 @@ require(knitr) require(ggplot2) require(data.table) require(tidyverse) +require(dplyr) ``` # Overview @@ -186,7 +187,7 @@ subplots <- divide_plot( For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners (keeping column names) but can also assign to each tree its subplot with the `tree_coord` : ```{r divide_plot_trees} -trees <- trees %>% relocate(xRel, yRel, .before = everything()) # to get the coordinates at the first two columns +trees <- relocate(trees, xRel, yRel, .before = everything()) # to get the coordinates at the first two columns subplots <- divide_plot( rel_coord = check_plot_trust_field$cornerCoord[c("Xrel","Yrel")], @@ -203,7 +204,7 @@ The function can handle as many plots as you wish, using the `plot_ID_corner` an ```{r divide_multiple_plots} multiple_plots <- rbind(check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")], check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")]) multiple_plots$plot_ID = rep(c("NB1","NB2"),e=4) -multiple_trees <- rbind(trees,trees %>% mutate(plot="NB2")) +multiple_trees <- rbind(trees, mutate(trees, plot="NB2")) mulitple_subplots <- divide_plot( rel_coord = multiple_plots[c("Xrel","Yrel")], From 234acece5336e0f9f4ef5e40f71998a276d81849 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 6 Nov 2024 16:39:00 +0100 Subject: [PATCH 22/25] Fixed vignette error with the use of relocate from dplyr ?? --- .../Vignette_spatialized_trees_and_forest_stands_metrics.Rmd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index fd51579..bd6318e 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -187,7 +187,7 @@ subplots <- divide_plot( For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners (keeping column names) but can also assign to each tree its subplot with the `tree_coord` : ```{r divide_plot_trees} -trees <- relocate(trees, xRel, yRel, .before = everything()) # to get the coordinates at the first two columns +trees <- dplyr::relocate(trees, xRel, yRel, .before = everything()) # to get the coordinates at the first two columns subplots <- divide_plot( rel_coord = check_plot_trust_field$cornerCoord[c("Xrel","Yrel")], From 8030720a391a181aa8c976c76f6da1ef20ae4bcb Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 6 Nov 2024 16:39:00 +0100 Subject: [PATCH 23/25] Fixed vignette error with the use of relocate from dplyr ?? --- ...Vignette_spatialized_trees_and_forest_stands_metrics.Rmd | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index fd51579..275afb0 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -187,7 +187,7 @@ subplots <- divide_plot( For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners (keeping column names) but can also assign to each tree its subplot with the `tree_coord` : ```{r divide_plot_trees} -trees <- relocate(trees, xRel, yRel, .before = everything()) # to get the coordinates at the first two columns +trees <- trees[,c("xRel","yRel","plot","D","WD","H")] subplots <- divide_plot( rel_coord = check_plot_trust_field$cornerCoord[c("Xrel","Yrel")], @@ -204,7 +204,9 @@ The function can handle as many plots as you wish, using the `plot_ID_corner` an ```{r divide_multiple_plots} multiple_plots <- rbind(check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")], check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")]) multiple_plots$plot_ID = rep(c("NB1","NB2"),e=4) -multiple_trees <- rbind(trees, mutate(trees, plot="NB2")) +trees_2 <- trees +trees_2$plot <- "NB2" +multiple_trees <- rbind(trees, trees_2) mulitple_subplots <- divide_plot( rel_coord = multiple_plots[c("Xrel","Yrel")], From 90b24e5453d350bcf57c8ad5817b1b300cf8a57b Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 7 Nov 2024 10:51:29 +0100 Subject: [PATCH 24/25] Added dplyr in suggests --- DESCRIPTION | 3 ++- R/divide_plot.R | 4 ++-- tests/testthat/test-dividePlot.R | 4 ++-- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 00eebfc..85ab2db 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -51,5 +51,6 @@ Suggests: curl, geodata, httr2, - pkgdown + pkgdown, + dplyr RoxygenNote: 7.3.2 diff --git a/R/divide_plot.R b/R/divide_plot.R index df11fc7..bba2d8f 100644 --- a/R/divide_plot.R +++ b/R/divide_plot.R @@ -16,7 +16,7 @@ #' @return If tree_coord isn't supplied, returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : #' - `plot_ID_corner`: If dealing with multiple plots : the plot code #' - `subplot_id`: The automatically generated subplot code, using the following rule : subplot_X_Y -#' - `x` and `y` (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. +#' - `x_rel` and `y_rel` (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. #' - `x_proj` and `y_proj` (or the column names provided by proj_coord) : if proj_coord is supplied, the projected X-axis and Y-axis coordinates of subplots corners #' #' If tree_coord is supplied, returns a list containing the previous data-frame and a data-frame containing for each tree : @@ -53,7 +53,7 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner # Parameters verification ---------------------------------------------------- if (is.matrix(rel_coord)) { - rel_coord <- data.frame(x=rel_coord[,1],y=rel_coord[,2]) + rel_coord <- data.frame(x_rel=rel_coord[,1],y_rel=rel_coord[,2]) } if (!is.null(proj_coord) && !is.data.frame(proj_coord)) { proj_coord <- data.frame(x_proj=proj_coord[,1], y_proj=proj_coord[,2]) diff --git a/tests/testthat/test-dividePlot.R b/tests/testthat/test-dividePlot.R index 538aa0d..d2b0b09 100644 --- a/tests/testthat/test-dividePlot.R +++ b/tests/testthat/test-dividePlot.R @@ -33,7 +33,7 @@ test_that("divide_plot on relative coordinates only", { # Test when rel_coord is a matrix (colnames aren't supplied) rel_coord <- matrix(c(0,100,0,100,0,0,100,100), ncol=2) subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50) - expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x=c(0,50,0,50), y=c(0,0,50,50))) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x_rel=c(0,50,0,50), y_rel=c(0,0,50,50))) # Test when rel_coord is a data.table rel_coord <- data.table(expand.grid(x_field = c(0, 100), y_field = c(0, 100))) @@ -182,4 +182,4 @@ test_that("divide_plot with tree coordinates", { expect_equal(names(multiple_subplots$tree_coord) , c("X","Y","plot_id","subplot_id","H","AGB")) -}) \ No newline at end of file +}) From f24f51f010167da52370ad8b32b39df6c9bed1ef Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 7 Nov 2024 10:51:29 +0100 Subject: [PATCH 25/25] Added dplyr in suggests Re-used dplyr functions and pipe in vignette --- DESCRIPTION | 3 ++- R/divide_plot.R | 4 ++-- tests/testthat/test-dividePlot.R | 4 ++-- ...Vignette_spatialized_trees_and_forest_stands_metrics.Rmd | 6 ++---- 4 files changed, 8 insertions(+), 9 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 00eebfc..85ab2db 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -51,5 +51,6 @@ Suggests: curl, geodata, httr2, - pkgdown + pkgdown, + dplyr RoxygenNote: 7.3.2 diff --git a/R/divide_plot.R b/R/divide_plot.R index df11fc7..bba2d8f 100644 --- a/R/divide_plot.R +++ b/R/divide_plot.R @@ -16,7 +16,7 @@ #' @return If tree_coord isn't supplied, returns a data-frame containing as many rows as there are corners corresponding to the subplots, and the following columns : #' - `plot_ID_corner`: If dealing with multiple plots : the plot code #' - `subplot_id`: The automatically generated subplot code, using the following rule : subplot_X_Y -#' - `x` and `y` (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. +#' - `x_rel` and `y_rel` (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. #' - `x_proj` and `y_proj` (or the column names provided by proj_coord) : if proj_coord is supplied, the projected X-axis and Y-axis coordinates of subplots corners #' #' If tree_coord is supplied, returns a list containing the previous data-frame and a data-frame containing for each tree : @@ -53,7 +53,7 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner # Parameters verification ---------------------------------------------------- if (is.matrix(rel_coord)) { - rel_coord <- data.frame(x=rel_coord[,1],y=rel_coord[,2]) + rel_coord <- data.frame(x_rel=rel_coord[,1],y_rel=rel_coord[,2]) } if (!is.null(proj_coord) && !is.data.frame(proj_coord)) { proj_coord <- data.frame(x_proj=proj_coord[,1], y_proj=proj_coord[,2]) diff --git a/tests/testthat/test-dividePlot.R b/tests/testthat/test-dividePlot.R index 538aa0d..d2b0b09 100644 --- a/tests/testthat/test-dividePlot.R +++ b/tests/testthat/test-dividePlot.R @@ -33,7 +33,7 @@ test_that("divide_plot on relative coordinates only", { # Test when rel_coord is a matrix (colnames aren't supplied) rel_coord <- matrix(c(0,100,0,100,0,0,100,100), ncol=2) subplots <- divide_plot(rel_coord = rel_coord, grid_size = 50) - expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x=c(0,50,0,50), y=c(0,0,50,50))) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x_rel=c(0,50,0,50), y_rel=c(0,0,50,50))) # Test when rel_coord is a data.table rel_coord <- data.table(expand.grid(x_field = c(0, 100), y_field = c(0, 100))) @@ -182,4 +182,4 @@ test_that("divide_plot with tree coordinates", { expect_equal(names(multiple_subplots$tree_coord) , c("X","Y","plot_id","subplot_id","H","AGB")) -}) \ No newline at end of file +}) diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index 275afb0..7d45177 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -187,7 +187,7 @@ subplots <- divide_plot( For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners (keeping column names) but can also assign to each tree its subplot with the `tree_coord` : ```{r divide_plot_trees} -trees <- trees[,c("xRel","yRel","plot","D","WD","H")] +trees <- trees %>% dplyr::relocate(xRel, yRel, .before = everything()) # to get the coordinates at the first two columns subplots <- divide_plot( rel_coord = check_plot_trust_field$cornerCoord[c("Xrel","Yrel")], @@ -204,9 +204,7 @@ The function can handle as many plots as you wish, using the `plot_ID_corner` an ```{r divide_multiple_plots} multiple_plots <- rbind(check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")], check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")]) multiple_plots$plot_ID = rep(c("NB1","NB2"),e=4) -trees_2 <- trees -trees_2$plot <- "NB2" -multiple_trees <- rbind(trees, trees_2) +multiple_trees <- rbind(trees, trees %>% dplyr::mutate(plot="NB2")) mulitple_subplots <- divide_plot( rel_coord = multiple_plots[c("Xrel","Yrel")],