From 8a93aeb89995add0e1d11e1fae586e193783bbda Mon Sep 17 00:00:00 2001 From: abailly Date: Fri, 15 Nov 2024 18:22:05 +0100 Subject: [PATCH 1/2] - Added check_plot_coord() which will replace checkPlotCoord() (added a deprecated warning in checkPlotCoord) - Replaced tree_coord argument by tree_df and tree_coords wich contains resp. the tree data-frame and the column names of the relative coordinates. --- R/checkPlotCoord.R | 118 ++++----- R/check_plot_coord.R | 334 +++++++++++++++++++++++++ tests/testthat/test-checkCoordPlot.R | 189 -------------- tests/testthat/test-check_coord_plot.R | 191 ++++++++++++++ 4 files changed, 580 insertions(+), 252 deletions(-) create mode 100644 R/check_plot_coord.R delete mode 100644 tests/testthat/test-checkCoordPlot.R create mode 100644 tests/testthat/test-check_coord_plot.R diff --git a/R/checkPlotCoord.R b/R/checkPlotCoord.R index 2a02190..e51154f 100644 --- a/R/checkPlotCoord.R +++ b/R/checkPlotCoord.R @@ -4,7 +4,7 @@ #' Quality check of plot corner and tree coordinates. #' #' @details -#' If trustGPScorners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the maxDist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. +#' If trustGPScorners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the maxDist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. Be aware that this projection only works if the plot, in the relative coordinates system, is rectangular (ie, has 4 right angles). #' #' If trustGPScorners is FALSE, corner coordinates in the projected coordinate system are calculated by a procrust analysis that preserves the shape and dimensions of the plot in the local coordinate system. Outlier corners are also identified sequentially and projected coordinates of the trees are calculated by applying the resulting procrust analysis. #' @@ -18,23 +18,23 @@ #' @param maxDist a numeric giving the maximum distance (in meters) above which GPS measurements should be considered outliers (default 15 m) #' @param rmOutliers a logical indicating if detected outliers are removed from the coordinate calculation #' @param drawPlot a logical indicating if the plot design should be displayed and returned -#' @param treeCoord a data frame containing at least the tree coordinates in the relative coordinates system (that of the field), with X and Y corresponding to the first and second columns respectively +#' @param treeCoord a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively #' #' @author Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY #' #' @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 #' -#' @importFrom data.table data.table := setnames +#' @importFrom data.table data.table := setnames %between% #' @importFrom sf st_multipoint st_polygon st_sfc -#' @importFrom ggplot2 ggplot aes geom_point geom_segment geom_polygon geom_text scale_shape_manual scale_color_manual ggtitle theme_minimal theme coord_equal arrow unit +#' @importFrom ggplot2 ggplot aes geom_point geom_segment geom_polygon geom_text scale_shape_manual scale_color_manual ggtitle theme_minimal theme coord_equal arrow unit element_blank #' #' @author Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY #' @@ -78,7 +78,10 @@ #' checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPScorners, cornerID=NULL, maxDist = 15, rmOutliers = TRUE, drawPlot = TRUE, treeCoord = NULL) { - # parameters verification ------------------------------------------------- + + warning("please, use check_plot_coord() instead of checkPlotCoord(). `checkPlotCoord' will be removed in the next version ") + + # Checking arguments ------------------------------------------------- if (is.null(longlat) && is.null(projCoord)) { stop("Give at least one set of coordinates: longlat or projCoord") @@ -143,17 +146,17 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc } if(trustGPScorners == TRUE) { - - if(nrow(projCoord)!= 4) { # if multiple measures of each corner, then do the mean of coordinates and look for outliers + + if(nrow(projCoord)!= 4) { # if multiple measures of each corner, then do the mean of coordinates and search for outliers cornerCoord <- data.table(cbind(projCoord[,1:2], relCoord[,1:2],cornerID=cornerID)) - setnames(cornerCoord, c("X","Y","Xrel","Yrel","cornerID")) + 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,22 +168,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 <- cornerCoord[ , c("Xmean","Ymean","Xrel","Yrel","cornerID")] + cornerCoord <- data.frame(cornerCoord[ , c("Xmean","Ymean","x_rel","y_rel","cornerID")]) cornerCoord <- unique(cornerCoord) - setnames(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.table(cbind(projCoord[,1:2], relCoord[,1:2])) - setnames(cornerCoord, c("X","Y","Xrel","Yrel")) + cornerCoord <- data.frame(cbind(projCoord[,1:2], relCoord[,1:2])) + colnames(cornerCoord) <- c("x_proj","y_proj","x_rel","y_rel") outliers <- data.frame() if(!is.null(cornerID)) cornerCoord$cornerID = cornerID } @@ -188,7 +190,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 @@ -200,7 +202,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 @@ -220,7 +222,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,47 +235,37 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc cornerProjCoord <- as.matrix(cornerRelCoord[,1:2]) %*% procrustRes$rotation cornerProjCoord <- sweep(cornerProjCoord, 2, procrustRes$translation, FUN = "+") - cornerCoord <- data.table(cbind(cornerProjCoord[,1:2], cornerRelCoord[,1:2])) - setnames(cornerCoord, c("X","Y","Xrel","Yrel")) - + cornerCoord <- data.frame(cbind(cornerProjCoord[,1:2], cornerRelCoord[,1:2])) + colnames(cornerCoord) <- c("x_proj","y_proj","x_rel","y_rel") + if(!is.null(cornerID)) cornerCoord$cornerID <- unique(cornerID) } # End trustGPScorners = "FALSE" - # Assign corner numbers in a clockwise direction ----------------------------- - m1 <- cornerCoord[ rank(Xrel) <= 2, ] - tmp1 <- m1[rank(Yrel),] - m2 <- cornerCoord[ rank(Xrel) > 2, ] - tmp2 <- m2[rank(Yrel),] - cornerCoord <- rbind(tmp1,tmp2) - cornerCoord$cornerNum <- c(1,2,4,3) - - # if(! is.null(originalCorner)) { - # # Shift the corner numbers to have corner 1 on the origin - # shift <- 5 - which(cornerCoord$cornerID == originalCorner) - # cornerCoord$cornerNum <- (cornerCoord$cornerNum + shift) %% 4 - # cornerCoord$cornerNum[cornerCoord$cornerNum == 0] <- 4 - # } - - cornerCoord <- cornerCoord[ order(cornerNum), ] + # Sort cornerCoord rows in a counter-clockwise direction --------------------- + centroid <- colMeans(cornerCoord[,c("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$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 <- bilinearInterpolation(relCoord = treeCoord[,1:2], - cornerCoord = cornerCoord[,c("X","Y","cornerNum")], - dimX = diff(range(cornerCoord$Xrel)), - dimY = diff(range(cornerCoord$Yrel)) - ) - colnames(treeProjCoord) <- c("X","Y") + treeProjCoord <- bilinear_interpolation(coord = treeCoord[,1:2], + 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.table(sweep(treeProjCoord, 2, procrustRes$translation, FUN = "+")) - colnames(treeProjCoord) <- c("X","Y") + treeProjCoord <- data.frame(sweep(treeProjCoord, 2, procrustRes$translation, FUN = "+")) + colnames(treeProjCoord) <- c("x_proj","y_proj") } } @@ -281,9 +273,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)) { @@ -291,20 +283,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())+ @@ -322,7 +314,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc } output <- list( - cornerCoord = data.frame(cornerCoord), + cornerCoord = cornerCoord, polygon = cornerPolygon, outliers = outliers ) if (drawPlot) { @@ -332,7 +324,7 @@ checkPlotCoord <- function(projCoord = NULL, longlat = NULL, relCoord, trustGPSc output$codeUTM <- codeUTM } if (!is.null(treeCoord)) { - output$treeProjCoord <- data.frame(treeProjCoord[,c("X","Y")]) + output$treeProjCoord <- treeProjCoord[,c("x_proj","y_proj")] } return(output) -} +} \ No newline at end of file diff --git a/R/check_plot_coord.R b/R/check_plot_coord.R new file mode 100644 index 0000000..04745c3 --- /dev/null +++ b/R/check_plot_coord.R @@ -0,0 +1,334 @@ +#' Check coordinates of plot corners and trees +#' +#' @description +#' Quality check of plot corner and tree coordinates. +#' +#' @details +#' If trust_GPS_corners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the max_dist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. Be aware that this projection only works if the plot, in the relative coordinates system, is rectangular (ie, has 4 right angles). +#' +#' If trust_GPS_corners is FALSE, corner coordinates in the projected coordinate system are calculated by a procrust analysis that preserves the shape and dimensions of the plot in the local coordinate system. Outlier corners are also identified sequentially and projected coordinates of the trees are calculated by applying the resulting procrust analysis. +#' +#' If longlat is supplied instead of proj_coord, the function will first convert the long/lat coordinates into UTM coordinates. An error may result if the parcel is located right between two UTM zones. In this case, the user has to convert himself his long/lat coordinates into any projected coordinates which have the same dimension than his local coordinates (in meters most of the time). +#' +#' @param proj_coord (optional, if longlat is not supplied) a data frame containing the plot corner coordinates in a projected coordinate system, with X and Y corresponding to the first and second columns respectively. +#' @param longlat (optional, if proj_coord is not supplied) a data frame containing the plot corner coordinates in longitude and latitude, with longitude and latitude corresponding to the first and second columns respectively +#' @param rel_coord a data frame containing the plot corner coordinates in the relative coordinates system (that of the field), with X and Y corresponding to the first and second columns respectively, and with the same row order than longlat or proj_coord +#' @param trust_GPS_corners a logical indicating whether or not you trust the GPS coordinates of the plot's corners. See details. +#' @param corner_ID a vector indicating the ID of the corners (e.g c("SE","SE",...)) in the case you have multiple measurements for each corner +#' @param max_dist a numeric giving the maximum distance (in meters) above which GPS measurements should be considered outliers (default 15 m) +#' @param rm_outliers a logical indicating if detected outliers are removed from the coordinate calculation +#' @param draw_plot a logical indicating if the plot design should be displayed and returned +#' @param tree_df a data frame containing the relative coordinates (field/local coordinates) of the trees and optional other tree metrics +#' @param tree_coords a character vector of length 2 containing the column names of the relative coordinates in tree_df +#' +#' @author Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY +#' +#' @return Returns a list including : +#' - `corner_coord`: a data frame containing the projected coordinates (x_proj and y_proj), the relative coordinates (x_rel and y_rel), and the ID (if corner_ID is supplied) of the 4 corners of the plot +#' - `polygon`: a spatial polygon +#' - `outliers`: a data frame containing the projected coordinates, the ID (if corner_ID is supplied) and the row number of GPS points considered outliers +#' - `plot_design`: if `draw_plot` is TRUE, a ggplot object corresponding to the design of the plot +#' - `UTM_code`: if `longlat` is supplied, a character containing the UTM code of the GPS coordinates +#' - `tree_proj_coord`: if `tree_df` is supplied, a data frame containing the coordinates of the trees in the projected coordinate system (x_proj and y_proj) +#' +#' @export +#' +#' @importFrom data.table data.table := setnames %between% +#' @importFrom sf st_multipoint st_polygon st_sfc +#' @importFrom ggplot2 ggplot aes geom_point geom_segment geom_polygon geom_text scale_shape_manual scale_color_manual ggtitle theme_minimal theme coord_equal arrow unit element_blank +#' +#' @author Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY +#' +#' @examples +#' proj_coord <- data.frame( +#' X = c( +#' runif(5, min = 9, max = 11), runif(5, min = 8, max = 12), +#' runif(5, min = 80, max = 120), runif(5, min = 90, max = 110) +#' ), +#' Y = c( +#' runif(5, min = 9, max = 11), runif(5, min = 80, max = 120), +#' runif(5, min = 8, max = 12), runif(5, min = 90, max = 110) +#' ) +#' ) +#' proj_coord <- proj_coord + 1000 +#' rel_coord <- data.frame( +#' X = c(rep(0, 10), rep(100, 10)), +#' Y = c(rep(c(rep(0, 5), rep(100, 5)), 2)) +#' ) +#' corner_ID <- rep(c("SW","NW","SE","NE"),e=5) +#' +#' aa <- check_plot_coord( +#' proj_coord = proj_coord, rel_coord = rel_coord, +#' trust_GPS_corners = TRUE, corner_ID = corner_ID, +#' rm_outliers = FALSE , draw_plot = FALSE +#' ) +#' bb <- check_plot_coord( +#' proj_coord = proj_coord, rel_coord = rel_coord, +#' trust_GPS_corners = FALSE, +#' rm_outliers = TRUE, max_dist = 10, +#' draw_plot = FALSE +#' ) +#' \donttest{ +#' check_plot_coord( +#' proj_coord = proj_coord, rel_coord = rel_coord, +#' trust_GPS_corners = TRUE, corner_ID = corner_ID, +#' rm_outliers = TRUE , max_dist = 10, +#' 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 ------------------------------------------------- + + if (is.null(longlat) && is.null(proj_coord)) { + stop("Give at least one set of coordinates: longlat or proj_coord") + } + if (!is.null(longlat) & !is.null(proj_coord)) { + stop("Give only one set of coordinates: longlat or proj_coord") + } + if((!is.matrix(rel_coord) & !is.data.frame(rel_coord))){ + stop("rel_coord must be a matrix or a data frame") + } + if(!is.null(longlat) && (!is.matrix(longlat) & !is.data.frame(longlat))){ + stop("longlat must be a matrix or a data frame") + } + if(!is.null(proj_coord) && (!is.matrix(proj_coord) & !is.data.frame(proj_coord))){ + stop("proj_coord must be a matrix or a data frame") + } + if(!is.null(tree_df) && (!is.matrix(tree_df) & !is.data.frame(tree_df))){ + stop("tree_df must be a matrix or a data frame") + } + if (!is.null(tree_df) && is.null(tree_coords)) { + stop("You must supply the column names corresponding to the relative coordinates of tree_df using the argument `tree_coords`") + } + if (!is.null(tree_df) && !any(tree_coords %in% names(tree_df))) { + stop("column names supplied by tree_coords are not found in tree_df") + } + if (is.null(trust_GPS_corners) | !is.logical(trust_GPS_corners)) { + stop("The trust_GPS_corners argument must be TRUE or FALSE") + } + if (length(max_dist) != 1) { + stop("The max_dist argument must be of length 1") + } + if ((!is.null(longlat) && any(dim(longlat) != dim(rel_coord))) || (!is.null(proj_coord) && any(dim(proj_coord) != dim(rel_coord)))) { + stop("GPS and relative coordinates are not of the same dimension") + } + if(sum(is.na(rel_coord))!= 0) { + stop("Missing values are detected in rel_coord. Please remove them and call the function again") + } + if(!is.null(longlat) && sum(is.na(longlat))!= 0) { + stop("Missing values are detected in longlat. Please remove them and call the function again") + } + if(!is.null(proj_coord) && sum(is.na(proj_coord))!= 0) { + stop("Missing values are detected in proj_coord. Please remove them and call the function again") + } + 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 ------------------------------------------------------------------- + + # Transform the geographic coordinates into UTM coordinates + if (!is.null(longlat)) { + proj_coord <- latlong2UTM(longlat) + UTM_code <- unique(proj_coord[, "codeUTM"]) + if(length(UTM_code)>1) { + stop("More than one UTM zone are detected. This may be due to an error in the long/lat coordinates, or if the parcel is located right between two UTM zones. In this case, please convert yourself your long/lat coordinates into any projected coordinates which have the same dimension than your local coordinates") + } + proj_coord <- proj_coord[, c("X", "Y")] + } + + if(trust_GPS_corners == TRUE) { + + if(nrow(proj_coord)!= 4) { # if multiple measures of each corner, then do the mean of coordinates and search for outliers + + corner_coord <- data.table(cbind(proj_coord[,1:2], rel_coord[,1:2],corner_ID=corner_ID)) + setnames(corner_coord, c("x_proj","y_proj","x_rel","y_rel","corner_ID")) + if (any(table(corner_coord$corner_ID) < 5)) { + warning("At least one corner has less than 5 measurements. We suggest using the argument trust_GPS_corners = FALSE") + } + corner_coord[ , c("Xmean","Ymean","row_number") := list(mean(x_proj),mean(y_proj),.I) , by=corner_ID] + corner_coord[ , outlier := ifelse( sqrt((x_proj - Xmean)^2 + (y_proj-Ymean)^2) > max_dist, T, F)] + outliers <- corner_coord[outlier==T , c("x_proj","y_proj","corner_ID","row_number")] + + if (any(corner_coord[,sum(!outlier),by=corner_ID][[2]]==0)) { + stop("All measurements for at least one corner are considered as outliers.\n + This may be because some coordinates have very large error associated.\n + Try to remove these very large error or reconsider the max_dist parameter by increasing the distance") + } + + if(rm_outliers & nrow(outliers)!=0) { + corner_coord <- corner_coord[outlier==F ,] + exist_outliers <- T + while(exist_outliers) { + corner_coord[ , outlier := ifelse( sqrt((x_proj - Xmean)^2 + (y_proj-Ymean)^2) > max_dist, T, F)] + outliers <- rbind(outliers , corner_coord[outlier==T , c("x_proj","y_proj","corner_ID","row_number")]) + corner_coord <- corner_coord[outlier==F ,] + if(nrow(corner_coord)+nrow(outliers) == nrow(proj_coord)) exist_outliers <- F + corner_coord[ , c("Xmean","Ymean") := list(mean(x_proj),mean(y_proj)) , by=corner_ID] + } + } + + corner_coord <- data.frame(corner_coord[ , c("Xmean","Ymean","x_rel","y_rel","corner_ID")]) + corner_coord <- unique(corner_coord) + colnames(corner_coord) <- c("x_proj", "y_proj","x_rel","y_rel","corner_ID") + + } else { # if exactly 1 measures for each corner + corner_coord <- data.frame(cbind(proj_coord[,1:2], rel_coord[,1:2])) + colnames(corner_coord) <- c("x_proj","y_proj","x_rel","y_rel") + outliers <- data.frame() + if(!is.null(corner_ID)) corner_coord$corner_ID = corner_ID + } + + } else { # trust_GPS_corners = "FALSE" + + corner_coord <- data.table(cbind(proj_coord[,1:2], rel_coord[,1:2])) + setnames(corner_coord, c("x_proj","y_proj","x_rel","y_rel")) + corner_coord[, row_number := .I] + + # Transformation of rel_coord to projected coordinates by a procrust analyses + procrust_res <- procrust(proj_coord, rel_coord) + procrust_coord <- as.matrix(rel_coord) %*% procrust_res$rotation + procrust_coord <- sweep(procrust_coord, 2, procrust_res$translation, FUN = "+") + + # Calculate the distances between the GNSS measurements and the procrust_coord + corner_coord[ , outlier := ifelse( + sqrt((procrust_coord[, 1] - proj_coord[, 1])^2 + (procrust_coord[, 2] - proj_coord[, 2])^2) > max_dist, T, F)] + + outliers <- corner_coord[outlier==T , c("x_proj","y_proj","row_number")] + + if (length(outliers)==nrow(proj_coord)){ + stop("All coordinates points are considered as outliers at the first stage.\n + This may be because some coordinates have very large error associated.\n + Try to remove these very large error or reconsider the max_dist parameter by increasing the distance") + } + + # retransform the rel_coord without the outliers + if (rm_outliers & nrow(outliers)>0) { + corner_coord <- corner_coord[outlier==F ,] + refineCoord <- TRUE + while(refineCoord){ + procrust_res <- procrust(proj_coord[-outliers$row_number, ], rel_coord[-outliers$row_number,]) + procrust_coord <- as.matrix(rel_coord[-outliers$row_number,]) %*% procrust_res$rotation + procrust_coord <- sweep(procrust_coord, 2, procrust_res$translation, FUN = "+") + corner_coord[ , outlier := ifelse( + sqrt((procrust_coord[, 1] - proj_coord[-outliers$row_number, 1])^2 + + (procrust_coord[, 2] - proj_coord[-outliers$row_number, 2])^2) > max_dist, T, F)] + + outliers <- rbind(outliers , corner_coord[outlier==T , c("x_proj","y_proj","row_number")]) + corner_coord <- corner_coord[outlier==F ,] + if(nrow(corner_coord)+nrow(outliers) == nrow(proj_coord)) refineCoord <- F + } + } + + # Create the data.table of corners to return the projected coordinate of the plot's corner + corner_rel_coord <- unique(rel_coord) + + # Project corner_rel_coord + corner_proj_coord <- as.matrix(corner_rel_coord[,1:2]) %*% procrust_res$rotation + corner_proj_coord <- sweep(corner_proj_coord, 2, procrust_res$translation, FUN = "+") + + corner_coord <- data.frame(cbind(corner_proj_coord[,1:2], corner_rel_coord[,1:2])) + colnames(corner_coord) <- c("x_proj","y_proj","x_rel","y_rel") + + if(!is.null(corner_ID)) corner_coord$corner_ID <- unique(corner_ID) + } # End trust_GPS_corners = "FALSE" + + # Sort corner_coord rows in a counter-clockwise direction --------------------- + centroid <- colMeans(corner_coord[,c("x_rel","y_rel")]) + angles <- base::atan2(corner_coord[["y_rel"]] - centroid[2], corner_coord[["x_rel"]] - centroid[1]) + corner_coord <- corner_coord[order(angles), ] + + # Create a polygon ----------------------------------------------------------- + corner_polygon <- st_multipoint(as.matrix(rbind(corner_coord[,c("x_proj","y_proj")], corner_coord[1,c("x_proj","y_proj")]))) + corner_polygon <- st_polygon(x = list(corner_polygon), dim = "XY") + corner_polygon <- st_sfc(list(corner_polygon)) + + # Calculate projected tree coordinates from relative tree coordinates -------- + if(!is.null(tree_df)) { + if(any(! tree_df[,tree_coords[1]] %between% range(corner_coord$x_rel)) | any(! tree_df[,tree_coords[2]] %between% range(corner_coord$y_rel))) { + warning("Be careful, one or more trees are not inside the plot defined by rel_coord") + } + if(trust_GPS_corners) { + tree_proj_coord <- bilinear_interpolation(coord = tree_df[,tree_coords], + from_corner_coord = corner_coord[,c("x_rel","y_rel")], + to_corner_coord = corner_coord[,c("x_proj","y_proj")], + ordered_corner = T) + colnames(tree_proj_coord) <- c("x_proj","y_proj") + } else { + tree_proj_coord <- as.matrix(tree_df[,tree_coords]) %*% procrust_res$rotation + tree_proj_coord <- data.frame(sweep(tree_proj_coord, 2, procrust_res$translation, FUN = "+")) + colnames(tree_proj_coord) <- c("x_proj","y_proj") + } + } + + # draw plot ------------------------------------------------------------------ + + if (draw_plot) { + proj_coord_plot <- proj_coord[,c(1,2)] + setnames(proj_coord_plot , c("x_proj","y_proj")) + proj_coord_plot$whatpoint <- "GPS measurements" + corner_coord_plot <- corner_coord[,c("x_proj","y_proj")] + corner_coord_plot$whatpoint <- "Reference corners" + proj_coord_plot <- rbind(proj_coord_plot,corner_coord_plot) + if(!is.null(tree_df)) { + tree_proj_coord$whatpoint <- "Trees" + proj_coord_plot <- rbind(proj_coord_plot,tree_proj_coord) + } + proj_coord_plot$whatpoint <- factor(proj_coord_plot$whatpoint, levels = c("GPS measurements","Outliers (discarded)","Reference corners","Trees")) + arrow_plot <- data.frame(X = rep(corner_coord$x_proj[1]), Y =rep(corner_coord$y_proj[1]), + Xend = c(corner_coord$x_proj[1]+(corner_coord$x_proj[4]-corner_coord$x_proj[1])/4,corner_coord$x_proj[1]+(corner_coord$x_proj[2]-corner_coord$x_proj[1])/4), + Yend = c(corner_coord$y_proj[1]+(corner_coord$y_proj[4]-corner_coord$y_proj[1])/4,corner_coord$y_proj[1]+(corner_coord$y_proj[2]-corner_coord$y_proj[1])/4)) + + if(exists("outliers") && nrow(outliers)!=0) { + proj_coord_plot$whatpoint[outliers$row_number] <- "Outliers (discarded)" + } + plot_design <- ggplot2::ggplot(data = proj_coord_plot) + + 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 = corner_polygon[[1]][[1]][,] , mapping = aes(x=x_proj,y=y_proj), colour="black",fill=NA,linewidth=1.2)+ + geom_segment(data=arrow_plot, mapping=aes(x=X,y=Y,xend=Xend,yend=Yend), arrow=arrow(length=unit(0.4,"cm")),linewidth=1, col="blue") + + geom_text(data = arrow_plot, 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())+ + coord_equal() + print(plot_design) + } + + # return --------------------------------------------------------------------- + + if (nrow(outliers) != 0 & !rm_outliers) { + warning( + "Be carefull, you may have GNSS measurement outliers. \n", + "Removing them may improve the georeferencing of your plot (see the rm_outliers argument)." + ) + } + + output <- list( + corner_coord = corner_coord, + polygon = corner_polygon, outliers = outliers + ) + if (draw_plot) { + output$plot_design <- plot_design + } + if (!is.null(longlat)) { + output$UTM_code <- UTM_code + } + if (!is.null(tree_df)) { + output$tree_proj_coord <- tree_proj_coord[,c("x_proj","y_proj")] + } + return(output) +} \ No newline at end of file diff --git a/tests/testthat/test-checkCoordPlot.R b/tests/testthat/test-checkCoordPlot.R deleted file mode 100644 index 1cff25b..0000000 --- a/tests/testthat/test-checkCoordPlot.R +++ /dev/null @@ -1,189 +0,0 @@ -set.seed(2) - -projCoord <- data.frame( - X = c(runif(5, min = 9, max = 11), runif(5, min = 8, max = 12), runif(5, min = 80, max = 120), runif(5, min = 90, max = 110)), - Y = c(runif(5, min = 9, max = 11), runif(5, min = 80, max = 120), runif(5, min = 8, max = 12), runif(5, min = 90, max = 110)) -) -projCoord$X <- projCoord$X + 200000 -projCoord$Y <- projCoord$Y + 9000000 - -relCoord <- data.frame( - X = c(rep(0, 10), rep(100, 10)), - Y = c(rep(c(rep(0, 5), rep(100, 5)), 2)) -) -cornerID <- rep(c("SW","NW","SE","NE"),e=5) - -context("Check plot coordinates") -test_that("checkCoordPlot error", { - expect_error(checkPlotCoord(), "Give at least one set of coordinates") - expect_error(checkPlotCoord(longlat=projCoord, projCoord=projCoord),"Give only one set of coordinates") - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=c(0,0)),"relCoord must be a matrix or a data frame") - expect_error(checkPlotCoord(longlat=c(0,0), relCoord=relCoord),"longlat must be a matrix or a data frame") - expect_error(checkPlotCoord(projCoord=c(0,0), relCoord=relCoord),"projCoord must be a matrix or a data frame") - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=relCoord, trustGPScorners=F, treeCoord = c(1,1)),"treeCoord must be a matrix or a data frame") - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=relCoord, trustGPScorners=NULL),"The trustGPScorners argument must be TRUE or FALSE") - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=relCoord, trustGPScorners=T, maxDist=c(10,20) ),"The maxDist argument must be of length 1") - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=rbind(relCoord, c(40, 40)),trustGPScorners=T),"same dimension") - wrongRelCoord <- relCoord ; wrongRelCoord[1,1] <- NA - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=wrongRelCoord, trustGPScorners=T,cornerID = cornerID),"Missing values are detected in relCoord. Please remove them and call the function again") - wrongProjCoord <- projCoord ; wrongProjCoord[1,1] <- NA - expect_error(checkPlotCoord(projCoord=wrongProjCoord, relCoord=relCoord, trustGPScorners=T,cornerID = cornerID),"Missing values are detected in projCoord. Please remove them and call the function again") - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=relCoord, trustGPScorners=T),"The argument cornerID is needed if trustGPScorners is TRUE and if multiple measurements of each corner have been realized") - wrongRelCoord[1,1] <- 10 - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=wrongRelCoord, trustGPScorners=T,cornerID = cornerID),"relCoord instead of 4") - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=relCoord, trustGPScorners=T,cornerID = c("a","b")),"cornerID must be the same length as the number of rows in relCoord") - wrongCornerID <- cornerID ; wrongCornerID[1] <- "SN" - expect_error(checkPlotCoord(projCoord=projCoord, relCoord=relCoord, trustGPScorners=T,cornerID = wrongCornerID),"cornerID instead of 4") -}) - - -test_that("checkPlotCoord outputs and outliers", { - # with max dist equal 10 - expect_warning( - checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F, rmOutliers = F, cornerID = cornerID, drawPlot = F, maxDist = 10),"Be carefull, you may have GNSS measurement outliers" - ) - outputs <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F, rmOutliers = T, cornerID = cornerID, drawPlot = F, maxDist = 10) - expect_is(outputs, "list") - expect_length(outputs, 3) - expect_equal(names(outputs), c("cornerCoord", "polygon", "outliers")) - expect_is(outputs$cornerCoord, "data.frame") - expect_is(outputs$polygon, "sfc_POLYGON") - expect_is(outputs$outliers, "data.frame") - - expect_equal(dim(outputs$outliers), c(10,3)) - expect_equal(dim(outputs$cornerCoord), c(4, 6)) - - # with max dist equal 25 there isn't outliers anymore - expect_failure(expect_warning( - checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F, rmOutliers = F, cornerID = cornerID, drawPlot = F, maxDist = 25), - "Be carefull, you may have GNSS measurement outliers" - )) - - # with rmOutliers = TRUE - expect_failure(expect_warning( - checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F, rmOutliers = T, cornerID = cornerID, drawPlot = F, maxDist = 25), - "Be carefull, you may have GNSS measurement outliers" - )) - outputs_2 <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F, rmOutliers = T, cornerID = cornerID, drawPlot = F, maxDist = 15) - - expect_failure(expect_equal(outputs$cornerCoord, outputs_2$cornerCoord)) - expect_failure(expect_equal(outputs$polygon, outputs_2$polygon)) -}) - - -test_that("checkPlotCoord in long lat", { - longlat <- as.data.frame(proj4::project(projCoord, - proj = "+proj=utm +zone=50 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs", - inverse = TRUE - )) - outputs_longlat <- checkPlotCoord(longlat = longlat, relCoord = relCoord, trustGPScorners = F,rmOutliers = T, cornerID = cornerID, drawPlot = F, maxDist = 10) - outputs_projCoord <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F,rmOutliers = T, cornerID = cornerID, drawPlot = F, maxDist = 10) - expect_is(outputs_longlat, "list") - expect_length(outputs_longlat, 4) - expect_equal(names(outputs_longlat), c("cornerCoord", "polygon", "outliers","codeUTM")) - expect_equal(outputs_longlat[1:3],outputs_projCoord) -}) - -test_that("checkPlotCoord, trustGPScorners", { - projCoord0 <- data.frame( - X = rep(c(10,10,100,100), e=7) + 200000, - Y = rep(c(10,100,100,10), e=7) + 9000000 - ) - projCoord <- projCoord0 - projCoord$X <- projCoord0$X + rep(sample(-3:3),4) - projCoord$Y <- projCoord0$Y + rep(sample(-3:3),4) - relCoord <- data.frame( - X = rep(c(0,0,95,95),e=7), - Y = rep(c(0,95,95,0),e=7) - ) - cornerID <- rep(c("SW","NW","SE","NE"),e=7) - - expect_warning( - checkPlotCoord(projCoord = projCoord[-(1:3),], relCoord = relCoord[-(1:3),], trustGPScorners = T, rmOutliers = T, cornerID = cornerID[-(1:3)], drawPlot = F),"At least one corner has less than 5 measurements. We suggest using the argument trustGPScorners = FALSE" - ) - res_trustGPScorners_T <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) - expect_equivalent(res_trustGPScorners_T$cornerCoord[,1:2] , unique(projCoord0)) - - res_trustGPScorners_F <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = F, rmOutliers = T, cornerID = cornerID, drawPlot = F) - expect_equal(res_trustGPScorners_T$cornerCoord[,3:6] , res_trustGPScorners_F$cornerCoord[,3:6]) - expect_equal(res_trustGPScorners_T$cornerCoord[,1:2], round(res_trustGPScorners_F$cornerCoord[,1:2],digits = -1)) - - # Test when there's only 4 measures and trustGPScorners = F - expect_equal(res_trustGPScorners_F , checkPlotCoord(projCoord = projCoord0[c(1,8,15,22),], relCoord = relCoord[c(1,8,15,22),], trustGPScorners = F, rmOutliers = T, drawPlot = F, cornerID = cornerID[c(1,8,15,22)])) -}) - -test_that("checkPlotCoord, origin corner", { - projCoord <- data.frame( - X = rep(c(0,0,100,100), e=7) + 200000 + rep(sample(-3:3),4), - Y = rep(c(0,100,100,0), e=7) + 9000000 + rep(sample(-3:3),4) - ) - cornerID <- rep(c("SW","NW","SE","NE"),e=7) - - # Test when the origin is not the South-East corner - relCoord_NE <- data.frame( - X = rep(c(100,100,0,0),e=7), - Y = rep(c(100,0,0,100),e=7) - ) - relCoord_SW <- data.frame( - X = rep(c(0,0,100,100),e=7), - Y = rep(c(0,100,100,0),e=7) - ) - res_NE <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord_NE, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) - res_SW <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord_SW, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) - - expect_equal(res_SW$cornerCoord[c("Xrel","Yrel","cornerNum")] , res_NE$cornerCoord[c("Xrel","Yrel","cornerNum")]) - expect_equivalent(res_SW$cornerCoord[c("X","Y","cornerID")] , res_NE$cornerCoord[c(3,4,1,2),c("X","Y","cornerID")]) - - # Test when the origin of the relative coordinates is not(0;0) - relCoord_SW <- relCoord_SW + 200 - res_SW_200 <- checkPlotCoord(projCoord = projCoord, relCoord = relCoord_SW, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F) - expect_equal(res_SW$cornerCoord[c("X","Y","cornerID","cornerNum")] , res_SW_200$cornerCoord[c("X","Y","cornerID","cornerNum")]) - expect_equal(res_SW$cornerCoord[c("Xrel","Yrel")]+200 , res_SW_200$cornerCoord[c("Xrel","Yrel")]) -}) - - -test_that("checkPlotCoord, tree coordinates", { - projCoord <- data.frame( - X = rep(c(0,0,100,100), e=7) + 200000 + rep(sample(-3:3),4), - Y = rep(c(0,100,100,0), e=7) + 9000000 + rep(sample(-3:3),4) - ) - relCoord <- data.frame( - X = rep(c(0,0,100,100),e=7), - Y = rep(c(0,100,100,0),e=7) - ) - cornerID <- rep(c("SW","NW","SE","NE"),e=7) - treeCoord <- data.frame( - X = c(10,20,30), - Y = 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") - expect_equal(dim(res$treeProjCoord), c(nrow(treeCoord),2)) - expect_equal(res$treeProjCoord , sweep(treeCoord, 2, c(200000,9000000), FUN = "+")) - treeCoord$AGB <- c(20,25,30) - expect_equal(res, checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = T, rmOutliers = T, cornerID = cornerID, drawPlot = F, treeCoord = treeCoord)) - -}) - - -test_that("checkPlotCoord, splashed corners", { - projCoord <- data.frame( - X = rep(c(0,0,100,100), e=4), - Y = rep(c(0,100,100,0), e=4) - ) - projCoord[1:4,] <- matrix(c(25,0,-25,0,0,25,0,-25), ncol=2, byrow = F) - projCoord$X = projCoord$X + 200000 - projCoord$Y = projCoord$Y + 9000000 - - relCoord <- data.frame( - X = rep(c(0,0,100,100),e=4), - Y = rep(c(0,100,100,0),e=4) - ) - cornerID <- rep(c("SW","NW","SE","NE"),e=4) - expect_error(suppressWarnings(checkPlotCoord(projCoord = projCoord, relCoord = relCoord, trustGPScorners = T, cornerID = cornerID, rmOutliers = T, drawPlot = F, maxDist = 5))) -}) - - - - - diff --git a/tests/testthat/test-check_coord_plot.R b/tests/testthat/test-check_coord_plot.R new file mode 100644 index 0000000..0996728 --- /dev/null +++ b/tests/testthat/test-check_coord_plot.R @@ -0,0 +1,191 @@ +set.seed(2) + +proj_coord <- data.frame( + X = c(runif(5, min = 9, max = 11), runif(5, min = 8, max = 12), runif(5, min = 80, max = 120), runif(5, min = 90, max = 110)), + Y = c(runif(5, min = 9, max = 11), runif(5, min = 80, max = 120), runif(5, min = 8, max = 12), runif(5, min = 90, max = 110)) +) +proj_coord$X <- proj_coord$X + 200000 +proj_coord$Y <- proj_coord$Y + 9000000 + +rel_coord <- data.frame( + X = c(rep(0, 10), rep(100, 10)), + Y = c(rep(c(rep(0, 5), rep(100, 5)), 2)) +) +corner_ID <- rep(c("SW","NW","SE","NE"),e=5) + +context("Check plot coordinates") +test_that("check_plot_coord error", { + expect_error(check_plot_coord(), "Give at least one set of coordinates") + expect_error(check_plot_coord(longlat=proj_coord, proj_coord=proj_coord),"Give only one set of coordinates") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=c(0,0)),"rel_coord must be a matrix or a data frame") + expect_error(check_plot_coord(longlat=c(0,0), rel_coord=rel_coord),"longlat must be a matrix or a data frame") + expect_error(check_plot_coord(proj_coord=c(0,0), rel_coord=rel_coord),"proj_coord must be a matrix or a data frame") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=F, tree_df = c(1,1)),"tree_df must be a matrix or a data frame") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=F, tree_df = data.frame(x=1,y=1)),"You must supply the column names corresponding to the relative coordinates of tree_df using the argument `tree_coords`") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=F,tree_df = data.frame(x=1,y=1), tree_coords = c("a","b")),"column names supplied by tree_coords are not found in tree_df") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=NULL),"The trust_GPS_corners argument must be TRUE or FALSE") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=T, max_dist=c(10,20) ),"The max_dist argument must be of length 1") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rbind(rel_coord, c(40, 40)),trust_GPS_corners=T),"same dimension") + wrong_rel_coord <- rel_coord ; wrong_rel_coord[1,1] <- NA + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=wrong_rel_coord, trust_GPS_corners=T,corner_ID = corner_ID),"Missing values are detected in rel_coord. Please remove them and call the function again") + wrong_proj_coord <- proj_coord ; wrong_proj_coord[1,1] <- NA + expect_error(check_plot_coord(proj_coord=wrong_proj_coord, rel_coord=rel_coord, trust_GPS_corners=T,corner_ID = corner_ID),"Missing values are detected in proj_coord. Please remove them and call the function again") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=T),"The argument corner_ID is needed if trust_GPS_corners is TRUE and if multiple measurements of each corner have been realized") + wrong_rel_coord[1,1] <- 10 + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=wrong_rel_coord, trust_GPS_corners=T,corner_ID = corner_ID),"rel_coord instead of 4") + expect_error(check_plot_coord(proj_coord=proj_coord, rel_coord=rel_coord, trust_GPS_corners=T,corner_ID = c("a","b")),"corner_ID must be the same length as the number of rows in rel_coord") + wrong_corner_ID <- corner_ID ; wrong_corner_ID[1] <- "SN" + 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( + check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = F, corner_ID = corner_ID, draw_plot = F, max_dist = 10),"Be carefull, you may have GNSS measurement outliers" + ) + outputs <- check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = T, corner_ID = corner_ID, draw_plot = F, max_dist = 10) + expect_is(outputs, "list") + expect_length(outputs, 3) + expect_equal(names(outputs), c("corner_coord", "polygon", "outliers")) + expect_is(outputs$corner_coord, "data.frame") + expect_is(outputs$polygon, "sfc_POLYGON") + expect_is(outputs$outliers, "data.frame") + + expect_equal(dim(outputs$outliers), c(10,3)) + expect_equal(dim(outputs$corner_coord), c(4, 5)) + + # with max dist equal 25 there isn't outliers anymore + expect_failure(expect_warning( + check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = F, corner_ID = corner_ID, draw_plot = F, max_dist = 25), + "Be carefull, you may have GNSS measurement outliers" + )) + + # with rm_outliers = TRUE + expect_failure(expect_warning( + check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = T, corner_ID = corner_ID, draw_plot = F, max_dist = 25), + "Be carefull, you may have GNSS measurement outliers" + )) + outputs_2 <- check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = T, corner_ID = corner_ID, draw_plot = F, max_dist = 15) + + expect_failure(expect_equal(outputs$corner_coord, outputs_2$corner_coord)) + expect_failure(expect_equal(outputs$polygon, outputs_2$polygon)) +}) + + +test_that("check_plot_coord in long lat", { + longlat <- as.data.frame(proj4::project(proj_coord, + proj = "+proj=utm +zone=50 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs", + inverse = TRUE + )) + outputs_longlat <- check_plot_coord(longlat = longlat, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = T, corner_ID = corner_ID, draw_plot = F, max_dist = 10) + outputs_proj_coord <- check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = T, corner_ID = corner_ID, draw_plot = F, max_dist = 10) + expect_is(outputs_longlat, "list") + expect_length(outputs_longlat, 4) + expect_equal(names(outputs_longlat), c("corner_coord", "polygon", "outliers","UTM_code")) + expect_equal(outputs_longlat[1:3],outputs_proj_coord) +}) + +test_that("check_plot_coord, trust_GPS_corners", { + proj_coord0 <- data.frame( + X = rep(c(10,10,100,100), e=7) + 200000, + Y = rep(c(10,100,100,10), e=7) + 9000000 + ) + proj_coord <- proj_coord0 + proj_coord$X <- proj_coord0$X + rep(sample(-3:3),4) + proj_coord$Y <- proj_coord0$Y + rep(sample(-3:3),4) + rel_coord <- data.frame( + X = rep(c(0,0,95,95),e=7), + Y = rep(c(0,95,95,0),e=7) + ) + corner_ID <- rep(c("SW","NW","SE","NE"),e=7) + + expect_warning( + check_plot_coord(proj_coord = proj_coord[-(1:3),], rel_coord = rel_coord[-(1:3),], trust_GPS_corners = T, rm_outliers = T, corner_ID = corner_ID[-(1:3)], draw_plot = F),"At least one corner has less than 5 measurements. We suggest using the argument trust_GPS_corners = FALSE" + ) + res_trust_GPS_corners_T <- 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) + expect_equivalent(res_trust_GPS_corners_T$corner_coord[,1:2] , unique(proj_coord0)[c(1,4,3,2),]) + + res_trust_GPS_corners_F <- check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord, trust_GPS_corners = F, rm_outliers = T, corner_ID = corner_ID, draw_plot = F) + expect_equivalent(res_trust_GPS_corners_T$corner_coord[,3:5] , res_trust_GPS_corners_F$corner_coord[,3:5]) + expect_equivalent(res_trust_GPS_corners_T$corner_coord[,1:2], round(res_trust_GPS_corners_F$corner_coord[,1:2],digits = -1)) + + # Test when there's only 4 measures and trust_GPS_corners = F + expect_equal(res_trust_GPS_corners_F , check_plot_coord(proj_coord = proj_coord0[c(1,8,15,22),], rel_coord = rel_coord[c(1,8,15,22),], trust_GPS_corners = F, rm_outliers = T, draw_plot = F, corner_ID = corner_ID[c(1,8,15,22)])) +}) + +test_that("check_plot_coord, origin corner", { + proj_coord <- data.frame( + X = rep(c(0,0,100,100), e=7) + 200000 + rep(sample(-3:3),4), + Y = rep(c(0,100,100,0), e=7) + 9000000 + rep(sample(-3:3),4) + ) + corner_ID <- rep(c("SW","NW","SE","NE"),e=7) + + # Test when the origin is not the South-East corner + rel_coord_NE <- data.frame( + X = rep(c(100,100,0,0),e=7), + Y = rep(c(100,0,0,100),e=7) + ) + rel_coord_SW <- data.frame( + X = rep(c(0,0,100,100),e=7), + Y = rep(c(0,100,100,0),e=7) + ) + res_NE <- check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord_NE, trust_GPS_corners = T, rm_outliers = T, corner_ID = corner_ID, draw_plot = F) + res_SW <- check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord_SW, trust_GPS_corners = T, rm_outliers = T, corner_ID = corner_ID, draw_plot = F) + + expect_equivalent(res_SW$corner_coord[c("x_rel","y_rel")] , res_NE$corner_coord[c("x_rel","y_rel")]) + expect_equivalent(res_SW$corner_coord[c("x_proj","y_proj","corner_ID")] , res_NE$corner_coord[c(3,4,1,2),c("x_proj","y_proj","corner_ID")]) + + # Test when the origin of the relative coordinates is not(0;0) + rel_coord_SW <- rel_coord_SW + 200 + res_SW_200 <- check_plot_coord(proj_coord = proj_coord, rel_coord = rel_coord_SW, trust_GPS_corners = T, rm_outliers = T, corner_ID = corner_ID, draw_plot = F) + expect_equal(res_SW$corner_coord[c("x_proj","y_proj","corner_ID")] , res_SW_200$corner_coord[c("x_proj","y_proj","corner_ID")]) + expect_equal(res_SW$corner_coord[c("x_rel","y_rel")]+200 , res_SW_200$corner_coord[c("x_rel","y_rel")]) +}) + + +test_that("check_plot_coord, tree coordinates", { + proj_coord <- data.frame( + X = rep(c(0,0,100,100), e=7) + 200000 + rep(sample(-3:3),4), + Y = rep(c(0,100,100,0), e=7) + 9000000 + rep(sample(-3:3),4) + ) + rel_coord <- data.frame( + X = rep(c(0,0,100,100),e=7), + Y = rep(c(0,100,100,0),e=7) + ) + corner_ID <- rep(c("SW","NW","SE","NE"),e=7) + tree_df <- data.frame( + X = c(10,20,30), + Y = c(10,20,30) + ) + 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")) + expect_is(res$treeproj_coord, "data.frame") + expect_equal(dim(res$treeproj_coord), c(nrow(tree_df),2)) + expect_equal(res$treeproj_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), + Y = rep(c(0,100,100,0), e=4) + ) + proj_coord[1:4,] <- matrix(c(25,0,-25,0,0,25,0,-25), ncol=2, byrow = F) + proj_coord$X = proj_coord$X + 200000 + proj_coord$Y = proj_coord$Y + 9000000 + + rel_coord <- data.frame( + X = rep(c(0,0,100,100),e=4), + Y = rep(c(0,100,100,0),e=4) + ) + 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 39fdf64737e60e696b393e2b19e96af74b9937f4 Mon Sep 17 00:00:00 2001 From: ArthurBailly Date: Mon, 18 Nov 2024 11:29:04 +0100 Subject: [PATCH 2/2] Added documentations --- NAMESPACE | 1 + man/checkPlotCoord.Rd | 6 +-- man/check_plot_coord.Rd | 103 ++++++++++++++++++++++++++++++++++++++++ man/divide_plot.Rd | 2 +- 4 files changed, 108 insertions(+), 4 deletions(-) create mode 100644 man/check_plot_coord.Rd diff --git a/NAMESPACE b/NAMESPACE index c585726..f153429 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,6 +7,7 @@ export(bilinear_interpolation) export(cacheManager) export(cachePath) export(checkPlotCoord) +export(check_plot_coord) export(computeAGB) export(computeE) export(computeFeldRegion) diff --git a/man/checkPlotCoord.Rd b/man/checkPlotCoord.Rd index 0b5a520..e07215a 100644 --- a/man/checkPlotCoord.Rd +++ b/man/checkPlotCoord.Rd @@ -33,17 +33,17 @@ checkPlotCoord( \item{drawPlot}{a logical indicating if the plot design should be displayed and returned} -\item{treeCoord}{(optional) a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively} +\item{treeCoord}{a data frame containing at least the relative tree coordinates (field/local coordinates), with X and Y corresponding to the first and second columns respectively} } \value{ Returns a list including : \itemize{ -\item \code{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 +\item \code{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 \item \code{polygon}: a spatial polygon \item \code{outliers}: a data frame containing the projected coordinates, the ID (if cornerID is supplied) and the row number of GPS points considered outliers \item \code{plotDesign}: if \code{drawPlot} is TRUE, a ggplot object corresponding to the design of the plot \item \code{codeUTM}: if \code{longlat} is supplied, a character containing the UTM code of the GPS coordinates -\item \code{treeProjCoord}: if \code{treeCoord} is supplied, a data frame containing the coordinates of the trees in the projected coordinate system +\item \code{treeProjCoord}: if \code{treeCoord} is supplied, a data frame containing the coordinates of the trees in the projected coordinate system (x_proj and y_proj) } } \description{ diff --git a/man/check_plot_coord.Rd b/man/check_plot_coord.Rd new file mode 100644 index 0000000..99e0f04 --- /dev/null +++ b/man/check_plot_coord.Rd @@ -0,0 +1,103 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/check_plot_coord.R +\name{check_plot_coord} +\alias{check_plot_coord} +\title{Check coordinates of plot corners and trees} +\usage{ +check_plot_coord( + 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 +) +} +\arguments{ +\item{proj_coord}{(optional, if longlat is not supplied) a data frame containing the plot corner coordinates in a projected coordinate system, with X and Y corresponding to the first and second columns respectively.} + +\item{longlat}{(optional, if proj_coord is not supplied) a data frame containing the plot corner coordinates in longitude and latitude, with longitude and latitude corresponding to the first and second columns respectively} + +\item{rel_coord}{a data frame containing the plot corner coordinates in the relative coordinates system (that of the field), with X and Y corresponding to the first and second columns respectively, and with the same row order than longlat or proj_coord} + +\item{trust_GPS_corners}{a logical indicating whether or not you trust the GPS coordinates of the plot's corners. See details.} + +\item{corner_ID}{a vector indicating the ID of the corners (e.g c("SE","SE",...)) in the case you have multiple measurements for each corner} + +\item{max_dist}{a numeric giving the maximum distance (in meters) above which GPS measurements should be considered outliers (default 15 m)} + +\item{rm_outliers}{a logical indicating if detected outliers are removed from the coordinate calculation} + +\item{draw_plot}{a logical indicating if the plot design should be displayed and returned} + +\item{tree_df}{a data frame containing the relative coordinates (field/local coordinates) of the trees and optional other tree metrics} + +\item{tree_coords}{a character vector of length 2 containing the column names of the relative coordinates in tree_df} +} +\value{ +Returns a list including : +\itemize{ +\item \code{corner_coord}: a data frame containing the projected coordinates (x_proj and y_proj), the relative coordinates (x_rel and y_rel), and the ID (if corner_ID is supplied) of the 4 corners of the plot +\item \code{polygon}: a spatial polygon +\item \code{outliers}: a data frame containing the projected coordinates, the ID (if corner_ID is supplied) and the row number of GPS points considered outliers +\item \code{plot_design}: if \code{draw_plot} is TRUE, a ggplot object corresponding to the design of the plot +\item \code{UTM_code}: if \code{longlat} is supplied, a character containing the UTM code of the GPS coordinates +\item \code{tree_proj_coord}: if \code{tree_df} is supplied, a data frame containing the coordinates of the trees in the projected coordinate system (x_proj and y_proj) +} +} +\description{ +Quality check of plot corner and tree coordinates. +} +\details{ +If trust_GPS_corners is TRUE, corner coordinates in the projected coordinate system are averaging by corner (if multiple measures) and outlier corners are identified sequentially using these averages and the max_dist argument. Then, projected coordinates of the trees are calculated from the local coordinates using a bilinear interpolation that follows the correspondence of the corners between these two coordinate systems. Be aware that this projection only works if the plot, in the relative coordinates system, is rectangular (ie, has 4 right angles). + +If trust_GPS_corners is FALSE, corner coordinates in the projected coordinate system are calculated by a procrust analysis that preserves the shape and dimensions of the plot in the local coordinate system. Outlier corners are also identified sequentially and projected coordinates of the trees are calculated by applying the resulting procrust analysis. + +If longlat is supplied instead of proj_coord, the function will first convert the long/lat coordinates into UTM coordinates. An error may result if the parcel is located right between two UTM zones. In this case, the user has to convert himself his long/lat coordinates into any projected coordinates which have the same dimension than his local coordinates (in meters most of the time). +} +\examples{ +proj_coord <- data.frame( + X = c( + runif(5, min = 9, max = 11), runif(5, min = 8, max = 12), + runif(5, min = 80, max = 120), runif(5, min = 90, max = 110) + ), + Y = c( + runif(5, min = 9, max = 11), runif(5, min = 80, max = 120), + runif(5, min = 8, max = 12), runif(5, min = 90, max = 110) + ) +) +proj_coord <- proj_coord + 1000 +rel_coord <- data.frame( + X = c(rep(0, 10), rep(100, 10)), + Y = c(rep(c(rep(0, 5), rep(100, 5)), 2)) +) +corner_ID <- rep(c("SW","NW","SE","NE"),e=5) + +aa <- check_plot_coord( + proj_coord = proj_coord, rel_coord = rel_coord, + trust_GPS_corners = TRUE, corner_ID = corner_ID, + rm_outliers = FALSE , draw_plot = FALSE +) +bb <- check_plot_coord( + proj_coord = proj_coord, rel_coord = rel_coord, + trust_GPS_corners = FALSE, + rm_outliers = TRUE, max_dist = 10, + draw_plot = FALSE +) +\donttest{ +check_plot_coord( + proj_coord = proj_coord, rel_coord = rel_coord, + trust_GPS_corners = TRUE, corner_ID = corner_ID, + rm_outliers = TRUE , max_dist = 10, + draw_plot = TRUE +) +} + +} +\author{ +Arthur PERE, Maxime REJOU-MECHAIN, Arthur BAILLY +} diff --git a/man/divide_plot.Rd b/man/divide_plot.Rd index 56c072a..dc7e9a0 100644 --- a/man/divide_plot.Rd +++ b/man/divide_plot.Rd @@ -31,7 +31,7 @@ If tree_coord isn't supplied, returns a data-frame containing as many rows as th \itemize{ \item \code{plot_ID_corner}: If dealing with multiple plots : the plot code \item \code{subplot_id}: The automatically generated subplot code, using the following rule : subplot_X_Y -\item \code{x} and \code{y} (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. +\item \code{x_rel} and \code{y_rel} (or the column names provided by rel_coord) : the relative X-axis and Y-axis coordinates of subplots corners. \item \code{x_proj} and \code{y_proj} (or the column names provided by proj_coord) : if proj_coord is supplied, the projected X-axis and Y-axis coordinates of subplots corners }