From 97046593abc1bfbadd356d678716452bf0473f3a Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 13 Nov 2024 15:08:59 +0100 Subject: [PATCH 01/16] In checkPlotCoord function : replaced returned coordinates column names from (Xrel;Yrel) to (x_rel,y_rel) and (X;Y)to (x_proj;y_proj) --- R/checkPlotCoord.R | 68 ++++++++++++++-------------- R/divide_plot.R | 65 +++++++++++++------------- tests/testthat/test-checkCoordPlot.R | 12 ++--- tests/testthat/test-dividePlot.R | 52 +++++++++++---------- 4 files changed, 103 insertions(+), 94 deletions(-) diff --git a/R/checkPlotCoord.R b/R/checkPlotCoord.R index 9dd5895..fe56760 100644 --- a/R/checkPlotCoord.R +++ b/R/checkPlotCoord.R @@ -23,12 +23,12 @@ #' @author Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY #' #' @return Returns a list including : -#' - `cornerCoord`: a data frame containing the projected and relative coordinates, the ID (if cornerID is supplied) and the number of the 4 corners of the plot +#' - `cornerCoord`: a data frame containing the projected coordinates (x_proj and y_proj), the relative coordinates (x_rel and y_rel), and the ID (if cornerID is supplied) of the 4 corners of the plot #' - `polygon`: a spatial polygon #' - `outliers`: a data frame containing the projected coordinates, the ID (if cornerID is supplied) and the row number of GPS points considered outliers #' - `plotDesign`: if `drawPlot` is TRUE, a ggplot object corresponding to the design of the plot #' - `codeUTM`: if `longlat` is supplied, a character containing the UTM code of the GPS coordinates -#' - `treeProjCoord`: if `treeCoord` is supplied, a data frame containing the coordinates of the trees in the projected coordinate system +#' - `treeProjCoord`: if `treeCoord` is supplied, a data frame containing the coordinates of the trees in the projected coordinate system (x_proj and y_proj) #' #' @export #' @@ -147,13 +147,13 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc 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")) + setnames(cornerCoord, c("x_proj","y_proj","x_rel","y_rel","cornerID")) if (any(table(cornerCoord$cornerID) < 5)) { warning("At least one corner has less than 5 measurements. We suggest using the argument trustGPScorners = FALSE") } - cornerCoord[ , c("Xmean","Ymean","nRow") := list(mean(X),mean(Y),.I) , by=cornerID] - cornerCoord[ , outlier := ifelse( sqrt((X - Xmean)^2 + (Y-Ymean)^2) > maxDist, T, F)] - outliers <- cornerCoord[outlier==T , c("X","Y","cornerID","nRow")] + cornerCoord[ , c("Xmean","Ymean","nRow") := list(mean(x_proj),mean(y_proj),.I) , by=cornerID] + cornerCoord[ , outlier := ifelse( sqrt((x_proj - Xmean)^2 + (y_proj-Ymean)^2) > maxDist, T, F)] + outliers <- cornerCoord[outlier==T , c("x_proj","y_proj","cornerID","nRow")] if (any(cornerCoord[,sum(!outlier),by=cornerID][[2]]==0)) { stop("All measurements for at least one corner are considered as outliers.\n @@ -165,21 +165,21 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc cornerCoord <- cornerCoord[outlier==F ,] existOutliers <- T while(existOutliers) { - cornerCoord[ , outlier := ifelse( sqrt((X - Xmean)^2 + (Y-Ymean)^2) > maxDist, T, F)] - outliers <- rbind(outliers , cornerCoord[outlier==T , c("X","Y","cornerID","nRow")]) + cornerCoord[ , outlier := ifelse( sqrt((x_proj - Xmean)^2 + (y_proj-Ymean)^2) > maxDist, T, F)] + outliers <- rbind(outliers , cornerCoord[outlier==T , c("x_proj","y_proj","cornerID","nRow")]) cornerCoord <- cornerCoord[outlier==F ,] if(nrow(cornerCoord)+nrow(outliers) == nrow(projCoord)) existOutliers <- F - cornerCoord[ , c("Xmean","Ymean") := list(mean(X),mean(Y)) , by=cornerID] + cornerCoord[ , c("Xmean","Ymean") := list(mean(x_proj),mean(y_proj)) , by=cornerID] } } - cornerCoord <- data.frame(cornerCoord[ , c("Xmean","Ymean","Xrel","Yrel","cornerID")]) + cornerCoord <- data.frame(cornerCoord[ , c("Xmean","Ymean","x_rel","y_rel","cornerID")]) cornerCoord <- unique(cornerCoord) - colnames(cornerCoord) <- c("X", "Y","Xrel","Yrel","cornerID") + colnames(cornerCoord) <- c("x_proj", "y_proj","x_rel","y_rel","cornerID") } else { # if exactly 1 measures for each corner cornerCoord <- data.frame(cbind(projCoord[,1:2], relCoord[,1:2])) - colnames(cornerCoord) <- c("X","Y","Xrel","Yrel") + colnames(cornerCoord) <- c("x_proj","y_proj","x_rel","y_rel") outliers <- data.frame() if(!is.null(cornerID)) cornerCoord$cornerID = cornerID } @@ -187,7 +187,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc } else { # trustGPScorners = "FALSE" cornerCoord <- data.table(cbind(projCoord[,1:2], relCoord[,1:2])) - setnames(cornerCoord, c("X","Y","Xrel","Yrel")) + setnames(cornerCoord, c("x_proj","y_proj","x_rel","y_rel")) cornerCoord[, nRow := .I] # Transformation of relCoord to projected coordinates by a procrust analyses @@ -199,7 +199,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc cornerCoord[ , outlier := ifelse( sqrt((procrustCoord[, 1] - projCoord[, 1])^2 + (procrustCoord[, 2] - projCoord[, 2])^2) > maxDist, T, F)] - outliers <- cornerCoord[outlier==T , c("X","Y","nRow")] + outliers <- cornerCoord[outlier==T , c("x_proj","y_proj","nRow")] if (length(outliers)==nrow(projCoord)){ stop("All coordinates points are considered as outliers at the first stage.\n @@ -219,7 +219,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc sqrt((procrustCoord[, 1] - projCoord[-outliers$nRow, 1])^2 + (procrustCoord[, 2] - projCoord[-outliers$nRow, 2])^2) > maxDist, T, F)] - outliers <- rbind(outliers , cornerCoord[outlier==T , c("X","Y","nRow")]) + outliers <- rbind(outliers , cornerCoord[outlier==T , c("x_proj","y_proj","nRow")]) cornerCoord <- cornerCoord[outlier==F ,] if(nrow(cornerCoord)+nrow(outliers) == nrow(projCoord)) refineCoord <- F } @@ -233,36 +233,36 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc cornerProjCoord <- sweep(cornerProjCoord, 2, procrustRes$translation, FUN = "+") cornerCoord <- data.frame(cbind(cornerProjCoord[,1:2], cornerRelCoord[,1:2])) - colnames(cornerCoord) <- c("X","Y","Xrel","Yrel") + colnames(cornerCoord) <- c("x_proj","y_proj","x_rel","y_rel") if(!is.null(cornerID)) cornerCoord$cornerID <- unique(cornerID) } # End trustGPScorners = "FALSE" # 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]) + centroid <- colMeans(cornerCoord[,c("x_rel","y_rel")]) + angles <- base::atan2(cornerCoord[["y_rel"]] - centroid[2], cornerCoord[["x_rel"]] - centroid[1]) cornerCoord <- cornerCoord[order(angles), ] # Create a polygon ----------------------------------------------------------- - cornerPolygon <- st_multipoint(as.matrix(rbind(cornerCoord[,c("X","Y")], cornerCoord[1,c("X","Y")]))) + cornerPolygon <- st_multipoint(as.matrix(rbind(cornerCoord[,c("x_proj","y_proj")], cornerCoord[1,c("x_proj","y_proj")]))) cornerPolygon <- st_polygon(x = list(cornerPolygon), dim = "XY") cornerPolygon <- st_sfc(list(cornerPolygon)) # 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))) { + if(any(! treeCoord[,1] %between% range(cornerCoord$x_rel)) | any(! treeCoord[,2] %between% range(cornerCoord$y_rel))) { warning("Be careful, one or more trees are not inside the plot defined by relCoord") } if(trustGPScorners) { 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") + from_corner_coord = cornerCoord[,c("x_rel","y_rel")], + to_corner_coord = cornerCoord[,c("x_proj","y_proj")], + ordered_corner = T) + colnames(treeProjCoord) <- c("x_proj","y_proj") } else { treeProjCoord <- as.matrix(treeCoord[,1:2]) %*% procrustRes$rotation treeProjCoord <- data.frame(sweep(treeProjCoord, 2, procrustRes$translation, FUN = "+")) - colnames(treeProjCoord) <- c("X","Y") + colnames(treeProjCoord) <- c("x_proj","y_proj") } } @@ -270,9 +270,9 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc if (drawPlot) { projCoordPlot <- projCoord[,c(1,2)] - setnames(projCoordPlot , c("X","Y")) + setnames(projCoordPlot , c("x_proj","y_proj")) projCoordPlot$whatpoint <- "GPS measurements" - cornerCoordPlot <- cornerCoord[,c("X","Y")] + cornerCoordPlot <- cornerCoord[,c("x_proj","y_proj")] cornerCoordPlot$whatpoint <- "Reference corners" projCoordPlot <- rbind(projCoordPlot,cornerCoordPlot) if(!is.null(treeCoord)) { @@ -280,20 +280,20 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc projCoordPlot <- rbind(projCoordPlot,treeProjCoord) } projCoordPlot$whatpoint <- factor(projCoordPlot$whatpoint, levels = c("GPS measurements","Outliers (discarded)","Reference corners","Trees")) - arrowPlot <- data.frame(X = rep(cornerCoord$X[1]), Y =rep(cornerCoord$Y[1]), - Xend = c(cornerCoord$X[1]+(cornerCoord$X[4]-cornerCoord$X[1])/4,cornerCoord$X[1]+(cornerCoord$X[2]-cornerCoord$X[1])/4), - Yend = c(cornerCoord$Y[1]+(cornerCoord$Y[4]-cornerCoord$Y[1])/4,cornerCoord$Y[1]+(cornerCoord$Y[2]-cornerCoord$Y[1])/4)) + arrowPlot <- data.frame(X = rep(cornerCoord$x_proj[1]), Y =rep(cornerCoord$y_proj[1]), + Xend = c(cornerCoord$x_proj[1]+(cornerCoord$x_proj[4]-cornerCoord$x_proj[1])/4,cornerCoord$x_proj[1]+(cornerCoord$x_proj[2]-cornerCoord$x_proj[1])/4), + Yend = c(cornerCoord$y_proj[1]+(cornerCoord$y_proj[4]-cornerCoord$y_proj[1])/4,cornerCoord$y_proj[1]+(cornerCoord$y_proj[2]-cornerCoord$y_proj[1])/4)) if(exists("outliers") && nrow(outliers)!=0) { projCoordPlot$whatpoint[outliers$nRow] <- "Outliers (discarded)" } plotDesign <- ggplot2::ggplot(data = projCoordPlot) + - geom_point(aes(x=X,y=Y,col=whatpoint,shape = whatpoint),size=2) + + geom_point(aes(x=x_proj,y=y_proj,col=whatpoint,shape = whatpoint),size=2) + scale_shape_manual(values=c(2,4,15,1), drop=FALSE) + scale_color_manual(values=c('black','red',"black","darkgreen"), drop = FALSE) + - geom_polygon(data = cornerPolygon[[1]][[1]][,] , mapping = aes(x=X,y=Y), colour="black",fill=NA,linewidth=1.2)+ + geom_polygon(data = cornerPolygon[[1]][[1]][,] , mapping = aes(x=x_proj,y=y_proj), colour="black",fill=NA,linewidth=1.2)+ geom_segment(data=arrowPlot, mapping=aes(x=X,y=Y,xend=Xend,yend=Yend), arrow=arrow(length=unit(0.4,"cm")),linewidth=1, col="blue") + - geom_text(data = arrowPlot, mapping = aes(x=Xend,y=Yend,label=c("Xrel","Yrel")), hjust=c(0,-0.3), vjust=c(-0.5,0), col="blue" ) + + geom_text(data = arrowPlot, mapping = aes(x=Xend,y=Yend,label=c("x_rel","y_rel")), hjust=c(0,-0.3), vjust=c(-0.5,0), col="blue" ) + ggtitle("Plot display") + theme_minimal() + theme(legend.title=element_blank())+ @@ -321,7 +321,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc output$codeUTM <- codeUTM } if (!is.null(treeCoord)) { - output$treeProjCoord <- treeProjCoord[,c("X","Y")] + output$treeProjCoord <- treeProjCoord[,c("x_proj","y_proj")] } return(output) } diff --git a/R/divide_plot.R b/R/divide_plot.R index bba2d8f..6634abf 100644 --- a/R/divide_plot.R +++ b/R/divide_plot.R @@ -6,8 +6,8 @@ #' @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 rel_coord a matrix or 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 matrix or 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 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 @@ -16,12 +16,12 @@ #' @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_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 +#' - `x_rel` and `y_rel` : the relative X-axis and Y-axis coordinates of subplots corners. +#' - `x_proj` and `y_proj` : 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 +#' - `x_rel` and `y_rel` : the relative X-axis and Y-axis coordinates +#' - `subplot_id` : the subplot associated #' - the rest of tree_coord columns supplied #' #' @export @@ -36,10 +36,10 @@ #' # 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) +#' subplots <- divide_plot(rel_coord, proj_coord, grid_size = 100) #' #' 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) +#' subplots <- divide_plot(rel_coord, proj_coord, grid_size = 100, tree_coord = tree_coord) #' #' # Dealing with multiple plots #' rel_coord <- rbind(rel_coord, rel_coord) @@ -51,12 +51,12 @@ #' divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner = NULL, tree_coord = NULL, plot_ID_tree = NULL) { - # Parameters verification ---------------------------------------------------- + # Checking parameters -------------------------------------------------------- if (is.matrix(rel_coord)) { - rel_coord <- data.frame(x_rel=rel_coord[,1],y_rel=rel_coord[,2]) + rel_coord <- data.frame(rel_coord) } if (!is.null(proj_coord) && !is.data.frame(proj_coord)) { - proj_coord <- data.frame(x_proj=proj_coord[,1], y_proj=proj_coord[,2]) + proj_coord <- data.frame(proj_coord) } 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") @@ -90,18 +90,17 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner } - # Formatting data ----------------------------------------------------------- + # Data processing ------------------------------------------------------------ if(length(grid_size)!=2) grid_size = rep(grid_size,2) - rel_coord <- as.data.table(rel_coord) + rel_coord <- data.table(x_rel = rel_coord[[1]], y_rel = rel_coord[[2]]) 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) { @@ -116,7 +115,9 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner } if(!is.null(proj_coord)) { - combine_coord <- cbind(rel_coord[,1:3] , proj_coord) + combine_coord <- cbind(rel_coord, + data.table(x_proj = proj_coord[[1]], y_proj = proj_coord[[2]]) + ) combine_coord <- combine_coord[, sort_rows(.SD), by = plot_ID_corner] } else { combine_coord <- rel_coord[, sort_rows(.SD), by = plot_ID_corner] @@ -128,32 +129,34 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner divide_plot_fct <- function(dat, grid_size) { # Controlling for a regular grid - if((max(dat[[2]])-min(dat[[2]])) %% grid_size[1] != 0) { + if((max(dat[["x_rel"]])-min(dat[["x_rel"]])) %% 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) { + if((max(dat[["y_rel"]])-min(dat[["y_rel"]])) %% 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]) + x_rel = seq(min(dat[["x_rel"]]), max(dat[["x_rel"]]), by = grid_size[1]), + y_rel = seq(min(dat[["y_rel"]]), max(dat[["y_rel"]]), 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 <- rbindlist(apply(plot_grid[x_rel < max(x_rel) & y_rel < max(y_rel),], 1, function(grid_dat) { + X <- as.numeric(grid_dat[["x_rel"]]) + Y <- as.numeric(grid_dat[["y_rel"]]) 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)] + (x_rel == X & y_rel == Y) | (x_rel == X + grid_size[1] & y_rel == Y) | (x_rel == X + grid_size[1] & y_rel == Y + grid_size[2]) | (x_rel == X & y_rel == Y + grid_size[2]), + .(subplot_id = paste(plot_ID_corner, (X-min(plot_grid$x_rel)) / grid_size[1], (Y-min(plot_grid$y_rel)) / grid_size[2], sep = "_"),x_rel, y_rel)] })) - setnames(plot_grid,2:3,names(combine_coord)[2:3]) + + # Sorting rows + plot_grid <- plot_grid[, sort_rows(.SD), by=subplot_id] # 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)) + plot_grid <- cbind(plot_grid,bilinear_interpolation(coord = plot_grid[,c("x_rel","y_rel")] , from_corner_coord = dat[,c("x_rel","y_rel")] , to_corner_coord = dat[,c("x_proj","y_proj")], ordered_corner = T)) } return(plot_grid) } @@ -166,16 +169,16 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner tree_coord <- as.data.table(tree_coord) if(!is.null(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 + setcolorder(tree_coord , "plot_id" , after = 2) } else{ tree_coord[,plot_id:=""] - setcolorder(tree_coord , "plot_id" , after = 2) # if tree_coord has more than 2 col + setcolorder(tree_coord , "plot_id" , after = 2) } invisible(lapply(split(sub_corner_coord, by = "subplot_id", keep.by = TRUE), function(dat) { 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]]), + get(names(tree_coord)[1]) %between% range(dat[["x_rel"]]) & + get(names(tree_coord)[2]) %between% range(dat[["y_rel"]]), subplot_id := dat$subplot_id[1] ] })) @@ -191,7 +194,7 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner } - # returns -------------------------------------------------------------------- + # Returns -------------------------------------------------------------------- if(is.null(plot_ID_corner)) { sub_corner_coord[ , c("subplot_id","plot_ID_corner") := list(paste0("subplot",subplot_id),NULL)] diff --git a/tests/testthat/test-checkCoordPlot.R b/tests/testthat/test-checkCoordPlot.R index 8579e9b..add2168 100644 --- a/tests/testthat/test-checkCoordPlot.R +++ b/tests/testthat/test-checkCoordPlot.R @@ -131,14 +131,14 @@ 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_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")]) + expect_equivalent(res_SW$cornerCoord[c("x_rel","y_rel")] , res_NE$cornerCoord[,c("x_rel","y_rel")]) + expect_equivalent(res_SW$cornerCoord[c("x_proj","y_proj","cornerID")] , res_NE$cornerCoord[c(3,4,1,2),c("x_proj","y_proj","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")] , res_SW_200$cornerCoord[c("X","Y","cornerID")]) - expect_equal(res_SW$cornerCoord[c("Xrel","Yrel")]+200 , res_SW_200$cornerCoord[c("Xrel","Yrel")]) + expect_equal(res_SW$cornerCoord[c("x_proj","y_proj","cornerID")] , res_SW_200$cornerCoord[c("x_proj","y_proj","cornerID")]) + expect_equal(res_SW$cornerCoord[c("x_rel","y_rel")]+200 , res_SW_200$cornerCoord[c("x_rel","y_rel")]) }) @@ -153,8 +153,8 @@ test_that("checkPlotCoord, tree coordinates", { ) cornerID <- rep(c("SW","NW","SE","NE"),e=7) treeCoord <- data.frame( - X = c(10,20,30), - Y = c(10,20,30) + x_proj = c(10,20,30), + y_proj = c(10,20,30) ) res <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F, treeCoord = treeCoord) expect_is(res$treeProjCoord, "data.frame") diff --git a/tests/testthat/test-dividePlot.R b/tests/testthat/test-dividePlot.R index d2b0b09..563ba4f 100644 --- a/tests/testthat/test-dividePlot.R +++ b/tests/testthat/test-dividePlot.R @@ -33,38 +33,38 @@ 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_rel=c(0,50,0,50), y_rel=c(0,0,50,50))) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x_rel=c(0,50,50,0), 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))) 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))) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x_rel=c(0,50,50,0), y_rel=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))) + expect_equal(subplots[1:4,] , data.frame(subplot_id=rep("subplot_0_0",4),x_rel=c(0,50,50,0), y_rel=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))) + expect_equivalent(rect_subplots[29:32,] , data.frame(subplot_id=rep("subplot_3_1",4),x_rel=c(75,100,100,75), y_rel=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))) + expect_equivalent(subplots[13:16,] , data.frame(subplot_id=rep("subplot_1_1",4),x=c(60,110,110,60), 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))) + expect_equivalent(multiple_subplots[29:32,] , data.frame(plotID=rep("plot2",4),subplot_id=rep("plot2_1_1",4),x_rel=c(60,110,110,60), 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))) + #ggplot(subplots,aes(x=x_rel,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,100,50), y=c(0,0,25,25))) }) test_that("divide_plot with projected coordinates", { @@ -76,47 +76,46 @@ test_that("divide_plot with projected coordinates", { 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() + #ggplot(subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x_rel,y=y_rel)) + 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)), + x_proj=c(128.3013,171.6025,146.6025,103.3013), + y_proj=c(178.3013,203.3013,246.6025,221.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")] , + expect_equivalent(subplots[13:16,c("subplot_id","x_proj","y_proj")] , 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)), + x_proj=c(128.3013,171.6025,146.6025,103.3013), + y_proj=c(178.3013,203.3013,246.6025,221.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")) + expect_equal(names(subplots) , c("subplot_id", "x_rel", "y_rel", "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() + #ggplot(rect_subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x_rel,y=y_rel)) + 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)), + data.frame(subplot_id=rep("subplot_1_2",4), x_proj=c(128.3013,171.6025,159.1025,115.8013), y_proj=c(178.3013,203.3013,224.9519,199.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() + #ggplot(rect_subplots,aes(x=x_proj,y=y_proj,label=subplot_id)) + geom_point() + geom_point(aes(x=x_rel,y=y_rel)) + 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)]), + x_proj=c(128.3013,171.6025,159.1025,115.8013)[c(3,4,1,2)], + y_proj=c(178.3013,203.3013,224.9519,199.9519)[c(3,4,1,2)]), tol=1e-5) # Test multiple plots @@ -132,10 +131,11 @@ 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_rel,y=y_rel)) + 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)), + x_rel=c(100,50), y_rel=c(100,100), + x_proj=c(286.6025,103.3013), y_proj=c(246.6025,135.0000)), tol=1e-5) }) @@ -147,10 +147,16 @@ test_that("divide_plot with tree coordinates", { # 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))) + # Test when there is no tree in a subplot + expect_equal(divide_plot(rel_coord = rel_coord_1, grid_size = 50, tree_coord = data.frame(x=1,y=1))$tree_coord , data.frame(x=1,y=1,subplot_id="subplot_0_0")) + + + # basic test 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 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) From 0684345f8760b68f75e8ba2ae412bfeb029bd9e0 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 14 Nov 2024 14:30:43 +0100 Subject: [PATCH 02/16] Commit in order to change computer, explanation will follow in next commit --- R/divide_plot.R | 7 +- R/subplot_summary.R | 134 ++++++++++++++++++ tests/testthat/test-subplot_summary.R | 7 + ...alized_trees_and_forest_stands_metrics.Rmd | 112 ++++++++++++--- 4 files changed, 240 insertions(+), 20 deletions(-) create mode 100644 R/subplot_summary.R create mode 100644 tests/testthat/test-subplot_summary.R diff --git a/R/divide_plot.R b/R/divide_plot.R index 6634abf..f45b0b1 100644 --- a/R/divide_plot.R +++ b/R/divide_plot.R @@ -164,17 +164,15 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner # 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 ------------------------------------------------ + # 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:=plot_ID_tree] - setcolorder(tree_coord , "plot_id" , after = 2) } else{ tree_coord[,plot_id:=""] - setcolorder(tree_coord , "plot_id" , after = 2) } - invisible(lapply(split(sub_corner_coord, by = "subplot_id", keep.by = TRUE), function(dat) { tree_coord[ plot_id == dat$plot_ID_corner[1] & get(names(tree_coord)[1]) %between% range(dat[["x_rel"]]) & @@ -189,6 +187,7 @@ divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, plot_ID_corner tree_coord[ , c("subplot_id","plot_id") := list(paste0("subplot",subplot_id),NULL)] setcolorder(tree_coord , "subplot_id" , after = 2) } else { + setcolorder(tree_coord , "plot_id" , after = 2) setcolorder(tree_coord , "subplot_id" , after = 3) } } diff --git a/R/subplot_summary.R b/R/subplot_summary.R new file mode 100644 index 0000000..1079c7b --- /dev/null +++ b/R/subplot_summary.R @@ -0,0 +1,134 @@ +#' Summarize and display tree information by subplot +#' +#' @description +#' +#' @details +#' +#' +#' @param subplots output of the function [divide_plot()] +#' @param value a character indicating the column in subplots$tree_ooord to display +#' @param fun the function to be applied +#' @param draw_plot a logical indicating wether the plot design should be displayed and returned +#' @param display_trees a logical indicating wether trees should be displayed +#' +#' @return a data frame where: +#' - `plot`: the code of the plot +#' +#' If the `subplot` is set, the output is a list with the previous data frame and a simple features (sf) geometry object. +#' +#' +#' @export +#' +#' @importFrom data.table data.table := first setDT +#' @importFrom grDevices terrain.colors +#' @importFrom graphics segments +#' +#' @examples +#' +#' +subplot_summary <- function(subplots, value=NULL, function_name="sum", draw_plot = TRUE, display_trees = FALSE) { + + # Checking parameters -------------------------------------------------------- + fun <- eval(parse(text = function_name)) + if(is.data.frame(subplots)) { + stop("subplots argument does\'nt contain any tree data frame. Use the divide_plot() function with a non-null tree_coord argument") + } + if (!is.list(subplots) && any(names(subplots)!=c("sub_corner_coord","tree_coord"))) { + stop("subplots argument must be the output of the divide_plot_function(), with a non-null tree_coord argument") + } + if(any(!c("x_proj","y_proj") %in% names(subplots$sub_corner_coord))) { + subplots$sub_corner_coord[,c("x_proj","y_proj")] <- subplots$sub_corner_coord[,c("x_rel","y_rel")] + #stop("projected coordinates are not found in subplots. Use the divide_plot() function with a non-null proj_coord argument (it can be proj_coord = rel_coord if you don\'t have any projected coordinates).") + } + if(!is.null(value) && !value %in% names(subplots$tree_coord)) { + stop("value is not a column name of subplots$tree_coord") + } + if(!is.function(fun)) { + stop(paste("your function",function_name, "is not a function")) + } + if(display_trees == T & draw_plot==F) { + stop("display_trees cannot be TRUE if draw_plot is FALSE") + } + if(display_trees == T & any(!c("x_proj","y_proj") %in% names(subplots$tree_coord))) { + stop("the tree_coord data-frame in subplots doesn\'t the columns x_proj and y_proj which are required to display the trees. Use the output of divide_plot() function with a non-null proj_coord argument (it can be proj_coord = rel_coord if you don\'t have any projected coordinates).") + } + + # Data processing ------------------------------------------------------------ + corner_dat <- data.table(subplots$sub_corner_coord) + tree_summary <- data.table(subplots$tree_coord)[, fun(get(value)) , by=c("subplot_id")] + tree_summary[,plot_id:=sapply(strsplit(tree_summary$subplot_id,"_"),function(x)x[[1]])]#Add plot_id to be able to loop on it + + sf_polygons <- do.call(rbind,lapply(split(corner_dat, by = "subplot_id"), function(data) { + + # creating polygon + mat <- data[, .(x_proj, y_proj)] + mat <- as.matrix(rbind(mat, mat[1, ])) + subplot_polygon <- sf::st_polygon(list(mat)) + + # value summarized by the function + value_fun <- tree_summary[subplot_id == unique(data$subplot_id), V1] + if(identical(value_fun,numeric(0))) value_fun <- 0 + + # creating data-frame to associate with the polygon + df_polygon <- data.frame(unique(data$subplot_id), + value_fun, + value_fun / sf::st_area(subplot_polygon) * 10000) + names(df_polygon) <- c("subplot_id", + paste(value, function_name, sep = "_"), + paste(value, function_name, "per_ha", sep = "_")) + + sf::st_sf(list(subplot_polygon),df_polygon) + })) + + if(draw_plot) { + + plot_list <- lapply(unique(tree_summary$plot_id), function(plot_id) { + + plot_design <- ggplot() + + geom_sf(data = sf_polygons[grep(plot_id,sf_polygons$subplot_id),], + mapping = aes(fill=.data[[names(sf_polygons)[3]]])) + + theme_minimal() + + scale_fill_gradientn(colours = rev(terrain.colors(20))) + + theme(legend.position="bottom", #legend.title = element_blank(), + legend.key.size = unit(0.8,"cm"), + axis.title = element_blank() + ) + if(display_trees) { + plot_design <- plot_design + + geom_point(data = subplots$tree_coord[grep(plot_id,subplots$tree_coord$subplot_id),], + mapping = aes(x = x_proj, y = y_proj), + shape = 3, size=0.5) + } + if(plot_id == "subplot") { + plot_design <- plot_design + ggtitle(paste(function_name,"of",value,"per ha")) + } else { + plot_design <- plot_design + ggtitle(paste(plot_id,":",function_name,"of",value,"per ha")) + + } + + print(plot_design) + plot_design + }) + + } + + setnames(tree_summary, "V1", paste(function_name,value,sep="_")) + tree_summary <- data.frame(tree_summary) + tree_summary[[paste(function_name,value,"per_ha",sep="_")]] <- sf_polygons[[3]] + + if(all(tree_summary$plot_id=="subplot")) { + tree_summary$plot_id <- NULL + plot_list <- plot_list[[1]] + } else{ + tree_summary <- tree_summary[c(3,1,2,4)] + } + + output <- list(tree_summary = tree_summary, polygon = sf_polygons) + + if(draw_plot) { + output$plot_design <- plot_list + } + + return(output) +} + diff --git a/tests/testthat/test-subplot_summary.R b/tests/testthat/test-subplot_summary.R new file mode 100644 index 0000000..9a9f8e2 --- /dev/null +++ b/tests/testthat/test-subplot_summary.R @@ -0,0 +1,7 @@ +context("subplot_summary") + +test_that("subplot_summary error", { + +}) + +# Tester quand il n'y a pas d'arbres dans un sous-plot \ 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 7d45177..7cd4ef7 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.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) @@ -126,7 +127,7 @@ check_plot_trust_field <- checkPlotCoord( We can see that the corners of the plot do not correspond to the GPS measurements. In fact, they correspond to the best compromise between the shape and dimensions of the plot and the GPS measurements. -## Tree coordinates in the projected/geographic coordinate system +## Tree coordinates in the projected/geographic coordinate system {#tree_coordinates} Tree coordinates are almost always measured in the relative (local/field) coordinate system. To retrieve them in the projected system, you can supply their relative coordinates using the `treeCoord` argument. @@ -153,7 +154,12 @@ Note the difference in corner and tree positions between the two situations. 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} +trees <- trees %>% mutate(check_plot_trust_GPS$treeProjCoord) +``` + + +You can also grep and modify the plot via the `$plotDesign` output which is a ggplot object. For example, to change the plot title : ```{r} plot_to_change <- check_plot_trust_GPS$plotDesign @@ -161,7 +167,7 @@ plot_to_change$labels$title <- "A nice title !" plot_to_change ``` -If you want to retrieve the GPS coordinates of the trees in a longitude/latitude format, see this line of code: +If you want to retrieve the GPS coordinates of the trees in a longitude/latitude format, use this code: ```{r} treeGPSCoord <- as.data.frame( proj4::project(check_plot_trust_GPS$treeProjCoord, proj = check_plot_trust_GPS$codeUTM, inverse = TRUE) ) @@ -174,43 +180,117 @@ WORK IN PROGRESS # Dividing plots -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) : +Dividing plots into several rectangular sub-plots is performed using the `divide_plot()` function. This function takes the corner relative coordinates of the plot to divide it into a grid (be aware that **the plot must be rectangular in the relative coordinates system**, 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) + rel_coord = check_plot_trust_GPS$cornerCoord[c("x_rel","y_rel")], + proj_coord = check_plot_trust_GPS$cornerCoord[c("x_proj","y_proj")], + grid_size = 25 # or c(50,50) ) + +kable(head(subplots, 10), digits = 1, row.names = FALSE, caption = "Head of the divide_plot() returns") ``` -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` : +If you don't have any projected coordinates, just let `proj_coord` = NULL. + +For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners but can also **assign to each tree its subplot** with the **`tree_coord` argument** : ```{r divide_plot_trees} 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")], - proj_coord = check_plot_trust_field$cornerCoord[c("X","Y")], #optional - grid_size = 50, # or c(50,50) + rel_coord = check_plot_trust_GPS$cornerCoord[c("x_rel","y_rel")], + proj_coord = check_plot_trust_GPS$cornerCoord[c("x_proj","y_proj")], + grid_size = 25, # 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 now returns a list containing : + +* `sub_corner_coord` : coordinates of subplot corners as previously + +* `tree_coord` : the tree data-frame with the subplot_id added in the 3rd column position + +```{r} +kable(head(subplots$tree_coord), digits = 1, row.names = FALSE, caption = "Head of the divide_plot()$tree_coord returns") +``` 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 <- rbind(check_plot_trust_GPS$cornerCoord[c("x_rel","y_rel","x_proj","y_proj")], check_plot_trust_GPS$cornerCoord[c("x_rel","y_rel","x_proj","y_proj")]) multiple_plots$plot_ID = rep(c("NB1","NB2"),e=4) multiple_trees <- rbind(trees, trees %>% dplyr::mutate(plot="NB2")) -mulitple_subplots <- divide_plot( - rel_coord = multiple_plots[c("Xrel","Yrel")], +multiple_subplots <- divide_plot( + rel_coord = multiple_plots[c("x_rel","y_rel")], + proj_coord = multiple_plots[c("x_proj","y_proj")], grid_size = 50, plot_ID_corner = multiple_plots$plot_ID, tree_coord = multiple_trees, plot_ID_tree = multiple_trees$plot ) ``` + +# Summarizing tree metrics at subplot level + +Once you've applied the `divide_plot()` function with a valid `tree_coord` argument, you can summarize any tree metric present in divide_plot()$tree_coord at the subplot level with the subplot_summary() function. + +```{r subplot_summary} +subplot_metric <- subplot_summary(subplots = subplots, + value = "H") # H is the tree metric + +subplot_metric$tree_summary +subplot_metric$polygon +subplot_metric$plot_design +``` + +By default, the function sums the metric per subplot but you can request any valid function with the `function_name` argument : + +```{r subplot_summary_sd} + +foo <- function(x) {return(mean(x,na.))} + +subplot_metric <- subplot_summary(subplots = subplots, + value = "H", #H is the tree metric + function_name = "foo") +subplot_metric <- subplot_summary(subplots = subplots, + value = "H", #H is the tree metric + function_name = "median") +``` + +The output of the function is a list containing : + +* `$tree_summary` : a summary of the metric per subplot + +* `$polygon` : same as `$tree_summary` with sf polygons in last column + +* `$plot_design` : (if draw_plot=T) a ggplot object that can easily be modified + + +** DEUX OPTIONS POUR LA SUITE ** + +* Soit l'utilisateur peut set display_trees = TRUE, et les arbres seront affichés mais dans ce cas, il faut que les colonnes x_proj et y_proj soit bien présentes dans son tree data-frame + +* Soit on laisse l'utilisateur afficher lui-même les arbres car c'est plutôt facile : + +```{r} +custom_plot <- subplot_metric$plot_design +custom_plot + + geom_point(data=trees, mapping = aes(x = x_proj, y = y_proj, size = D), shape=1) + + ggtitle("The title you want") +``` + + +```{r subplot_summary_display_trees} +multiple_subplot_metric <- subplot_summary( + subplots = multiple_subplots, + value = "H", display_trees = T) + + +``` + + + From 053b8a28ca0c7d3e487c8a0f50ef02144ab9e054 Mon Sep 17 00:00:00 2001 From: abailly Date: Fri, 15 Nov 2024 10:26:48 +0100 Subject: [PATCH 03/16] Removed temporary test --- ...Vignette_spatialized_trees_and_forest_stands_metrics.Rmd | 6 ------ 1 file changed, 6 deletions(-) diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index 7cd4ef7..6cc3115 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -250,12 +250,6 @@ subplot_metric$plot_design By default, the function sums the metric per subplot but you can request any valid function with the `function_name` argument : ```{r subplot_summary_sd} - -foo <- function(x) {return(mean(x,na.))} - -subplot_metric <- subplot_summary(subplots = subplots, - value = "H", #H is the tree metric - function_name = "foo") subplot_metric <- subplot_summary(subplots = subplots, value = "H", #H is the tree metric function_name = "median") From 105123c0947cdbc3cf35eda4862f690c4d31a227 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Mon, 18 Nov 2024 18:04:48 +0100 Subject: [PATCH 04/16] Improved subplot_summary documentation and vignette associated --- R/subplot_summary.R | 101 ++++++++---------- ...alized_trees_and_forest_stands_metrics.Rmd | 80 +++++++++----- 2 files changed, 97 insertions(+), 84 deletions(-) diff --git a/R/subplot_summary.R b/R/subplot_summary.R index 1079c7b..d35c6f9 100644 --- a/R/subplot_summary.R +++ b/R/subplot_summary.R @@ -1,109 +1,97 @@ -#' Summarize and display tree information by subplot +#' Summarise and display tree information by subplot #' #' @description +#' After applying the [divide_plot()] function, this function summarises the desired tree metric by sub-plot and displays the plot representation. #' -#' @details -#' -#' -#' @param subplots output of the function [divide_plot()] -#' @param value a character indicating the column in subplots$tree_ooord to display -#' @param fun the function to be applied +#' @param subplots output of the [divide_plot()] function +#' @param value a character indicating the column in subplots$tree_df to display #' @param draw_plot a logical indicating wether the plot design should be displayed and returned -#' @param display_trees a logical indicating wether trees should be displayed -#' -#' @return a data frame where: -#' - `plot`: the code of the plot -#' -#' If the `subplot` is set, the output is a list with the previous data frame and a simple features (sf) geometry object. +#' @param fun the function to be applied +#' @param ... optional arguments to fun #' +#' @return a list containg the following elements : +#' - `tree_summary` : a summary of the metric per subplot +#' - `polygon` : an sf object : simple feature collection of the subplot's polygon +#' - `plot_design` : (if draw_plot=T) a ggplot object (or a list of ggplot objects) that can easily be modified #' #' @export #' -#' @importFrom data.table data.table := first setDT +#' @importFrom data.table data.table := #' @importFrom grDevices terrain.colors -#' @importFrom graphics segments +#' @importFrom ggplot2 ggplot aes geom_sf theme_minimal scale_fill_gradientn theme ggtitle #' #' @examples #' #' -subplot_summary <- function(subplots, value=NULL, function_name="sum", draw_plot = TRUE, display_trees = FALSE) { +subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, ...) { # Checking parameters -------------------------------------------------------- - fun <- eval(parse(text = function_name)) if(is.data.frame(subplots)) { - stop("subplots argument does\'nt contain any tree data frame. Use the divide_plot() function with a non-null tree_coord argument") + stop("subplots argument does\'nt contain any tree data frame. Use the divide_plot() function with a non-null tree_df argument") } - if (!is.list(subplots) && any(names(subplots)!=c("sub_corner_coord","tree_coord"))) { - stop("subplots argument must be the output of the divide_plot_function(), with a non-null tree_coord argument") + if (!is.list(subplots) && any(names(subplots)!=c("sub_corner_coord","tree_df"))) { + stop("subplots argument must be the output of the divide_plot_function(), with a non-null tree_df argument") } if(any(!c("x_proj","y_proj") %in% names(subplots$sub_corner_coord))) { subplots$sub_corner_coord[,c("x_proj","y_proj")] <- subplots$sub_corner_coord[,c("x_rel","y_rel")] - #stop("projected coordinates are not found in subplots. Use the divide_plot() function with a non-null proj_coord argument (it can be proj_coord = rel_coord if you don\'t have any projected coordinates).") + message("projected coordinates are not found in subplots items (no x_proj and y_proj colnames found), tree metric will be summarised in the relative coordinate system") } - if(!is.null(value) && !value %in% names(subplots$tree_coord)) { - stop("value is not a column name of subplots$tree_coord") + if(!is.null(value) && !value %in% names(subplots$tree_df)) { + stop("value is not a column name of subplots$tree_df") } if(!is.function(fun)) { stop(paste("your function",function_name, "is not a function")) - } - if(display_trees == T & draw_plot==F) { - stop("display_trees cannot be TRUE if draw_plot is FALSE") - } - if(display_trees == T & any(!c("x_proj","y_proj") %in% names(subplots$tree_coord))) { - stop("the tree_coord data-frame in subplots doesn\'t the columns x_proj and y_proj which are required to display the trees. Use the output of divide_plot() function with a non-null proj_coord argument (it can be proj_coord = rel_coord if you don\'t have any projected coordinates).") + } else { + fun <- match.fun(fun) } # Data processing ------------------------------------------------------------ corner_dat <- data.table(subplots$sub_corner_coord) - tree_summary <- data.table(subplots$tree_coord)[, fun(get(value)) , by=c("subplot_id")] + tree_summary <- data.table(subplots$tree_df)[, fun(get(value), ...) , by=c("subplot_id")] tree_summary[,plot_id:=sapply(strsplit(tree_summary$subplot_id,"_"),function(x)x[[1]])]#Add plot_id to be able to loop on it - sf_polygons <- do.call(rbind,lapply(split(corner_dat, by = "subplot_id"), function(data) { + sf_polygons <- do.call(rbind,lapply(split(corner_dat, by = "subplot_id"), function(dat) { # creating polygon - mat <- data[, .(x_proj, y_proj)] + mat <- dat[, .(x_proj, y_proj)] mat <- as.matrix(rbind(mat, mat[1, ])) subplot_polygon <- sf::st_polygon(list(mat)) - # value summarized by the function - value_fun <- tree_summary[subplot_id == unique(data$subplot_id), V1] + # value summarised by the function + value_fun <- tree_summary[subplot_id == unique(dat$subplot_id), V1] if(identical(value_fun,numeric(0))) value_fun <- 0 # creating data-frame to associate with the polygon - df_polygon <- data.frame(unique(data$subplot_id), + df_polygon <- data.frame(unique(dat$subplot_id), value_fun, value_fun / sf::st_area(subplot_polygon) * 10000) names(df_polygon) <- c("subplot_id", - paste(value, function_name, sep = "_"), - paste(value, function_name, "per_ha", sep = "_")) + paste(value, "summary", sep = "_"), + paste(value, "summary", "per_ha", sep = "_")) - sf::st_sf(list(subplot_polygon),df_polygon) + sf_polygon <- sf::st_sf(list(subplot_polygon),df_polygon) + sf::st_geometry(sf_polygon) <- "sf_subplot_polygon" + sf_polygon })) + # Plot the plot(s) ----------------------------------------------------------- + if(draw_plot) { plot_list <- lapply(unique(tree_summary$plot_id), function(plot_id) { - plot_design <- ggplot() + - geom_sf(data = sf_polygons[grep(plot_id,sf_polygons$subplot_id),], - mapping = aes(fill=.data[[names(sf_polygons)[3]]])) + + plot_design <- ggplot(sf_polygons[grep(plot_id,sf_polygons$subplot_id),]) + + geom_sf(mapping = aes(fill=.data[[names(sf_polygons)[3]]])) + theme_minimal() + scale_fill_gradientn(colours = rev(terrain.colors(20))) + theme(legend.position="bottom", #legend.title = element_blank(), legend.key.size = unit(0.8,"cm"), axis.title = element_blank() ) - if(display_trees) { - plot_design <- plot_design + - geom_point(data = subplots$tree_coord[grep(plot_id,subplots$tree_coord$subplot_id),], - mapping = aes(x = x_proj, y = y_proj), - shape = 3, size=0.5) - } if(plot_id == "subplot") { - plot_design <- plot_design + ggtitle(paste(function_name,"of",value,"per ha")) + plot_design <- plot_design + ggtitle(paste("Summary of",value,"per ha")) } else { - plot_design <- plot_design + ggtitle(paste(plot_id,":",function_name,"of",value,"per ha")) - + plot_design <- plot_design + ggtitle(paste(plot_id,": summary of",value,"per ha")) } print(plot_design) @@ -112,15 +100,16 @@ subplot_summary <- function(subplots, value=NULL, function_name="sum", draw_plot } - setnames(tree_summary, "V1", paste(function_name,value,sep="_")) + # Adjusting column names and adding summary per hectare + setnames(tree_summary, "V1", paste(value,"summary",sep="_")) tree_summary <- data.frame(tree_summary) - tree_summary[[paste(function_name,value,"per_ha",sep="_")]] <- sf_polygons[[3]] + tree_summary[[paste(value,"summary per_ha",sep="_")]] <- sf_polygons[[3]] - if(all(tree_summary$plot_id=="subplot")) { - tree_summary$plot_id <- NULL - plot_list <- plot_list[[1]] + if(all(tree_summary$plot_id=="subplot")) { # If just one plot : + tree_summary$plot_id <- NULL # delete plot_id column + plot_list <- plot_list[[1]] # unlist the output } else{ - tree_summary <- tree_summary[c(3,1,2,4)] + tree_summary <- tree_summary[,c(3,1,2,4)] } output <- list(tree_summary = tree_summary, polygon = sf_polygons) diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index 1a27a3a..833685b 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -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) @@ -43,7 +42,7 @@ BIOMASS enables users to manage their plots by : * dividing plots into subplots -* summarizing the desired information at subplot level +* summarising the desired information at subplot level # Required data @@ -195,9 +194,11 @@ kable(head(subplots, 10), digits = 1, row.names = FALSE, caption = "Head of the If you don't have any projected coordinates, just let `proj_coord` = NULL. -For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners but can also **assign to each tree its subplot** with the **`tree_df` and `tree_coords` arguments** : +For the purposes of summarising and representing subplots (coming in the next section), the function returns the coordinates of subplot corners but can also **assign to each tree its subplot** with the **`tree_df` and `tree_coords` arguments** : ```{r divide_plot_trees} +trees$AGB <- computeAGB(trees$D, trees$WD, H = trees$H) + subplots <- divide_plot( rel_coord = check_plot_trust_GPS$corner_coord[c("x_rel","y_rel")], proj_coord = check_plot_trust_GPS$corner_coord[c("x_proj","y_proj")], @@ -229,59 +230,82 @@ multiple_subplots <- divide_plot( proj_coord = multiple_plots[c("x_proj","y_proj")], grid_size = 50, corner_plot_ID = multiple_plots$plot_ID, - tree_coord = multiple_trees, + tree_df = multiple_trees, + tree_coords = c("xRel","yRel"), tree_plot_ID = multiple_trees$plot ) ``` -# Summarizing tree metrics at subplot level +# Summarising tree metrics at subplot level -Once you've applied the `divide_plot()` function with a valid `tree_coord` argument, you can summarize any tree metric present in divide_plot()$tree_coord at the subplot level with the subplot_summary() function. +Once you've applied the `divide_plot()` function with a valid `tree_coord` argument, you can summarise any tree metric present in divide_plot()$tree_df at the subplot level with the `subplot_summary()` function. ```{r subplot_summary} -subplot_metric <- subplot_summary(subplots = subplots, - value = "H") # H is the tree metric +# AGB metric was added before applying the divide_plot() function, but it could have been added afterwards with : +# subplots$tree_df$AGB <- computeAGB(trees$D, trees$WD, H = trees$H) -subplot_metric$tree_summary -subplot_metric$polygon -subplot_metric$plot_design +subplot_metric <- subplot_summary(subplots = subplots, + value = "AGB") # AGB is the tree metric ``` -By default, the function sums the metric per subplot but you can request any valid function with the `function_name` argument : +By default, the function sums the metric per subplot but you can request any valid function with the `fun` argument : -```{r subplot_summary_sd} +```{r subplot_summary_quantile} subplot_metric <- subplot_summary(subplots = subplots, - value = "H", #H is the tree metric - function_name = "median") + value = "AGB", + fun = quantile, probs = 0.9) ``` The output of the function is a list containing : * `$tree_summary` : a summary of the metric per subplot -* `$polygon` : same as `$tree_summary` with sf polygons in last column +* `$polygon` : an sf object containing a simple feature collection of the subplot's polygon * `$plot_design` : (if draw_plot=T) a ggplot object that can easily be modified +## Customizing the ggplot -** DEUX OPTIONS POUR LA SUITE ** - -* Soit l'utilisateur peut set display_trees = TRUE, et les arbres seront affichés mais dans ce cas, il faut que les colonnes x_proj et y_proj soit bien présentes dans son tree data-frame - -* Soit on laisse l'utilisateur afficher lui-même les arbres car c'est plutôt facile : +Here are some examples to custom the ggplot of the `subplot_summary()` function : -```{r} +```{r customize_plot} +subplot_metric <- subplot_summary(subplots = subplots, + value = "AGB") custom_plot <- subplot_metric$plot_design +# Change the title : +custom_plot + + labs(title = "Sum of AGB per hectare") +# Change the legend text : +custom_plot + + labs(fill="Sum of AGB per hectare") +# Display trees : custom_plot + - geom_point(data=trees, mapping = aes(x = x_proj, y = y_proj, size = D), shape=1) + - ggtitle("The title you want") + geom_point(data=trees, mapping = aes(x = x_proj, y = y_proj), shape=3) +# Display trees with diameter as size point : +custom_plot + + geom_point(data=trees, mapping = aes(x = x_proj, y = y_proj, size = D), shape=1) + + labs(size = "Tree diameter (m)", fill = "Sum of AGB per hectare") +# Display trees with diameter as size point and height as color point (and a smaller legend on the right) : +custom_plot + + geom_point(data=trees, mapping = aes(x = x_proj, y = y_proj, size = D, col = H), shape=1) + + labs(size = "Tree diameter (m)", col = "Tree size (m)", fill = "Sum of AGB per hectare") + + scale_color_gradientn(colours = c("white","black")) + + theme(legend.position = "right", legend.key.size = unit(0.5, 'cm')) +``` + +## Exporting as a shapefile + +```{r exporting_shapefile} +subplot_polygons <- subplot_metric$polygon +# sf::st_crs(subplot_polygons) +# sf::st_write(obj = , dsn = "../../Documents/Nourragues.shp") ``` +## Dealing with multiple plots + ```{r subplot_summary_display_trees} multiple_subplot_metric <- subplot_summary( subplots = multiple_subplots, - value = "H", display_trees = T) - - + value = "AGB") ``` From f7374c522e8148e689d81239beaa483cb526dcd9 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Tue, 19 Nov 2024 15:32:06 +0100 Subject: [PATCH 05/16] - Fixed NA issue when trees are not assigned to a subplot - Named plot_list output --- R/subplot_summary.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/subplot_summary.R b/R/subplot_summary.R index d35c6f9..006b848 100644 --- a/R/subplot_summary.R +++ b/R/subplot_summary.R @@ -47,7 +47,7 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, # Data processing ------------------------------------------------------------ corner_dat <- data.table(subplots$sub_corner_coord) - tree_summary <- data.table(subplots$tree_df)[, fun(get(value), ...) , by=c("subplot_id")] + tree_summary <- data.table(subplots$tree_df)[!is.na(subplot_id), fun(get(value), ...) , by=c("subplot_id")] tree_summary[,plot_id:=sapply(strsplit(tree_summary$subplot_id,"_"),function(x)x[[1]])]#Add plot_id to be able to loop on it sf_polygons <- do.call(rbind,lapply(split(corner_dat, by = "subplot_id"), function(dat) { @@ -97,6 +97,7 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, print(plot_design) plot_design }) + names(plot_list) <- unique(tree_summary$plot_id) } From dd4e1e1707621e111c0294cd0f012daf665107c9 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Tue, 19 Nov 2024 16:05:00 +0100 Subject: [PATCH 06/16] Fixed issue wrong attribution of tree_metric_summary per hecate Added examples --- R/divide_plot.R | 2 +- R/subplot_summary.R | 17 +++++++++++++---- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/R/divide_plot.R b/R/divide_plot.R index eba548f..eac88c4 100644 --- a/R/divide_plot.R +++ b/R/divide_plot.R @@ -37,7 +37,7 @@ #' #' # Squared plot and projected coordinates associated #' 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)) +#' proj_coord <- data.frame(x_proj = c(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) #' subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100) #' #' # Assigning trees to subplots diff --git a/R/subplot_summary.R b/R/subplot_summary.R index 006b848..730e9bb 100644 --- a/R/subplot_summary.R +++ b/R/subplot_summary.R @@ -1,11 +1,11 @@ #' Summarise and display tree information by subplot #' #' @description -#' After applying the [divide_plot()] function, this function summarises the desired tree metric by sub-plot and displays the plot representation. +#' After applying the [divide_plot()] function, this function summarises with any defined function the desired tree metric by sub-plot and displays the plot representation. #' #' @param subplots output of the [divide_plot()] function #' @param value a character indicating the column in subplots$tree_df to display -#' @param draw_plot a logical indicating wether the plot design should be displayed and returned +#' @param draw_plot a logical indicating whether the plot design should be displayed and returned #' @param fun the function to be applied #' @param ... optional arguments to fun #' @@ -15,12 +15,20 @@ #' - `plot_design` : (if draw_plot=T) a ggplot object (or a list of ggplot objects) that can easily be modified #' #' @export +#' +#' @author Arthur Bailly #' #' @importFrom data.table data.table := #' @importFrom grDevices terrain.colors #' @importFrom ggplot2 ggplot aes geom_sf theme_minimal scale_fill_gradientn theme ggtitle #' #' @examples +#' 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(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) +#' tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200), metric = rnorm(50,10,5)) +#' subplots <- BIOMASS::divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) +#' subplot_summary(subplots , value = "metric", draw_plot = FALSE) # sum summary +#' subplot_summary(subplots , value = "metric", draw_plot = FALSE, fun = median, probs=0.9) # sum summary #' #' subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, ...) { @@ -103,8 +111,9 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, # Adjusting column names and adding summary per hectare setnames(tree_summary, "V1", paste(value,"summary",sep="_")) - tree_summary <- data.frame(tree_summary) - tree_summary[[paste(value,"summary per_ha",sep="_")]] <- sf_polygons[[3]] + tree_summary <- data.frame(tree_summary[order(subplot_id)]) + + tree_summary[[paste(value,"summary per_ha",sep="_")]] <- sf_polygons[[3]][match(x = tree_summary$subplot_id , table = sf_polygons[[1]])] if(all(tree_summary$plot_id=="subplot")) { # If just one plot : tree_summary$plot_id <- NULL # delete plot_id column From 102d00cd645ee8c459a6f2491f6a54d20de616d4 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Tue, 19 Nov 2024 17:58:18 +0100 Subject: [PATCH 07/16] - Added subplot_summary tests - Improved stop messages - Amended documentation --- NAMESPACE | 3 ++ R/subplot_summary.R | 26 ++++++++----- man/divide_plot.Rd | 6 +-- man/subplot_summary.Rd | 43 +++++++++++++++++++++ tests/testthat/test-subplot_summary.R | 55 ++++++++++++++++++++++++++- 5 files changed, 119 insertions(+), 14 deletions(-) create mode 100644 man/subplot_summary.Rd diff --git a/NAMESPACE b/NAMESPACE index f153429..1108b37 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -24,6 +24,7 @@ export(modelHD) export(numberCorner) export(predictHeight) export(retrieveH) +export(subplot_summary) export(summaryByPlot) importFrom(data.table,"%between%") importFrom(data.table,"%chin%") @@ -51,12 +52,14 @@ importFrom(ggplot2,element_text) importFrom(ggplot2,geom_point) importFrom(ggplot2,geom_polygon) importFrom(ggplot2,geom_segment) +importFrom(ggplot2,geom_sf) 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_fill_gradientn) importFrom(ggplot2,scale_shape_manual) importFrom(ggplot2,scale_x_continuous) importFrom(ggplot2,scale_y_continuous) diff --git a/R/subplot_summary.R b/R/subplot_summary.R index 730e9bb..625c559 100644 --- a/R/subplot_summary.R +++ b/R/subplot_summary.R @@ -35,20 +35,20 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, # Checking parameters -------------------------------------------------------- if(is.data.frame(subplots)) { - stop("subplots argument does\'nt contain any tree data frame. Use the divide_plot() function with a non-null tree_df argument") + stop("subplots argument does'nt contain any tree data frame. Use the divide_plot function with a non-null tree_df argument") } - if (!is.list(subplots) && any(names(subplots)!=c("sub_corner_coord","tree_df"))) { - stop("subplots argument must be the output of the divide_plot_function(), with a non-null tree_df argument") + if (is.list(subplots) && (is.null(names(subplots)) || any(names(subplots)!=c("sub_corner_coord","tree_df")))) { + stop("subplots argument must be the output of the divide_plot_function, with a non-null tree_df argument") } if(any(!c("x_proj","y_proj") %in% names(subplots$sub_corner_coord))) { subplots$sub_corner_coord[,c("x_proj","y_proj")] <- subplots$sub_corner_coord[,c("x_rel","y_rel")] - message("projected coordinates are not found in subplots items (no x_proj and y_proj colnames found), tree metric will be summarised in the relative coordinate system") + message("projected coordinates are not found in sub_corner_coord$subplots, tree metric will be summarised in the relative coordinate system") } - if(!is.null(value) && !value %in% names(subplots$tree_df)) { + if(is.null(value) || !value %in% names(subplots$tree_df)) { stop("value is not a column name of subplots$tree_df") } if(!is.function(fun)) { - stop(paste("your function",function_name, "is not a function")) + stop("the function supplied using `fun =` is not a function") } else { fun <- match.fun(fun) } @@ -56,7 +56,10 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, # Data processing ------------------------------------------------------------ corner_dat <- data.table(subplots$sub_corner_coord) tree_summary <- data.table(subplots$tree_df)[!is.na(subplot_id), fun(get(value), ...) , by=c("subplot_id")] - tree_summary[,plot_id:=sapply(strsplit(tree_summary$subplot_id,"_"),function(x)x[[1]])]#Add plot_id to be able to loop on it + if(any(duplicated(tree_summary$subplot_id))) { + stop("the function supplied using `fun` must return a single value") + } + tree_summary[,plot_id:=sapply(strsplit(tree_summary$subplot_id,"_"),function(x)x[[1]])] #Add plot_id to be able to loop on it sf_polygons <- do.call(rbind,lapply(split(corner_dat, by = "subplot_id"), function(dat) { @@ -75,12 +78,13 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, value_fun / sf::st_area(subplot_polygon) * 10000) names(df_polygon) <- c("subplot_id", paste(value, "summary", sep = "_"), - paste(value, "summary", "per_ha", sep = "_")) + paste(value, "summary_per_ha", sep = "_")) sf_polygon <- sf::st_sf(list(subplot_polygon),df_polygon) sf::st_geometry(sf_polygon) <- "sf_subplot_polygon" sf_polygon })) + sf_polygons <- sf_polygons[order(sf_polygons$subplot_id),] # Plot the plot(s) ----------------------------------------------------------- @@ -113,11 +117,13 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, setnames(tree_summary, "V1", paste(value,"summary",sep="_")) tree_summary <- data.frame(tree_summary[order(subplot_id)]) - tree_summary[[paste(value,"summary per_ha",sep="_")]] <- sf_polygons[[3]][match(x = tree_summary$subplot_id , table = sf_polygons[[1]])] + tree_summary[[paste(value,"summary_per_ha",sep="_")]] <- sf_polygons[[3]][match(x = tree_summary$subplot_id , table = sf_polygons[[1]])] if(all(tree_summary$plot_id=="subplot")) { # If just one plot : tree_summary$plot_id <- NULL # delete plot_id column - plot_list <- plot_list[[1]] # unlist the output + if(draw_plot) { + plot_list <- plot_list[[1]] # unlist the output + } } else{ tree_summary <- tree_summary[,c(3,1,2,4)] } diff --git a/man/divide_plot.Rd b/man/divide_plot.Rd index 3b175c0..88773fc 100644 --- a/man/divide_plot.Rd +++ b/man/divide_plot.Rd @@ -17,9 +17,9 @@ divide_plot( ) } \arguments{ -\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{rel_coord}{a matrix or 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}{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{proj_coord}{a matrix or 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, grid cells will be considered as squares.} @@ -64,7 +64,7 @@ subplots <- divide_plot(rel_coord = rel_coord, grid_size = c(100,50)) # Squared plot and projected coordinates associated 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)) +proj_coord <- data.frame(x_proj = c(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100) # Assigning trees to subplots diff --git a/man/subplot_summary.Rd b/man/subplot_summary.Rd new file mode 100644 index 0000000..d6caca9 --- /dev/null +++ b/man/subplot_summary.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/subplot_summary.R +\name{subplot_summary} +\alias{subplot_summary} +\title{Summarise and display tree information by subplot} +\usage{ +subplot_summary(subplots, value = NULL, draw_plot = TRUE, fun = sum, ...) +} +\arguments{ +\item{subplots}{output of the \code{\link[=divide_plot]{divide_plot()}} function} + +\item{value}{a character indicating the column in subplots$tree_df to display} + +\item{draw_plot}{a logical indicating whether the plot design should be displayed and returned} + +\item{fun}{the function to be applied} + +\item{...}{optional arguments to fun} +} +\value{ +a list containg the following elements : +\itemize{ +\item \code{tree_summary} : a summary of the metric per subplot +\item \code{polygon} : an sf object : simple feature collection of the subplot's polygon +\item \code{plot_design} : (if draw_plot=T) a ggplot object (or a list of ggplot objects) that can easily be modified +} +} +\description{ +After applying the \code{\link[=divide_plot]{divide_plot()}} function, this function summarises with any defined function the desired tree metric by sub-plot and displays the plot representation. +} +\examples{ +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(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) +tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200), metric = rnorm(50,10,5)) +subplots <- BIOMASS::divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) +subplot_summary(subplots , value = "metric", draw_plot = FALSE) # sum summary +subplot_summary(subplots , value = "metric", draw_plot = FALSE, fun = median, probs=0.9) # sum summary + + +} +\author{ +Arthur Bailly +} diff --git a/tests/testthat/test-subplot_summary.R b/tests/testthat/test-subplot_summary.R index 9a9f8e2..a14a07e 100644 --- a/tests/testthat/test-subplot_summary.R +++ b/tests/testthat/test-subplot_summary.R @@ -1,7 +1,60 @@ context("subplot_summary") test_that("subplot_summary error", { + rel_coord <- data.frame(x_rel = c(0, 200, 0, 200), y_rel = c(0, 0, 200, 200)) + subplots <- divide_plot(rel_coord, grid_size = 50) + expect_error(subplot_summary(subplots), "subplots argument does'nt contain any tree data frame. Use the divide_plot function with a non-null tree_df argument") + expect_error(subplot_summary(subplots = list(subplots)), "subplots argument must be the output of the divide_plot_function, with a non-null tree_df argument") + expect_error(subplot_summary(subplots = list(sub_corner_coord = subplots)), "subplots argument must be the output of the divide_plot_function, with a non-null tree_df argument") + + tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200), metric = rnorm(50,10,5)) + subplots <- divide_plot(rel_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) + expect_message(subplot_summary(subplots, value = "metric", draw_plot = F), "projected coordinates are not found in sub_corner_coord$subplots, tree metric will be summarised in the relative coordinate system", fixed=TRUE) + + proj_coord <- data.frame(x_proj = c(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) + subplots <- divide_plot(rel_coord, proj_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) + expect_error(subplot_summary(subplots, draw_plot = F) , "value is not a column name of subplots$tree_df", fixed=TRUE) + expect_error(subplot_summary(subplots, value = "toto", draw_plot = F) , "value is not a column name of subplots$tree_df", fixed=TRUE) + + expect_error(subplot_summary(subplots, value = "metric", draw_plot = F, fun = "quantile") , "the function supplied using `fun =` is not a function", fixed=TRUE) + expect_error(subplot_summary(subplots, value = "metric", draw_plot = F, fun = quantile) , "the function supplied using `fun` must return a single value", fixed=TRUE) }) -# Tester quand il n'y a pas d'arbres dans un sous-plot \ No newline at end of file +test_that("subplot_summary", { + set.seed(52) + rel_coord <- data.frame(x_rel = c(0, 200, 0, 200), y_rel = c(0, 0, 200, 200)) + subplots <- divide_plot(rel_coord, grid_size = 50) + tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200), metric = rnorm(50,10,5)) + subplots <- divide_plot(rel_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) + + # Test without proj_coord + res <- suppressMessages(subplot_summary(subplots, value = "metric", draw_plot = F)) + expect_equivalent(res$tree_summary[1:2,] , data.frame(subplot_id = c("subplot_0_0","subplot_0_1"), + metric_summary = c(12.460178,43.79396), + metric_summary_per_ha = c(49.84071, 175.175823)) , tol=1e-5) + expect_equal(res$tree_summary, data.frame(subplot_id=res$polygon[[1]], metric_summary = res$polygon[[2]], metric_summary_per_ha = res$polygon[[3]])) + expect_equivalent(sf::st_polygon(list(matrix(c(0, 0, 50, 0,50, 50, 0, 50, 0, 0), ncol = 2, byrow = T))), res$polygon[[4]][[1]]) + #subplot_summary(subplots, value = "metric", draw_plot = T) + + # Test with proj_coord + proj_coord <- data.frame(x_proj = c(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) + subplots <- divide_plot(rel_coord, proj_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) + res <- subplot_summary(subplots, value = "metric", draw_plot = F) + expect_equivalent(res$tree_summary[1:2,] , data.frame(subplot_id = c("subplot_0_0","subplot_0_1"), + metric_summary = c(12.460178,43.793956), + metric_summary_per_ha = c(49.929337, 175.487313))) + expect_equivalent(sf::st_polygon(list(matrix(c(210,210,253.25,235,228.25,278.25,185,253.25,210,210), ncol = 2, byrow = T))), res$polygon[[4]][[1]]) + #subplot_summary(subplots, value = "metric", draw_plot = T) + + # Test when there isn't a tree in a subplot + tree_df <- tree_df[-6,] + subplots_less_tree <- divide_plot(rel_coord, proj_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) + res_less_tree <- subplot_summary(subplots_less_tree, value = "metric", draw_plot = F) + expect_equivalent(res$tree_summary[-7,] , res_less_tree$tree_summary) + expect_equal(res$polygon$sf_subplot_polygon, res_less_tree$polygon$sf_subplot_polygon) + + # Test with quantile function + res_quantile <- subplot_summary(subplots_less_tree, value = "metric", draw_plot = F, fun = quantile, probs=0.75) + expect_equal(res_quantile$tree_summary[res_quantile$tree_summary$subplot_id=="subplot_3_2","metric_summary"],16.255414,tol=1e-5) +}) From f5a17877da73275751b2362b3ab009a3278d0342 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Wed, 20 Nov 2024 14:37:00 +0100 Subject: [PATCH 08/16] Added version 2.2 news decription in NEWS --- NEWS | 20 +++++++++++++++++++ ...alized_trees_and_forest_stands_metrics.Rmd | 12 ++++------- 2 files changed, 24 insertions(+), 8 deletions(-) diff --git a/NEWS b/NEWS index df7e76f..a56772f 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,23 @@ +* version 2.2 +- Modified bilinear interpolation to generalized bilinear interpolation that can deal with rectangular plot and non (0;0) origin corner plot +- Added check_plot_coord function which will replace correctCoordGPS() : + - calculation of projected coordinates of plot corners depend on trust_GPS_corners argument : if TRUE, a bilinear interpolation is applied (after averaging plot corners and removing outliers) to retrieve tree projected coordinates, if FALSE a procrust analyses is applied as before + - rangeX and rangeY are no longer needed : they are calculated using relative coordinates dimensions of the plot corners + - a ggplot is returned as output instead of a plot of the plot + - calculation of the projected coordinates of the trees is done, using tree_df and tree_coords arguments +- Added divide_plot function which will replace cutPlot() + - No corner numerotation is needed anymore, neither dimX and dimY (the dimensions of the plot) + - plot division is made on relative coordinates and then transformed in projected coordinates (if supplied) with a bilinear interpolation + - trees are attributed to each subplot using tree_df and tree_coords agruments + - Message warning if the grid dimensions don't match the plot dimensions + - Used of the argument grid_tol to control the plot area which is included in the grid + - The grid can be plot-centred using centred_grid argument +- Added subplot_summary which will replace summaryByPlot + - the function can take any tree metric using the value argument + - any valid function can be applied as a summary + - a ggplot is returned as output instead of just plot the plot division +- Added "Spatialized trees and forest stands metrics with BIOMASS" vignette which will replace the "Manage trees and plot coordinates with BIOMASS" vignette + * 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 diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index 833685b..77d4117 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -128,7 +128,7 @@ check_plot_trust_field <- check_plot_coord( We can see that the corners of the plot do not correspond to the GPS measurements. In fact, they correspond to the best compromise between the shape and dimensions of the plot and the GPS measurements. -## Tree coordinates in the projected/geographic coordinate system {#tree_coordinates} +## Retrieving tree projected coordinates {#tree_coordinates} Tree coordinates are almost always measured in the relative (local/field) coordinate system. To retrieve them in the projected system, you can supply their relative coordinates using the `tree_df` and `tree_coords` arguments. @@ -238,7 +238,7 @@ multiple_subplots <- divide_plot( # Summarising tree metrics at subplot level -Once you've applied the `divide_plot()` function with a valid `tree_coord` argument, you can summarise any tree metric present in divide_plot()$tree_df at the subplot level with the `subplot_summary()` function. +Once you've applied the `divide_plot()` function with a `tree_df` argument, you can summarise any tree metric present in divide_plot()$tree_df at the subplot level with the `subplot_summary()` function. ```{r subplot_summary} # AGB metric was added before applying the divide_plot() function, but it could have been added afterwards with : @@ -272,12 +272,9 @@ Here are some examples to custom the ggplot of the `subplot_summary()` function subplot_metric <- subplot_summary(subplots = subplots, value = "AGB") custom_plot <- subplot_metric$plot_design -# Change the title : +# Change the title and legend : custom_plot + - labs(title = "Sum of AGB per hectare") -# Change the legend text : -custom_plot + - labs(fill="Sum of AGB per hectare") + labs(title = "Nourragues plot" , fill="Sum of AGB per hectare") # Display trees : custom_plot + geom_point(data=trees, mapping = aes(x = x_proj, y = y_proj), shape=3) @@ -298,7 +295,6 @@ custom_plot + ```{r exporting_shapefile} subplot_polygons <- subplot_metric$polygon # sf::st_crs(subplot_polygons) -# sf::st_write(obj = , dsn = "../../Documents/Nourragues.shp") ``` From d4fc4dd3a2b7af84f625afd0dfeeae6ffbd03c54 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 21 Nov 2024 16:41:20 +0100 Subject: [PATCH 09/16] suppressed a deprecated stop message in check_plot_coord --- R/check_plot_coord.R | 3 --- 1 file changed, 3 deletions(-) diff --git a/R/check_plot_coord.R b/R/check_plot_coord.R index 6225fea..d818d05 100644 --- a/R/check_plot_coord.R +++ b/R/check_plot_coord.R @@ -126,9 +126,6 @@ check_plot_coord <- function(proj_coord = NULL, longlat = NULL, rel_coord, trust if (trust_GPS_corners==T && (!is.null(proj_coord)) && nrow(proj_coord)!=4 & is.null(corner_ID)) { stop("The argument corner_ID is needed if trust_GPS_corners is TRUE and if multiple measurements of each corner have been realized") } - if (nrow(unique(rel_coord)) != 4) { - stop( paste(nrow(unique(rel_coord)),"unique corners are detected in rel_coord instead of 4")) - } if (!is.null(corner_ID) && length(corner_ID) != nrow(rel_coord)) { stop("corner_ID must be the same length as the number of rows in rel_coord") } From c73ead01b65efa997d0c7bc50407ba3a1b6c6c5e Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 21 Nov 2024 16:41:20 +0100 Subject: [PATCH 10/16] suppressed another deprecated stop message in check_plot_coord --- R/check_plot_coord.R | 7 ------- 1 file changed, 7 deletions(-) diff --git a/R/check_plot_coord.R b/R/check_plot_coord.R index 6225fea..9e46439 100644 --- a/R/check_plot_coord.R +++ b/R/check_plot_coord.R @@ -126,16 +126,9 @@ check_plot_coord <- function(proj_coord = NULL, longlat = NULL, rel_coord, trust if (trust_GPS_corners==T && (!is.null(proj_coord)) && nrow(proj_coord)!=4 & is.null(corner_ID)) { stop("The argument corner_ID is needed if trust_GPS_corners is TRUE and if multiple measurements of each corner have been realized") } - if (nrow(unique(rel_coord)) != 4) { - stop( paste(nrow(unique(rel_coord)),"unique corners are detected in rel_coord instead of 4")) - } if (!is.null(corner_ID) && length(corner_ID) != nrow(rel_coord)) { stop("corner_ID must be the same length as the number of rows in rel_coord") } - if (!is.null(corner_ID) && length(unique(corner_ID)) !=4) { - stop(paste(length((unique(corner_ID))),"unique corners are detected in corner_ID instead of 4")) - } - # function ------------------------------------------------------------------- From 4cb39e26655aacd6e68bb8655a091674ea6c52ec Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Thu, 28 Nov 2024 17:20:06 +0100 Subject: [PATCH 11/16] Minor modification to pass NOTES and WARNINGS Uneval BIOMASS vignette because of non-working tnrs api --- NAMESPACE | 2 + R/bilinear_interpolation.R | 2 +- R/check_plot_coord.R | 1 - R/correctTaxo.R | 2 +- R/divide_plot.R | 8 +- R/globals.R | 19 +++ R/pkgname.R | 136 ------------------ R/subplot_summary.R | 10 +- readme.md => README.md | 0 man/BIOMASS-package.Rd | 136 ------------------ man/check_plot_coord.Rd | 1 - man/divide_plot.Rd | 7 +- man/subplot_summary.Rd | 10 +- vignettes/BIOMASS.Rmd | 4 +- ...alized_trees_and_forest_stands_metrics.Rmd | 2 - 15 files changed, 50 insertions(+), 290 deletions(-) create mode 100644 R/globals.R rename readme.md => README.md (100%) diff --git a/NAMESPACE b/NAMESPACE index 1108b37..27470ad 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -87,6 +87,7 @@ importFrom(sf,st_sfc) importFrom(stats,SSmicmen) importFrom(stats,as.formula) importFrom(stats,coef) +importFrom(stats,dist) importFrom(stats,formula) importFrom(stats,lm) importFrom(stats,median) @@ -105,5 +106,6 @@ importFrom(terra,vect) importFrom(utils,askYesNo) importFrom(utils,data) importFrom(utils,download.file) +importFrom(utils,globalVariables) importFrom(utils,head) importFrom(utils,unzip) diff --git a/R/bilinear_interpolation.R b/R/bilinear_interpolation.R index e451d3d..d6d22a8 100644 --- a/R/bilinear_interpolation.R +++ b/R/bilinear_interpolation.R @@ -66,7 +66,7 @@ bilinear_interpolation = function(coord, from_corner_coord, to_corner_coord, ord } # 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)) { + if(!all(abs(stats::dist(rbind(from_corner_coord[,1:2],centroid))[c(4,7,9,10)] - mean(stats::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") } diff --git a/R/check_plot_coord.R b/R/check_plot_coord.R index 6225fea..4a11083 100644 --- a/R/check_plot_coord.R +++ b/R/check_plot_coord.R @@ -76,7 +76,6 @@ #' draw_plot = TRUE #' ) #' } -#' check_plot_coord <- function(proj_coord = NULL, longlat = NULL, rel_coord, trust_GPS_corners, corner_ID=NULL, max_dist = 15, rm_outliers = TRUE, draw_plot = TRUE, tree_df = NULL, tree_coords=NULL) { # Checking arguments ------------------------------------------------- diff --git a/R/correctTaxo.R b/R/correctTaxo.R index 36069a0..1db1532 100644 --- a/R/correctTaxo.R +++ b/R/correctTaxo.R @@ -45,7 +45,7 @@ if (getRversion() >= "2.15.1") { #' @author Ariane TANGUY, Arthur PERE, Maxime REJOU-MECHAIN, Guillaume CORNU #' #' @examples -#' \donttest{ +#' \dontrun{ #' correctTaxo(genus = "Astrocarium", species = "standleanum") #' correctTaxo(genus = "Astrocarium standleanum") #' } diff --git a/R/divide_plot.R b/R/divide_plot.R index eac88c4..2ad4791 100644 --- a/R/divide_plot.R +++ b/R/divide_plot.R @@ -29,6 +29,7 @@ #' @export #' @author Arthur PERE, Arthur BAILLY #' @importFrom data.table data.table := setcolorder +#' @importFrom stats dist #' @examples #' #' # Rectangular plot and grid cells @@ -42,7 +43,8 @@ #' #' # Assigning trees to subplots #' tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200)) -#' subplots <- divide_plot(rel_coord, proj_coord, 100, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) +#' subplots <- divide_plot(rel_coord, proj_coord, 100, +#' tree_df = tree_df, tree_coords = c("x_tree","y_tree")) #' subplots$sub_corner_coord #' subplots$tree_df #' @@ -56,7 +58,9 @@ #' tree_df <- rbind(tree_df, data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200))) #' corner_plot_ID <- rep(c("plot1","plot2"), e=4) #' tree_plot_ID <- rep(c("plot1","plot2"), e=50) -#' subplots <- divide_plot(rel_coord, proj_coord, 100, tree_df, c("x_tree","y_tree"), corner_plot_ID, tree_plot_ID) +#' subplots <- divide_plot(rel_coord, proj_coord, 100, +#' tree_df, c("x_tree","y_tree"), +#' corner_plot_ID, tree_plot_ID) #' divide_plot <- function(rel_coord, proj_coord = NULL, grid_size, tree_df = NULL, tree_coords = NULL, corner_plot_ID = NULL, tree_plot_ID = NULL, grid_tol = 0.1, centred_grid = F) { diff --git a/R/globals.R b/R/globals.R new file mode 100644 index 0000000..a009081 --- /dev/null +++ b/R/globals.R @@ -0,0 +1,19 @@ +#' @importFrom utils globalVariables + + +if(getRversion() >= "2.15.1") { + # modelHD function + utils::globalVariables(c("x","y")) + + # check_plot_coord function + utils::globalVariables(c("x_proj","y_proj","outlier","Xmean","Ymean","row_number","whatpoint","Xend","Yend")) + + # checkPlotCoord function + utils::globalVariables("nRow") + + # divide_plot function + utils::globalVariables(c("x_rel","y_rel","subplot_id","plot_id")) + + # subplot_summary function + utils::globalVariables(c("subplot_id","plot_id",".data")) +} \ No newline at end of file diff --git a/R/pkgname.R b/R/pkgname.R index e5d0ac8..700233b 100644 --- a/R/pkgname.R +++ b/R/pkgname.R @@ -3,142 +3,6 @@ #' An R Package for estimating above-ground biomass and its uncertainty in tropical forests. #' Methods in Ecology and Evolution, 8(9), 1163-1167. #' -#' -#' @examples -#' \donttest{ -#' library(BIOMASS) -#' -#' # Dataset containing plot inventory data from Karnataka, India (Ramesh et al. 2010) -#' data(KarnatakaForest) -#' str(KarnatakaForest) -#' -#' # Dataset containing height and diameter measurements from two 1-ha plots -#' # established in the lowland rainforest of French Guiana, at the Nouragues -#' # Ecological Research Station -#' data(NouraguesHD) -#' str(NouraguesHD) -#' -#' ############################################################################# -#' # WOOD DENSITY -#' -#' # 1-RETRIEVE AND CORRECT TAXONOMY -#' -#' # Checking typos in taxonomy -#' Taxo <- correctTaxo(genus = KarnatakaForest$genus, species = KarnatakaForest$species) -#' KarnatakaForest$genusCorr <- Taxo$genusCorrected -#' KarnatakaForest$speciesCorr <- Taxo$speciesCorrected -#' -#' # Retrieving APG III Families and Orders from Genus names -#' APG <- getTaxonomy(KarnatakaForest$genusCorr, findOrder = TRUE) -#' KarnatakaForest$familyAPG <- APG$family -#' KarnatakaForest$orderAPG <- APG$order -#' -#' # 2-RETRIEVE WOOD DENSITY -#' dataWD <- getWoodDensity( -#' genus = KarnatakaForest$genusCorr, -#' species = KarnatakaForest$speciesCorr, -#' stand = KarnatakaForest$plotID -#' ) -#' -#' ############################################################################# -#' # TREE HEIGHT -#' -#' # Compare different local H-D models -#' modelHD( -#' D = NouraguesHD$D, H = NouraguesHD$H, -#' drawGraph = TRUE, useWeight = TRUE -#' ) -#' -#' # Compute the local H-D model with the lowest RSE -#' HDmodel <- modelHD( -#' D = NouraguesHD$D, H = NouraguesHD$H, -#' method = "log2", useWeight = TRUE -#' ) -#' -#' # Compute plot-specific H-D models -#' HDmodelPerPlot <- modelHD(NouraguesHD$D, NouraguesHD$H, -#' method = "weibull", -#' useWeight = TRUE, plot = NouraguesHD$plotId -#' ) -#' -#' RSEmodels <- sapply(HDmodelPerPlot, function(x) x$RSE) -#' Coeffmodels <- lapply(HDmodelPerPlot, function(x) x$coefficients) -#' -#' # Retrieve height data from a local HD model -#' dataHlocal <- retrieveH(D = KarnatakaForest$D, model = HDmodel) -#' -#' # Retrieve height data from a Feldpaush et al. (2012) averaged model -#' dataHfeld <- retrieveH(D = KarnatakaForest$D, region = "SEAsia") -#' -#' # Retrieve height data from Chave et al. (2012) equation 6 -#' dataHchave <- retrieveH( -#' D = KarnatakaForest$D, -#' coord = cbind(KarnatakaForest$long, KarnatakaForest$lat) -#' ) -#' -#' ############################################################################# -#' # AGB CALCULATION -#' -#' KarnatakaForest$WD <- dataWD$meanWD -#' KarnatakaForest$H <- dataHlocal$H -#' KarnatakaForest$Hfeld <- dataHfeld$H -#' -#' # Compute AGB(Mg) per tree -#' AGBtree <- computeAGB( -#' D = KarnatakaForest$D, WD = KarnatakaForest$WD, -#' H = KarnatakaForest$H -#' ) -#' -#' # Compute AGB(Mg) per plot -#' AGBplot <- summaryByPlot(AGBtree, KarnatakaForest$plotId) -#' -#' # Compute AGB(Mg) per tree without height information (Eq. 7 from Chave et al. (2014)) -#' AGBplotChave <- summaryByPlot( -#' computeAGB( -#' D = KarnatakaForest$D, WD = KarnatakaForest$WD, -#' coord = KarnatakaForest[, c("long", "lat")] -#' ), -#' plot = KarnatakaForest$plotId -#' ) -#' -#' # Compute AGB(Mg) per tree with Feldpausch et al. (2012) regional H-D model -#' AGBplotFeld <- summaryByPlot( -#' computeAGB( -#' D = KarnatakaForest$D, WD = KarnatakaForest$WD, -#' H = KarnatakaForest$Hfeld -#' ), -#' plot = KarnatakaForest$plotId -#' ) -#' -#' ############################################################################# -#' # PROPAGATING ERRORS -#' -#' KarnatakaForest$sdWD <- dataWD$sdWD -#' KarnatakaForest$HfeldRSE <- dataHfeld$RSE -#' -#' # Per plot using the local HD model constructed above (modelHD) -#' resultMC <- AGBmonteCarlo( -#' D = KarnatakaForest$D, WD = KarnatakaForest$WD, errWD = KarnatakaForest$sdWD, -#' HDmodel = HDmodel, Dpropag = "chave2004" -#' ) -#' resMC <- summaryByPlot(resultMC$AGB_simu, KarnatakaForest$plotId) -#' -#' # Per plot using the Feldpaush regional HD averaged model -#' AGBmonteCarlo( -#' D = KarnatakaForest$D, WD = KarnatakaForest$WD, -#' errWD = KarnatakaForest$sdWD, H = KarnatakaForest$Hfeld, -#' errH = KarnatakaForest$HfeldRSE, Dpropag = "chave2004" -#' ) -#' resMC <- summaryByPlot(resultMC$AGB_simu, KarnatakaForest$plotId) -#' -#' # Per plot using Chave et al. (2014) Equation 7 -#' resultMC <- AGBmonteCarlo( -#' D = KarnatakaForest$D, WD = KarnatakaForest$WD, errWD = KarnatakaForest$sdWD, -#' coord = KarnatakaForest[, c("long", "lat")], -#' Dpropag = "chave2004" -#' ) -#' resMC <- summaryByPlot(resultMC$AGB_simu, KarnatakaForest$plotId) -#' } #' @keywords internal "_PACKAGE" diff --git a/R/subplot_summary.R b/R/subplot_summary.R index 625c559..435b33d 100644 --- a/R/subplot_summary.R +++ b/R/subplot_summary.R @@ -26,9 +26,13 @@ #' 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(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) #' tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200), metric = rnorm(50,10,5)) -#' subplots <- BIOMASS::divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) -#' subplot_summary(subplots , value = "metric", draw_plot = FALSE) # sum summary -#' subplot_summary(subplots , value = "metric", draw_plot = FALSE, fun = median, probs=0.9) # sum summary +#' subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 50, +#' tree_df = tree_df, tree_coords = c("x_tree","y_tree")) +#' # Sum summary (by default) +#' subplot_summary(subplots , value = "metric", draw_plot = FALSE) +#' # 9th quantile summary (for example) +#' subplot_summary(subplots , value = "metric", draw_plot = FALSE, +#' fun = quantile, probs=0.9) #' #' subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, ...) { diff --git a/readme.md b/README.md similarity index 100% rename from readme.md rename to README.md diff --git a/man/BIOMASS-package.Rd b/man/BIOMASS-package.Rd index a55af88..2515063 100644 --- a/man/BIOMASS-package.Rd +++ b/man/BIOMASS-package.Rd @@ -8,142 +8,6 @@ \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, (3) construct height-diameter models, (4) manage tree and plot coordinates, (5) estimate the aboveground biomass/carbon at the stand level with associated uncertainty. To cite 'BIOMASS', please use citation("BIOMASS"). See more in the article of Réjou-Méchain et al. (2017) \doi{10.1111/2041-210X.12753}. } -\examples{ -\donttest{ -library(BIOMASS) - -# Dataset containing plot inventory data from Karnataka, India (Ramesh et al. 2010) -data(KarnatakaForest) -str(KarnatakaForest) - -# Dataset containing height and diameter measurements from two 1-ha plots -# established in the lowland rainforest of French Guiana, at the Nouragues -# Ecological Research Station -data(NouraguesHD) -str(NouraguesHD) - -############################################################################# -# WOOD DENSITY - -# 1-RETRIEVE AND CORRECT TAXONOMY - -# Checking typos in taxonomy -Taxo <- correctTaxo(genus = KarnatakaForest$genus, species = KarnatakaForest$species) -KarnatakaForest$genusCorr <- Taxo$genusCorrected -KarnatakaForest$speciesCorr <- Taxo$speciesCorrected - -# Retrieving APG III Families and Orders from Genus names -APG <- getTaxonomy(KarnatakaForest$genusCorr, findOrder = TRUE) -KarnatakaForest$familyAPG <- APG$family -KarnatakaForest$orderAPG <- APG$order - -# 2-RETRIEVE WOOD DENSITY -dataWD <- getWoodDensity( - genus = KarnatakaForest$genusCorr, - species = KarnatakaForest$speciesCorr, - stand = KarnatakaForest$plotID -) - -############################################################################# -# TREE HEIGHT - -# Compare different local H-D models -modelHD( - D = NouraguesHD$D, H = NouraguesHD$H, - drawGraph = TRUE, useWeight = TRUE -) - -# Compute the local H-D model with the lowest RSE -HDmodel <- modelHD( - D = NouraguesHD$D, H = NouraguesHD$H, - method = "log2", useWeight = TRUE -) - -# Compute plot-specific H-D models -HDmodelPerPlot <- modelHD(NouraguesHD$D, NouraguesHD$H, - method = "weibull", - useWeight = TRUE, plot = NouraguesHD$plotId -) - -RSEmodels <- sapply(HDmodelPerPlot, function(x) x$RSE) -Coeffmodels <- lapply(HDmodelPerPlot, function(x) x$coefficients) - -# Retrieve height data from a local HD model -dataHlocal <- retrieveH(D = KarnatakaForest$D, model = HDmodel) - -# Retrieve height data from a Feldpaush et al. (2012) averaged model -dataHfeld <- retrieveH(D = KarnatakaForest$D, region = "SEAsia") - -# Retrieve height data from Chave et al. (2012) equation 6 -dataHchave <- retrieveH( - D = KarnatakaForest$D, - coord = cbind(KarnatakaForest$long, KarnatakaForest$lat) -) - -############################################################################# -# AGB CALCULATION - -KarnatakaForest$WD <- dataWD$meanWD -KarnatakaForest$H <- dataHlocal$H -KarnatakaForest$Hfeld <- dataHfeld$H - -# Compute AGB(Mg) per tree -AGBtree <- computeAGB( - D = KarnatakaForest$D, WD = KarnatakaForest$WD, - H = KarnatakaForest$H -) - -# Compute AGB(Mg) per plot -AGBplot <- summaryByPlot(AGBtree, KarnatakaForest$plotId) - -# Compute AGB(Mg) per tree without height information (Eq. 7 from Chave et al. (2014)) -AGBplotChave <- summaryByPlot( - computeAGB( - D = KarnatakaForest$D, WD = KarnatakaForest$WD, - coord = KarnatakaForest[, c("long", "lat")] - ), - plot = KarnatakaForest$plotId -) - -# Compute AGB(Mg) per tree with Feldpausch et al. (2012) regional H-D model -AGBplotFeld <- summaryByPlot( - computeAGB( - D = KarnatakaForest$D, WD = KarnatakaForest$WD, - H = KarnatakaForest$Hfeld - ), - plot = KarnatakaForest$plotId -) - -############################################################################# -# PROPAGATING ERRORS - -KarnatakaForest$sdWD <- dataWD$sdWD -KarnatakaForest$HfeldRSE <- dataHfeld$RSE - -# Per plot using the local HD model constructed above (modelHD) -resultMC <- AGBmonteCarlo( - D = KarnatakaForest$D, WD = KarnatakaForest$WD, errWD = KarnatakaForest$sdWD, - HDmodel = HDmodel, Dpropag = "chave2004" -) -resMC <- summaryByPlot(resultMC$AGB_simu, KarnatakaForest$plotId) - -# Per plot using the Feldpaush regional HD averaged model -AGBmonteCarlo( - D = KarnatakaForest$D, WD = KarnatakaForest$WD, - errWD = KarnatakaForest$sdWD, H = KarnatakaForest$Hfeld, - errH = KarnatakaForest$HfeldRSE, Dpropag = "chave2004" -) -resMC <- summaryByPlot(resultMC$AGB_simu, KarnatakaForest$plotId) - -# Per plot using Chave et al. (2014) Equation 7 -resultMC <- AGBmonteCarlo( - D = KarnatakaForest$D, WD = KarnatakaForest$WD, errWD = KarnatakaForest$sdWD, - coord = KarnatakaForest[, c("long", "lat")], - Dpropag = "chave2004" -) -resMC <- summaryByPlot(resultMC$AGB_simu, KarnatakaForest$plotId) -} -} \references{ Réjou-Méchain M., Tanguy A., Piponiot C., Chave J., Hérault B. (2017). BIOMASS : An R Package for estimating above-ground biomass and its uncertainty in tropical forests. diff --git a/man/check_plot_coord.Rd b/man/check_plot_coord.Rd index 99e0f04..4fac915 100644 --- a/man/check_plot_coord.Rd +++ b/man/check_plot_coord.Rd @@ -96,7 +96,6 @@ check_plot_coord( draw_plot = TRUE ) } - } \author{ Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY diff --git a/man/divide_plot.Rd b/man/divide_plot.Rd index 88773fc..ea3028e 100644 --- a/man/divide_plot.Rd +++ b/man/divide_plot.Rd @@ -69,7 +69,8 @@ subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100) # Assigning trees to subplots tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200)) -subplots <- divide_plot(rel_coord, proj_coord, 100, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) +subplots <- divide_plot(rel_coord, proj_coord, 100, + tree_df = tree_df, tree_coords = c("x_tree","y_tree")) subplots$sub_corner_coord subplots$tree_df @@ -83,7 +84,9 @@ proj_coord <- rbind(proj_coord, proj_coord + 200) tree_df <- rbind(tree_df, data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200))) corner_plot_ID <- rep(c("plot1","plot2"), e=4) tree_plot_ID <- rep(c("plot1","plot2"), e=50) -subplots <- divide_plot(rel_coord, proj_coord, 100, tree_df, c("x_tree","y_tree"), corner_plot_ID, tree_plot_ID) +subplots <- divide_plot(rel_coord, proj_coord, 100, + tree_df, c("x_tree","y_tree"), + corner_plot_ID, tree_plot_ID) } \author{ diff --git a/man/subplot_summary.Rd b/man/subplot_summary.Rd index d6caca9..6101625 100644 --- a/man/subplot_summary.Rd +++ b/man/subplot_summary.Rd @@ -32,9 +32,13 @@ After applying the \code{\link[=divide_plot]{divide_plot()}} function, this func 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(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200), metric = rnorm(50,10,5)) -subplots <- BIOMASS::divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 50, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) -subplot_summary(subplots , value = "metric", draw_plot = FALSE) # sum summary -subplot_summary(subplots , value = "metric", draw_plot = FALSE, fun = median, probs=0.9) # sum summary +subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 50, + tree_df = tree_df, tree_coords = c("x_tree","y_tree")) +# Sum summary (by default) +subplot_summary(subplots , value = "metric", draw_plot = FALSE) +# 9th quantile summary (for example) +subplot_summary(subplots , value = "metric", draw_plot = FALSE, + fun = quantile, probs=0.9) } diff --git a/vignettes/BIOMASS.Rmd b/vignettes/BIOMASS.Rmd index 356356a..0986757 100644 --- a/vignettes/BIOMASS.Rmd +++ b/vignettes/BIOMASS.Rmd @@ -12,8 +12,8 @@ vignette: > --- ```{r setup, include=FALSE} -knitr::opts_chunk$set(echo = TRUE, cache = TRUE) -test <- TRUE +knitr::opts_chunk$set(echo = TRUE, cache = TRUE, eval=FALSE) +test <- FALSE CACHE <- TRUE require(knitr) require(BIOMASS) diff --git a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd index 77d4117..07358e4 100644 --- a/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd +++ b/vignettes/Vignette_spatialized_trees_and_forest_stands_metrics.Rmd @@ -24,8 +24,6 @@ require(BIOMASS) require(knitr) require(ggplot2) require(data.table) -require(tidyverse) -require(dplyr) ``` /!\\ VIGNETTE IN PROGRESS /!\\ From f745444108542130969f25afcf1fd8fb9112327f Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Fri, 29 Nov 2024 10:42:50 +0100 Subject: [PATCH 12/16] Replaced pkgname.R by BIOMASS-package.R Moved .onLoad and .onAttach codes in zzz.R --- R/BIOMASS-package.R | 9 +++++++++ R/{pkgname.R => zzz.R} | 10 +--------- 2 files changed, 10 insertions(+), 9 deletions(-) create mode 100644 R/BIOMASS-package.R rename R/{pkgname.R => zzz.R} (66%) diff --git a/R/BIOMASS-package.R b/R/BIOMASS-package.R new file mode 100644 index 0000000..7b37f30 --- /dev/null +++ b/R/BIOMASS-package.R @@ -0,0 +1,9 @@ +#' @keywords internal +"_PACKAGE" + +# The following block is used by usethis to automatically manage +# roxygen namespace tags. Modify with care! +## usethis namespace: start +## usethis namespace: end + +NULL \ No newline at end of file diff --git a/R/pkgname.R b/R/zzz.R similarity index 66% rename from R/pkgname.R rename to R/zzz.R index 700233b..eb64ab9 100644 --- a/R/pkgname.R +++ b/R/zzz.R @@ -1,18 +1,10 @@ -#' @references -#' Réjou-Méchain M., Tanguy A., Piponiot C., Chave J., Hérault B. (2017). BIOMASS : -#' An R Package for estimating above-ground biomass and its uncertainty in tropical forests. -#' Methods in Ecology and Evolution, 8(9), 1163-1167. -#' -#' @keywords internal -"_PACKAGE" - .onLoad <- function(libname, pkgname) { } .onAttach <- function(libname, pkgname) { basePath <- cachePath() - + if(attr(basePath, "source")=="temp") { packageStartupMessage( "Using temporary cache", From c9b665fc98f25dd45f35e47210adeff69a6d66b7 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Fri, 29 Nov 2024 10:43:54 +0100 Subject: [PATCH 13/16] Fixed tests issues --- R/check_plot_coord.R | 6 ++++++ tests/testthat/test-check_coord_plot.R | 8 -------- 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/R/check_plot_coord.R b/R/check_plot_coord.R index 55956c1..60114b2 100644 --- a/R/check_plot_coord.R +++ b/R/check_plot_coord.R @@ -125,9 +125,15 @@ check_plot_coord <- function(proj_coord = NULL, longlat = NULL, rel_coord, trust if (trust_GPS_corners==T && (!is.null(proj_coord)) && nrow(proj_coord)!=4 & is.null(corner_ID)) { stop("The argument corner_ID is needed if trust_GPS_corners is TRUE and if multiple measurements of each corner have been realized") } + if (nrow(unique(rel_coord)) != 4) { + stop( paste(nrow(unique(rel_coord)),"unique corners are detected in rel_coord instead of 4")) + } if (!is.null(corner_ID) && length(corner_ID) != nrow(rel_coord)) { stop("corner_ID must be the same length as the number of rows in rel_coord") } + if (!is.null(corner_ID) && length(unique(corner_ID)) !=4) { + stop(paste(length((unique(corner_ID))),"unique corners are detected in corner_ID instead of 4")) + } # function ------------------------------------------------------------------- diff --git a/tests/testthat/test-check_coord_plot.R b/tests/testthat/test-check_coord_plot.R index 93e2771..87b6822 100644 --- a/tests/testthat/test-check_coord_plot.R +++ b/tests/testthat/test-check_coord_plot.R @@ -38,7 +38,6 @@ test_that("check_plot_coord error", { expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=T,corner_ID = wrong_corner_ID),"corner_ID instead of 4") }) - test_that("check_plot_coord outputs and outliers", { # with max dist equal 10 expect_warning( @@ -164,10 +163,8 @@ test_that("check_plot_coord, tree coordinates", { expect_equivalent(res$tree_proj_coord , sweep(tree_df, 2, c(200000,9000000), FUN = "+")) tree_df$AGB <- c(20,25,30) expect_equal(res, check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = T, rm_outliers = T, corner_ID = corner_ID, draw_plot = F, tree_df = tree_df, tree_coords = c("X","Y"))) - }) - test_that("check_plot_coord, splashed corners", { proj_coord <- data.frame( X = rep(c(0,0,100,100), e=4), @@ -184,8 +181,3 @@ test_that("check_plot_coord, splashed corners", { corner_ID <- rep(c("SW","NW","SE","NE"),e=4) expect_error(suppressWarnings(check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = T, corner_ID = corner_ID, rm_outliers = T, draw_plot = F, max_dist = 5))) }) - - - - - From da33b448a604f8238e7e9b65be47a4ff8250154e Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Fri, 29 Nov 2024 12:03:08 +0100 Subject: [PATCH 14/16] changed donttest to dontrun for correctTaxo --- man/BIOMASS-package.Rd | 7 +------ man/correctTaxo.Rd | 2 +- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/man/BIOMASS-package.Rd b/man/BIOMASS-package.Rd index 2515063..6394c7b 100644 --- a/man/BIOMASS-package.Rd +++ b/man/BIOMASS-package.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/pkgname.R +% Please edit documentation in R/BIOMASS-package.R \docType{package} \name{BIOMASS-package} \alias{BIOMASS} @@ -8,11 +8,6 @@ \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, (3) construct height-diameter models, (4) manage tree and plot coordinates, (5) estimate the aboveground biomass/carbon at the stand level with associated uncertainty. To cite 'BIOMASS', please use citation("BIOMASS"). See more in the article of Réjou-Méchain et al. (2017) \doi{10.1111/2041-210X.12753}. } -\references{ -Réjou-Méchain M., Tanguy A., Piponiot C., Chave J., Hérault B. (2017). BIOMASS : -An R Package for estimating above-ground biomass and its uncertainty in tropical forests. -Methods in Ecology and Evolution, 8(9), 1163-1167. -} \seealso{ Useful links: \itemize{ diff --git a/man/correctTaxo.Rd b/man/correctTaxo.Rd index 8e22d4d..bd12d60 100644 --- a/man/correctTaxo.Rd +++ b/man/correctTaxo.Rd @@ -57,7 +57,7 @@ Cache path discovery protocol } \examples{ -\donttest{ +\dontrun{ correctTaxo(genus = "Astrocarium", species = "standleanum") correctTaxo(genus = "Astrocarium standleanum") } From 5e7ccff0008a2158c31f4aa76a83772588d5f327 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Fri, 29 Nov 2024 15:14:41 +0100 Subject: [PATCH 15/16] Added subplot_summary tests, including ggplot tests based on snaps (vdiffr::expect_doppelganger) --- DESCRIPTION | 1 + .../subplot_summary/disp-subplot-summary.svg | 86 +++++++++++++++++++ tests/testthat/test-subplot_summary.R | 23 +++++ 3 files changed, 110 insertions(+) create mode 100644 tests/testthat/_snaps/subplot_summary/disp-subplot-summary.svg diff --git a/DESCRIPTION b/DESCRIPTION index 85ab2db..83d4e4a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -48,6 +48,7 @@ Suggests: rmarkdown, prettydoc, testthat, + vdiffr, curl, geodata, httr2, diff --git a/tests/testthat/_snaps/subplot_summary/disp-subplot-summary.svg b/tests/testthat/_snaps/subplot_summary/disp-subplot-summary.svg new file mode 100644 index 0000000..0aed408 --- /dev/null +++ b/tests/testthat/_snaps/subplot_summary/disp-subplot-summary.svg @@ -0,0 +1,86 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +200 +250 +300 +350 +400 +450 +100 +150 +200 +250 +300 +350 +metric_summary_per_ha + + + + + + + + + + + +100 +110 +120 +130 +140 +Summary of metric per ha + + diff --git a/tests/testthat/test-subplot_summary.R b/tests/testthat/test-subplot_summary.R index a14a07e..2c460aa 100644 --- a/tests/testthat/test-subplot_summary.R +++ b/tests/testthat/test-subplot_summary.R @@ -57,4 +57,27 @@ test_that("subplot_summary", { # Test with quantile function res_quantile <- subplot_summary(subplots_less_tree, value = "metric", draw_plot = F, fun = quantile, probs=0.75) expect_equal(res_quantile$tree_summary[res_quantile$tree_summary$subplot_id=="subplot_3_2","metric_summary"],16.255414,tol=1e-5) + + # Test multiple plots + multiple_rel_coord <- rbind(rel_coord , rel_coord) + multiple_proj_coord <- rbind(proj_coord,proj_coord) + multiple_rel_coord$plotID <- rep(c("plot1","plot2"),e=4) + multiple_tree_df <- rbind(tree_df,tree_df) + multiple_tree_df$plot_id <- rep(c("plot1","plot2"),e=nrow(tree_df)) + multiple_subplots <- divide_plot(rel_coord = multiple_rel_coord[,c("x_rel","y_rel")], grid_size = 50, proj_coord = multiple_proj_coord, corner_plot_ID = multiple_rel_coord$plotID, + tree_df = multiple_tree_df , tree_coords = c("x_tree","y_tree"), tree_plot_ID = multiple_tree_df$plot_id) + res_multiple <- subplot_summary(multiple_subplots, value = "metric", draw_plot = F) + expect_equivalent(res_multiple$tree_summary[1:15,c("metric_summary","metric_summary_per_ha")] , res_less_tree$tree_summary[c("metric_summary","metric_summary_per_ha")]) + expect_equivalent(res_multiple$tree_summary$subplot_id[1:5],c("plot1_0_0","plot1_0_1","plot1_0_2","plot1_0_3","plot1_1_0")) +}) + + +test_that("ggplot_subplot_summary", { + set.seed(52) + 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(210, 383, 110, 283), y_proj = c(210, 310, 383, 483)) + subplots <- divide_plot(rel_coord, grid_size = 50) + tree_df <- data.frame(x_tree = runif(50,0,200), y_tree = runif(50,0,200), metric = rnorm(50,10,5)) + subplots <- divide_plot(rel_coord, proj_coord = proj_coord, grid_size = 100, tree_df = tree_df, tree_coords = c("x_tree","y_tree")) + vdiffr::expect_doppelganger("disp-subplot-summary", subplot_summary(subplots, value = "metric", draw_plot = T)$plot_design) }) From 18cec9596eef08ff8ee1be6c8d98af096242c02f Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Fri, 29 Nov 2024 16:30:13 +0100 Subject: [PATCH 16/16] Deal with plot names containing an underscore Invisible the ggplot return --- R/subplot_summary.R | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/R/subplot_summary.R b/R/subplot_summary.R index 435b33d..e4483dc 100644 --- a/R/subplot_summary.R +++ b/R/subplot_summary.R @@ -63,7 +63,9 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, if(any(duplicated(tree_summary$subplot_id))) { stop("the function supplied using `fun` must return a single value") } - tree_summary[,plot_id:=sapply(strsplit(tree_summary$subplot_id,"_"),function(x)x[[1]])] #Add plot_id to be able to loop on it + tree_summary[, plot_id:= #Add plot_id to be able to loop on it + sapply(strsplit(tree_summary$subplot_id,"_"), + function(x) paste(x[-((length(x)-1) : length(x))],collapse="_"))] # instead of function(x) x[1] in case plot_id contains any "_" sf_polygons <- do.call(rbind,lapply(split(corner_dat, by = "subplot_id"), function(dat) { @@ -111,7 +113,7 @@ subplot_summary <- function(subplots, value = NULL, draw_plot = TRUE, fun = sum, } print(plot_design) - plot_design + return(invisible(plot_design)) }) names(plot_list) <- unique(tree_summary$plot_id)