From 8bcb7b8109195c3a7fc91b2cd3cb3bb05f94d1d9 Mon Sep 17 00:00:00 2001 From: Vincyane Badouard Date: Mon, 16 Aug 2021 16:20:38 -0300 Subject: [PATCH] createcanopy and treefromthesky functions #38 #39 --- NAMESPACE | 2 + R/Crowns_draft.Rmd | 72 ++++++++++++++++++ R/createcanopy.R | 105 +++++++++++++++++++++++++++ R/treefelling.R | 57 +++++++++------ man/createcanopy.Rd | 36 +++++++++ man/st_tree.Rd | 8 +- man/treefromthesky.Rd | 31 ++++++++ tests/testthat/test-createcanopy.R | 26 +++++++ tests/testthat/test-st_tree.R | 20 +++-- tests/testthat/test-treefromthesky.R | 26 +++++++ 10 files changed, 352 insertions(+), 31 deletions(-) create mode 100644 R/Crowns_draft.Rmd create mode 100644 R/createcanopy.R create mode 100644 man/createcanopy.Rd create mode 100644 man/treefromthesky.Rd create mode 100644 tests/testthat/test-createcanopy.R create mode 100644 tests/testthat/test-treefromthesky.R diff --git a/NAMESPACE b/NAMESPACE index 708f90f2..6b5dece4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,7 @@ export("%>%") export(ONFGuyafortaxojoin) export(addtreedim) +export(createcanopy) export(directionalfellingsuccessdef) export(exploitablefuelwoodvolume) export(futurereserve) @@ -16,6 +17,7 @@ export(selected) export(st_tree) export(timberharvestedvolume) export(treefelling) +export(treefromthesky) export(treeselection) importFrom(BIOMASS,computeAGB) importFrom(BIOMASS,getWoodDensity) diff --git a/R/Crowns_draft.Rmd b/R/Crowns_draft.Rmd new file mode 100644 index 00000000..c1171245 --- /dev/null +++ b/R/Crowns_draft.Rmd @@ -0,0 +1,72 @@ +--- +title: "Crowns_draft" +author: "Vincyane Badouard" +date: "16/08/2021" +output: html_document +--- + +```{r} +inventory <- addtreedim(inventorycheckformat(Paracou6_2016)) +dat <- inventory[679,] +``` + + +```{r} +# The crown +treefromthesky <- function(dat){ + + Crown <- dat %>% + mutate(xCrown = Xutm, # X centroid + yCrown = Yutm, # Y ventroid + exCrown = CrownDiameter/2, + eyCrown = CrownHeight/2) %>% + st_as_sf(coords = c("xCrown", "yCrown")) # ellipse centroid coordinates + Crown <- st_ellipse(Crown, Crown$exCrown, Crown$eyCrown) # create the ellipse + +} +``` + +```{r} +library(ggplot2) +g <- ggplot() + + geom_sf(data = st_as_sf(inventory, coords = c("Xutm", "Yutm")), aes(label = ScientificName)) + + geom_sf(data = Crown, fill = "forestgreen") # trees polygons + +plotly::ggplotly(g) +``` + +```{r} +Canopy <- inventory %>% + group_by(idTree) %>% # for each tree + do(Crowns = # inform geometry. # Filling a column from a function whose input is a table + treefromthesky(.) %>% + st_as_text()) %>% # as text to easy join with a non spacial table + tidyr::unnest(Crowns) # here to pass from list to character + +inventory <- left_join(inventory, Canopy, by = "idTree") # join spatial filtered inventory and non spatial complete inventory +``` + +```{r} +#The small ones first so that they are behind the big ones on the plot +inventory <- arrange(inventory, TreeHeight) + +ggplot() + + geom_sf(data = getgeometry(inventory, Crowns), + aes(alpha = TreeHeight), # label = paste(idTree, Species), + fill = "forestgreen") +``` + +# Plotly version +```{r} +#The small ones first so that they are behind the big ones on the plot +inventory <- arrange(inventory, TreeHeight) %>% + filter(TreeHeight >30) + +g <- ggplot() + + geom_sf(data = getgeometry(inventory, Crowns), + aes(label = paste(idTree, Species), alpha = 1/2), + fill = "forestgreen") # trees polygons + +plotly::ggplotly(g) +``` + diff --git a/R/createcanopy.R b/R/createcanopy.R new file mode 100644 index 00000000..5758a3a6 --- /dev/null +++ b/R/createcanopy.R @@ -0,0 +1,105 @@ +#' createcanopy +#' +#' @param inventory (data.frame) +#' +#' @return a dataframe with a column 'Crowns' containing the ellipses +#' (sfc_POLYGON) as trees crown, with their diameter and height filled in, +#' representing trees from the sky. +#' @export +#' +#' @importFrom dplyr group_by do left_join +#' @importFrom sf st_as_text +#' @importFrom tidyr unnest +#' +#' @examples +#' data(Paracou6_2016) +#' Paracou6_2016 <- dplyr::slice(Paracou6_2016, 1:10) +#' +#' inventory <- addtreedim(inventorycheckformat(Paracou6_2016)) +#' +#' inventory <- createcanopy(inventory) +#' +#' # The small ones first so that they are behind the big ones on the plot +#' inventory <- dplyr::arrange(inventory, TreeHeight) +#' library(ggplot2) +#' ggplot() + +#' geom_sf(data = getgeometry(inventory, Crowns), +#' aes(alpha = TreeHeight), +#' fill = "forestgreen") +#' +createcanopy <- function(inventory){ + + # Arguments check + + if(!inherits(inventory, "data.frame")) + stop("The 'inventory' argument of the 'createcanopy' function must be a data.frame") + + # Global variables + idTree <- Crowns <- . <- NULL + + # Function + + Canopy <- inventory %>% + + group_by(idTree) %>% # for each tree + + do(Crowns = # inform geometry. # do: filling a column from a function whose input is a table + treefromthesky(.) %>% + + st_as_text()) %>% # as text to easy join with a non spacial table + tidyr::unnest(Crowns) # here to pass from list to character + + inventory <- left_join(inventory, Canopy, by = "idTree") # join the column 'Crowns' to the inventory + +} + +#' treefromthesky +#' +#' @param dat 1 row data.frame with columns: +#' Xutm, Yutm (Tree coordinates), CrownDiameter, CrownHeight +#' +#' @return an ellipse (sfc_POLYGON) as a crown, with its diameter and height +#' filled in, representing the tree from the sky. +#' @export +#' +#' @importFrom dplyr mutate +#' @importFrom sf st_as_sf +#' @importFrom nngeo st_ellipse +#' +#' @examples +#' inventory <- addtreedim(inventorycheckformat(Paracou6_2016)) +#' dat <- inventory[679,] +#' +#' Crown <- treefromthesky(dat) +#' +#' library(ggplot2) +#' ggplot() + +#' geom_sf(data = st_as_sf(inventory, coords = c("Xutm", "Yutm"))) + +#' geom_sf(data = Crown, fill = "forestgreen") # trees polygons +#' +treefromthesky <- function(dat){ + + # Arguments check + + if(!inherits(dat, "data.frame")) + stop("The 'dat' argument of the 'treefromthesky' function must be a data.frame") + + if(nrow(dat)!=1) + stop("the data.frame given in the 'dat' argument + of the 'treefromthesky' function must contain only 1 row") + + # Global variables + Xutm <- Yutm <- CrownDiameter <- CrownHeight <- NULL + + # Function + + Crown <- dat %>% + mutate(xCrown = Xutm, # X centroid + yCrown = Yutm, # Y ventroid + exCrown = CrownDiameter/2, # crown radius + eyCrown = CrownHeight/2) %>% # half crown height + st_as_sf(coords = c("xCrown", "yCrown")) # ellipse centroid coordinates + + Crown <- st_ellipse(Crown, Crown$exCrown, Crown$eyCrown) # create the ellipse + +} diff --git a/R/treefelling.R b/R/treefelling.R index a5b97c72..c2058123 100644 --- a/R/treefelling.R +++ b/R/treefelling.R @@ -17,10 +17,7 @@ #' "DamageTreesPoints", "DeadTreesPoints" vectors #' @export #' -#' @importFrom dplyr filter -#' @importFrom dplyr group_by -#' @importFrom dplyr do -#' @importFrom dplyr left_join +#' @importFrom dplyr group_by do left_join #' @importFrom sf st_as_text #' @importFrom tidyr unnest #' @@ -401,12 +398,16 @@ rotatepolygon <- function( #' grapple (see the \code{GrappleLength} argument of the #' \code{\link{loggingparameters}} function). In other cases, the fall will be #' made from the base of the tree towards the trail. The orientation of the -#' fall succeeds or fails according to a Bernoulli law where the probability of success -#' is by default 60%, and can be changed with the \code{advancedloggingparameters} argument. +#' fall succeeds or fails according to a Bernoulli law where the probability +#' of success is by default 60%, and can be changed with the +#' \code{advancedloggingparameters} argument. +#' +#' @param dat 1 row data.frame with columns: Xutm, Yutm, CrownDiameter, +#' CrownHeight, DBH, TrunkHeight, TreeHeight, TreeFellingOrientationSuccess #' -#' @param dat (dataframe) #' @param MainTrail (sfg) #' @param ScndTrail (sfg) +#' #' @param fuel Fuel wood exploitation: no exploitation = "0", damage #' exploitation in fuelwood = "1", exploitation of hollow trees and damage in #' fuelwood = "2" @@ -499,6 +500,10 @@ st_tree <- function( if(!inherits(dat, "data.frame")) stop("The 'dat' argument of the 'st_tree' function must be data.frame") + if(nrow(dat)!=1) + stop("the data.frame given in the 'dat' argument + of the 'st_tree' function must contain only 1 row") + if (!any(fuel == "0" || fuel == "1"|| fuel == "2"|| is.null(fuel))) stop("The 'fuel' argument of the 'st_tree' function must be '0', '1', '2' or NULL") @@ -574,10 +579,11 @@ st_tree <- function( )) } - # To direct only to avoid damage to future and reserve trees. Winching: Foot before. + + # # To direct only to avoid damage to future and reserve trees. Winching: Foot before. # if (directionalfelling == "1"&& (fuel !="1" || fuel !="2")) { # if(dat$TreeFellingOrientationSuccess == "1"){ - + # # }else{ # else random felling # RandomAngle <- as.numeric(sample(c(0:359.9), size = 1)) # A <- st_difference(st_union( @@ -588,16 +594,25 @@ st_tree <- function( # # } + + # Scenarios with track orientation: + # Compute the last angle of the right-angled triangle (see vignette figure) + # TreefallOrientation is between the mimimun (30°) and the maximum (45°) angle + TreefallOrientation <- as.numeric(sample(c(advancedloggingparameters$MinTreefallOrientation: + advancedloggingparameters$MaxTreefallOrientation), size = 1)) + OppAng <- 180-(90 + TreefallOrientation) + + # To direct !!to avoid damage to future and reserve trees!! + track orientation. Winching: Foot before. if(directionalfelling == "2"&& (fuel !="1" || fuel !="2")){ if(dat$TreeFellingOrientationSuccess == "1"){ A <- st_difference(st_union( - rotatepolygon(Trunk, angle = 240 + theta, fixed = Foot), # turned trunk - rotatepolygon(Crown, angle = 240 + theta, fixed = Foot) # turned crown + rotatepolygon(Trunk, angle = 180 + OppAng + theta, fixed = Foot), # turned trunk + rotatepolygon(Crown, angle = 180 + OppAng + theta, fixed = Foot) # turned crown )) B <- st_difference(st_union( - rotatepolygon(Trunk, angle = 120 + theta, fixed = Foot), # turned trunk - rotatepolygon(Crown, angle = 120 + theta, fixed = Foot) # turned crown + rotatepolygon(Trunk, angle = 180 - OppAng + theta, fixed = Foot), # turned trunk + rotatepolygon(Crown, angle = 180 - OppAng + theta, fixed = Foot) # turned crown )) }else{ # else random felling RandomAngle <- as.numeric(sample(c(0:359.9), size = 1)) @@ -621,21 +636,21 @@ st_tree <- function( if(TrailDist <= advancedloggingparameters$GrappleLength){ # <= 6m (= grapple length) -> winching by grapple -> crown to trail A <- st_difference(st_union( - rotatepolygon(Trunk, angle = theta + 60 , fixed = Foot), # turned trunk - rotatepolygon(Crown, angle = theta + 60 , fixed = Foot) # turned crown + rotatepolygon(Trunk, angle = theta + OppAng , fixed = Foot), # turned trunk + rotatepolygon(Crown, angle = theta + OppAng , fixed = Foot) # turned crown )) B <- st_difference(st_union( - rotatepolygon(Trunk, angle = 300 + theta, fixed = Foot), # turned trunk - rotatepolygon(Crown, angle = 300 + theta, fixed = Foot) # turned crown + rotatepolygon(Trunk, angle = 360 - OppAng + theta, fixed = Foot), # turned trunk + rotatepolygon(Crown, angle = 360 - OppAng + theta, fixed = Foot) # turned crown )) } else { # > 6m -> winching by cable -> foot to trail A <- st_difference(st_union( - rotatepolygon(Trunk, angle = 240 + theta, fixed = Foot), # turned trunk - rotatepolygon(Crown, angle = 240 + theta, fixed = Foot) # turned crown + rotatepolygon(Trunk, angle = 180 + OppAng + theta, fixed = Foot), # turned trunk + rotatepolygon(Crown, angle = 180 + OppAng + theta, fixed = Foot) # turned crown )) B <- st_difference(st_union( - rotatepolygon(Trunk, angle = 120 + theta, fixed = Foot), # turned trunk - rotatepolygon(Crown, angle = 120 + theta, fixed = Foot) # turned crown + rotatepolygon(Trunk, angle = 180 - OppAng + theta, fixed = Foot), # turned trunk + rotatepolygon(Crown, angle = 180 - OppAng + theta, fixed = Foot) # turned crown )) } }else{ # else random felling diff --git a/man/createcanopy.Rd b/man/createcanopy.Rd new file mode 100644 index 00000000..d87b9af9 --- /dev/null +++ b/man/createcanopy.Rd @@ -0,0 +1,36 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/createcanopy.R +\name{createcanopy} +\alias{createcanopy} +\title{createcanopy} +\usage{ +createcanopy(inventory) +} +\arguments{ +\item{inventory}{(data.frame)} +} +\value{ +a dataframe with a column 'Crowns' containing the ellipses +(sfc_POLYGON) as trees crown, with their diameter and height filled in, +representing trees from the sky. +} +\description{ +createcanopy +} +\examples{ +data(Paracou6_2016) +Paracou6_2016 <- dplyr::slice(Paracou6_2016, 1:10) + +inventory <- addtreedim(inventorycheckformat(Paracou6_2016)) + +inventory <- createcanopy(inventory) + +# The small ones first so that they are behind the big ones on the plot +inventory <- dplyr::arrange(inventory, TreeHeight) +library(ggplot2) +ggplot() + + geom_sf(data = getgeometry(inventory, Crowns), + aes(alpha = TreeHeight), + fill = "forestgreen") + +} diff --git a/man/st_tree.Rd b/man/st_tree.Rd index 0960ba47..1112f090 100644 --- a/man/st_tree.Rd +++ b/man/st_tree.Rd @@ -14,7 +14,8 @@ st_tree( ) } \arguments{ -\item{dat}{(dataframe)} +\item{dat}{1 row data.frame with columns: Xutm, Yutm, CrownDiameter, +CrownHeight, DBH, TrunkHeight, TreeHeight, TreeFellingOrientationSuccess} \item{fuel}{Fuel wood exploitation: no exploitation = "0", damage exploitation in fuelwood = "1", exploitation of hollow trees and damage in @@ -41,8 +42,9 @@ crowns will be directed towards the trail if they can be accessed with a grapple (see the \code{GrappleLength} argument of the \code{\link{loggingparameters}} function). In other cases, the fall will be made from the base of the tree towards the trail. The orientation of the -fall succeeds or fails according to a Bernoulli law where the probability of success -is by default 60\%, and can be changed with the \code{advancedloggingparameters} argument. +fall succeeds or fails according to a Bernoulli law where the probability +of success is by default 60\%, and can be changed with the +\code{advancedloggingparameters} argument. } \examples{ MainTrail <- sf::st_linestring(matrix(c(286400, 583130, diff --git a/man/treefromthesky.Rd b/man/treefromthesky.Rd new file mode 100644 index 00000000..8935f723 --- /dev/null +++ b/man/treefromthesky.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/createcanopy.R +\name{treefromthesky} +\alias{treefromthesky} +\title{treefromthesky} +\usage{ +treefromthesky(dat) +} +\arguments{ +\item{dat}{1 row data.frame with columns: +Xutm, Yutm (Tree coordinates), CrownDiameter, CrownHeight} +} +\value{ +an ellipse (sfc_POLYGON) as a crown, with its diameter and height +filled in, representing the tree from the sky. +} +\description{ +treefromthesky +} +\examples{ +inventory <- addtreedim(inventorycheckformat(Paracou6_2016)) +dat <- inventory[679,] + +Crown <- treefromthesky(dat) + +library(ggplot2) +ggplot() + + geom_sf(data = st_as_sf(inventory, coords = c("Xutm", "Yutm"))) + + geom_sf(data = Crown, fill = "forestgreen") # trees polygons + +} diff --git a/tests/testthat/test-createcanopy.R b/tests/testthat/test-createcanopy.R new file mode 100644 index 00000000..3a6529ca --- /dev/null +++ b/tests/testthat/test-createcanopy.R @@ -0,0 +1,26 @@ +test_that("createcanopy", { + + # Test data + data(Paracou6_2016) + Paracou6_2016 <- dplyr::slice(Paracou6_2016, 1:10) + inventory <- addtreedim(inventorycheckformat(Paracou6_2016)) + + MatrixInventory <- as.matrix(Paracou6_2016) + + Rslt <- createcanopy(inventory) + + # Check the function arguments + expect_error(createcanopy(MatrixInventory), + regexp = "The 'inventory' argument of the 'createcanopy' function must be a data.frame") + + # Check output + expect_s3_class(Rslt, "data.frame") + expect_type(Rslt$Crowns, "character") + + }) + + +# ggplot() + +# geom_sf(data = getgeometry(Rslt, Crowns), +# aes(alpha = TreeHeight), +# fill = "forestgreen") diff --git a/tests/testthat/test-st_tree.R b/tests/testthat/test-st_tree.R index f40c311d..d7b102f3 100644 --- a/tests/testthat/test-st_tree.R +++ b/tests/testthat/test-st_tree.R @@ -48,24 +48,28 @@ test_that("st_tree", { advancedloggingparameters = loggingparameters()), regexp = "The 'dat' argument of the 'st_tree' function must be data.frame") - expect_error(st_tree(Paracou6_2016, directionalfelling = "2", + expect_error(st_tree(Paracou6_2016), + regexp = "the data.frame given in the 'dat' argument + of the 'st_tree' function must contain only 1 row") + + expect_error(st_tree(dat, directionalfelling = "2", advancedloggingparameters = loggingparameters(), fuel = TRUE), regexp = "The 'fuel' argument of the 'st_tree' function must be '0', '1', '2' or NULL") - expect_error(st_tree(Paracou6_2016, fuel = "2", + expect_error(st_tree(dat, fuel = "2", advancedloggingparameters = loggingparameters(), directionalfelling = TRUE), regexp = "The 'directionalfelling' argument of the 'st_tree' function must be '0', '1', '2' or NULL") - expect_error(st_tree(Paracou6_2016, + expect_error(st_tree(dat, fuel = "2", directionalfelling = "2", advancedloggingparameters = loggingparameters(), MainTrail = st_as_text(MainTrail), ScndTrail = st_as_text(ScndTrail)), regexp = "The 'MainTrail' and 'ScndTrail' arguments of the 'st_tree' function must be sfg") - expect_error(st_tree(Paracou6_2016, fuel = "2", + expect_error(st_tree(dat, fuel = "2", directionalfelling = "2", MainTrail = MainTrail, ScndTrail = ScndTrail, advancedloggingparameters = as.matrix(loggingparameters())), regexp = "The 'advancedloggingparameters' argument of the 'st_tree' function must be a list") @@ -88,10 +92,12 @@ test_that("st_tree", { Arrival <- st_point(as.numeric(unlist( # sfc to sfg sf::st_centroid(Rslt$A)))) + Orientation <- as.numeric(matlib::angle(c(Rslt$TrailPt[1] - Rslt$TrailPt[1], (Rslt$TrailPt[2]+10) - Rslt$TrailPt[2]), + c(Arrival[1] - Rslt$Foot[1], Arrival[2] - Rslt$Foot[2]), + degree = TRUE)) + - expect_equal(as.numeric(matlib::angle(c(Rslt$TrailPt[1] - Rslt$TrailPt[1], (Rslt$TrailPt[2]+10) - Rslt$TrailPt[2]), - c(Arrival[1] - Rslt$Foot[1], Arrival[2] - Rslt$Foot[2]), - degree = TRUE)), 30) + expect_true(Orientation >= 30 & Orientation <= 45) }) diff --git a/tests/testthat/test-treefromthesky.R b/tests/testthat/test-treefromthesky.R new file mode 100644 index 00000000..0e99e317 --- /dev/null +++ b/tests/testthat/test-treefromthesky.R @@ -0,0 +1,26 @@ +test_that("treefromthesky", { + + # Test data + data(Paracou6_2016) + inventory <- addtreedim(inventorycheckformat(Paracou6_2016)) + + dat <- inventory[1,] + MatrixInventory <- as.matrix(dat) + + Rslt <- treefromthesky(dat) + + + + # Check the function arguments + expect_error(treefromthesky(MatrixInventory), + regexp = "The 'dat' argument of the 'treefromthesky' function must be a data.frame") + + expect_error(treefromthesky(inventory), + regexp = "the data.frame given in the 'dat' argument + of the 'treefromthesky' function must contain only 1 row") + + + # Results class: + expect_s3_class(Rslt, "sfc_POLYGON") + + })