Spatialized trees and forest stands metrics with BIOMASS

Overview

BIOMASS enables users to manage their plots by :

calculating the projected/geographic coordinates of the plot’s corners and trees from the relative coordinates (or local coordinates, ie, those of the field)
visualizing the plots
validating plot’s corner and tree coordinates with LiDAR data
dividing plots into subplots
summarizing the desired information at subplot level

Required data

Two data frames are required for the rest :

a data frame of corner coordinates, containing at least :
- the name of the plots if there are several plots
- the coordinates of the plot’s corners in the geographic or projected coordinate system (the GPS coordinates or another projected coordinates)
- the coordinates of the plot’s corners in the relative coordinate system (the local or field coordinates)

In this vignette, we will use the Nouragues dataset as an exemple :

cornerCoord <- read.csv(system.file("external", "Coord.csv",package = "BIOMASS", mustWork = TRUE))

kable(head(cornerCoord), digits = 5, row.names = FALSE, caption = "Head of the table cornerCoord")

Head of the table cornerCoord
Plot	Corners	Lat	Long	xRel	xProj	yProj
NB1	NB1_SE	4.06692	52.68883	100	687482.0	449720.8
NB1	NB1_SE	4.06694	52.68883	100	687480.9	449722.6
NB1	NB1_SE	4.06694	52.68884	100	687482.8	449722.2
NB1	NB1_SE	4.06692	52.68882	100	687480.7	449720.7
NB1	NB1_SE	4.06695	52.68883	100	687481.5	449723.5
NB1	NB1_SE	4.06695	52.68884	100	687483.1	449723.5

Note that both Lat/Long and xProj/yProj coordinates are included in this dataset but only one of these coordinate systems is needed.

a data frame of trees coordinates, containing at least :
- the name of the plots if there are several plots
- the tree coordinates in the relative coordinate system (the local/field one)
- the desired information about trees, such as diameter, wood density, height, etc. (see BIOMASS vignette)

trees <- read.csv(system.file("external", "NouraguesPlot.csv",package = "BIOMASS", mustWork = TRUE))

kable(head(trees), digits = 3, row.names = FALSE, caption = "Head of the table trees")

Head of the table trees
plot	xRel	yRel	D	WD	H
NB1	1.30	4.7	11.459	0.643	12
NB1	2.65	4.3	11.618	0.580	16
NB1	4.20	6.9	83.875	0.591	40
NB1	5.90	4.7	14.961	0.568	18
NB1	6.40	4.1	36.765	0.530	27
NB1	13.50	2.3	13.528	0.409	20

Checking plot’s coordinates

The user is faced with two situations :

The GPS coordinates of the plot corners are considered very accurate or enough measurements have been made to be confident in the accuracy of their mean. In this case, the shape of the plot measured on the field will follow the GPS coordinates of the plot corners when projected into the projected/geographic coordinate system. See 3.1.1
Too few measurements of the GPS coordinates of plot corners have been taken and/or are not reliable. In this case, the shape of the plot measured on the field is considered to be accurate and the GPS corner coordinates will be recalculated to fit the shape and dimensions of the plot (the projected coordinates fit the relative coordinates). See 3.1.2

In both case, the use of the checkPlotCoord() function is recommended as a first step.

Checking the corners of the plot

The checkPlotCoord function handles both situations using the trustGPScorners argument (= TRUE or FALSE).

You can give either GPS coordinates (with the longlat argument) or another projected coordinates (with the projCoord argument) for the corner coordinates.

If we rely on the GPS coordinates of the corners

If enough coordinates have been recorded for each corner (for more information, see the CEOS good practices protocol, section A.1.3.1 ), you will need to provide the corner names via the cornerID argument. In this case, each coordinates will be averaging by corner, resulting in 4 reference coordinates. The function can also detect and remove GPS outliers using the rmOutliers and maxDist arguments.

If only 4 GPS measurements have been taken with a high degree of accuracy (by a geometer, for example), or if you have averaged your measurements by yourself, you can supply these 4 GPS coordinates to the function.

check_plot_trust_GPS <- checkPlotCoord(
  longlat = cornerCoord[c("Long", "Lat")],
  # OR, if exists : projCoord = cornerCoord[c("xProj", "yProj")], 
  relCoord = cornerCoord[c("xRel", "yRel")],
  trustGPScorners = T, cornerID = cornerCoord$Corners,
  drawPlot = TRUE,
  maxDist = 10, rmOutliers = TRUE )
#> Loading required namespace: proj4

The two blue arrows represent the relative coordinate system when projected into the projected coordinate system.

If we rely on the shape of the plot measured on the field

Let’s degrade the data to simulate the fact that we only have 8 GPS coordinates that we don’t trust.

degradedCornerCoord <- cornerCoord[c(1:2,11:12,21:22,31:32),]
  
check_plot_trust_field <- checkPlotCoord(
  longlat = degradedCornerCoord[, c("Long", "Lat")],
  # OR projCoord = cornerCoord[c("xProj", "yProj")], 
  relCoord = degradedCornerCoord[, c("xRel", "yRel")],
  trustGPScorners = F,
  drawPlot = TRUE )

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 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.

check_plot_trust_GPS <- checkPlotCoord(
  longlat = cornerCoord[c("Long", "Lat")],
  # OR projCoord = cornerCoord[c("xProj", "yProj")], 
  relCoord = cornerCoord[c("xRel", "yRel")],
  trustGPScorners = T, cornerID = cornerCoord$Corners,
  drawPlot = TRUE,
  maxDist = 10, rmOutliers = TRUE,
  treeCoord = trees[c("xRel","yRel")] )


check_plot_trust_field <- checkPlotCoord(
  longlat = degradedCornerCoord[, c("Long", "Lat")],
  # OR projCoord = cornerCoord[c("xProj", "yProj")], 
  relCoord = degradedCornerCoord[, c("xRel", "yRel")],
  trustGPScorners = F,
  drawPlot = TRUE,
  treeCoord = trees[c("xRel","yRel")] )

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 :

plot_to_change <- check_plot_trust_GPS$plotDesign
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:

treeGPSCoord <- as.data.frame( proj4::project(check_plot_trust_GPS$treeProjCoord, proj = check_plot_trust_GPS$codeUTM, inverse = TRUE) )

Integrating LiDAR data

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) :

subplots <- divide_plot(
  rel_coord = check_plot_trust_field$cornerCoord[c("Xrel","Yrel")],
  proj_coord = check_plot_trust_field$cornerCoord[c("X","Y")], #optional
  grid_size = 50 # or c(50,50)
  )

For the purposes of summarizing and representing subplots (coming in the next section), the function returns the coordinates of subplot corners (keeping column names) but can also assign to each tree its subplot with the tree_coord :

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)
  tree_coord = trees
  )

Note that tree_coord can contain any number of columns, but the first two must be the coordinates of the trees in the relative system.

The function can handle as many plots as you wish, using the plot_ID_corner and plot_ID_tree arguments :

multiple_plots <- rbind(check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")], check_plot_trust_GPS$cornerCoord[c("Xrel","Yrel")])
multiple_plots$plot_ID = rep(c("NB1","NB2"),e=4)
multiple_trees <- rbind(trees, trees %>% dplyr::mutate(plot="NB2"))

mulitple_subplots <- divide_plot(
  rel_coord = multiple_plots[c("Xrel","Yrel")],
  grid_size = 50, 
  plot_ID_corner = multiple_plots$plot_ID,
  tree_coord = multiple_trees,
  plot_ID_tree = multiple_trees$plot
)

Arthur Bailly

2024-11-07