diff --git a/R/area-based.1.data.preparation.Rmd b/R/area-based.1.data.preparation.Rmd index 3d0ceee5938b00a86fcf73147b00376401042ab1..73461884de0777812b1d7572bb51300afe020fbc 100644 --- a/R/area-based.1.data.preparation.Rmd +++ b/R/area-based.1.data.preparation.Rmd @@ -242,7 +242,7 @@ Terrain statistics can be extracted from the cloud of ground points with altitud ```{r terrainStats, include=TRUE, warning=FALSE, message=FALSE, fig.width = 4, fig.height = 6} # use mapply to apply points2terrainStats function to each point cloud # while providing the coordinates of each center -terrain.stats <- mapply(function(x, y) +metrics.terrain <- mapply(function(x, y) { lidaRtRee::points2terrainStats(x, plots[y, c("X", "Y")], p.radius) }, @@ -250,13 +250,13 @@ x = llas.ground, y = as.list(1:length(llas.ground)), SIMPLIFY = FALSE) # bind results in data.frame -terrain.stats <- do.call(rbind, terrain.stats) +metrics.terrain <- do.call(rbind, metrics.terrain) # # compute terrain stats without specifying centre -terrain.stats2 <- lapply(llas.ground, lidaRtRee::points2terrainStats) -terrain.stats2 <- do.call(rbind, terrain.stats2) +metrics.terrain2 <- lapply(llas.ground, lidaRtRee::points2terrainStats) +metrics.terrain2 <- do.call(rbind, metrics.terrain2) # display results for comparison -head(cbind(terrain.stats, terrain.stats2), n=3) +head(cbind(metrics.terrain, metrics.terrain2), n=3) ``` ## Check field inventory data @@ -529,7 +529,7 @@ The following lines save the data required for the [area-based model calibration ```{r save.plots, include=TRUE} # save(plots, file="../data/aba.model/output/plots.rda") # save(llas.height, file="../data/aba.model/output/llas.height.rda", compress="bzip2") -# save(terrain.stats, file="../data/aba.model/output/terrain.stats.rda") +# save(metrics.terrain, file="../data/aba.model/output/metrics.terrain.rda") ``` diff --git a/R/area-based.2.model.calibration.Rmd b/R/area-based.2.model.calibration.Rmd index 542069cca22b06d211586e16fa380533ca9f2dd7..71319ea336c9cf2f0eb74f1041a73174e61ad7ba 100644 --- a/R/area-based.2.model.calibration.Rmd +++ b/R/area-based.2.model.calibration.Rmd @@ -61,15 +61,15 @@ The first lines of the terrain statistics are displayed hereafter. ```{r loadTerrainStats, echo = FALSE} # terrain statistics previously computed with (non-normalized) ground points inside each plot extent -load("../data/aba.model/output/terrain.stats.rda") -head(terrain.stats[, 1:3], n=3) +load("../data/aba.model/output/metrics.terrain.rda") +head(metrics.terrain[, 1:3], n=3) ``` The following lines ensure that the plots are ordered in the same way in the three data objects. ```{r orderData, include=TRUE} llas.height <- llas.height[plots$plotId] -terrain.stats <- terrain.stats[plots$plotId,] +metrics.terrain <- metrics.terrain[plots$plotId,] ``` # ALS metrics computation @@ -87,8 +87,8 @@ Point cloud metrics are computed with the function `lidaRtRee::cloudMetrics`, wh # define function for later use aba.pointMetricsFUN <- ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2) # apply function on each point cloud in list -point.metrics <- lidaRtRee::cloudMetrics(llas.height, aba.pointMetricsFUN) -round(head(point.metrics[, 1:8], n = 3),2) +metrics.points <- lidaRtRee::cloudMetrics(llas.height, aba.pointMetricsFUN) +round(head(metrics.points[, 1:8], n = 3),2) ``` ## Tree metrics @@ -101,14 +101,14 @@ aba.resCHM <- 0.5 # specify plot radius to exclude trees located outside plots plot.radius <- 15 # compute tree metrics -tree.metrics <- lidaRtRee::cloudTreeMetrics(llas.height, plots[, c("X", "Y")], +metrics.tree <- lidaRtRee::cloudTreeMetrics(llas.height, plots[, c("X", "Y")], plot.radius, res = aba.resCHM, func=function(x) { lidaRtRee::stdTreeMetrics(x, area.ha=pi*plot.radius^2/10000) }) -round(head(tree.metrics[, 1:5], n = 3), 2) +round(head(metrics.tree[, 1:5], n = 3), 2) ``` ## Other metrics @@ -116,9 +116,9 @@ round(head(tree.metrics[, 1:5], n = 3), 2) In case terrain metrics have been computed from the cloud of ground points only, they can also be added as variables, and so do other environmental variables which might be relevant in modeling. ```{r bindMetrics, include=TRUE} -metrics <- cbind(point.metrics[plots$plotId, ], - tree.metrics[plots$plotId, ], - terrain.stats[plots$plotId, 1:3]) +metrics <- cbind(metrics.points[plots$plotId, ], + metrics.tree[plots$plotId, ], + metrics.terrain[plots$plotId, 1:3]) ``` # Model calibration @@ -147,7 +147,7 @@ In this example, only tree metrics are selected in the basal area prediction mod ```{r modelPlot, include=TRUE, fig.height = 4.5, fig.width = 8} # check correlation between errors and other variables -round(cor(cbind(model.ABA$values$residual, plots[subsample, c("G.m2.ha","N.ha","D.mean.cm")], terrain.stats[subsample, 1:3])), 2)[1,] +round(cor(cbind(model.ABA$values$residual, plots[subsample, c("G.m2.ha","N.ha","D.mean.cm")], metrics.terrain[subsample, 1:3])), 2)[1,] # significance of correlation value cor.test(model.ABA$values$residual, plots[subsample, variable]) # plot predicted VS field values @@ -158,20 +158,20 @@ abline(h = 0, lty = 2) ``` In case only point cloud metrics are used as potential inputs, the errors are hardly better distributed. Coloring points by ownership shows that plots located in private forests have the largest basal area values which tend to be under-estimated. -```{r point.metricsOnly, include=TRUE, fig.height = 4.5, fig.width = 8} -model.ABA.point.metrics <- lidaRtRee::ABAmodel(plots[subsample,variable], point.metrics[subsample,], transform="boxcox", nmax=4, xy = plots[subsample, c("X", "Y")]) +```{r metrics.pointsOnly, include=TRUE, fig.height = 4.5, fig.width = 8} +model.ABA.metrics.points <- lidaRtRee::ABAmodel(plots[subsample,variable], metrics.points[subsample,], transform="boxcox", nmax=4, xy = plots[subsample, c("X", "Y")]) # renames outputs -row.names(model.ABA.point.metrics$stats) <- names(model.ABA.point.metrics$model) <- variable -# model.ABA.point.metrics$model[[variable]] -model.ABA.point.metrics$stats -# cor.test(model.ABA.point.metrics$values$residual, plots[subsample, variable]) +row.names(model.ABA.metrics.points$stats) <- names(model.ABA.metrics.points$model) <- variable +# model.ABA.metrics.points$model[[variable]] +model.ABA.metrics.points$stats +# cor.test(model.ABA.metrics.points$values$residual, plots[subsample, variable]) par(mfrow=c(1,2)) # plot predicted VS field values -lidaRtRee::ABAmodelPlot(model.ABA.point.metrics, main = variable, +lidaRtRee::ABAmodelPlot(model.ABA.metrics.points, main = variable, col = ifelse(plots$stratum == "public", "green", "blue")) legend("topleft", c("public", "private"), col = c("green", "blue"), pch = 1) plot(plots[subsample, c("G.m2.ha")], - model.ABA.point.metrics$values$residual, + model.ABA.metrics.points$values$residual, ylab = "Prediction errors", xlab = "Field values", col = ifelse(plots$stratum == "public", "green", "blue")) abline(h = 0, lty = 2) @@ -279,7 +279,7 @@ for (i in levels(plots[,strat])) subsample <- which(plots[,strat]==i) if (length(subsample)>0) { - model.ABA.stratified.none[[i]] <- lidaRtRee::ABAmodel(plots[subsample, variable], point.metrics[subsample,], transform="none", xy = plots[subsample,c("X", "Y")]) + model.ABA.stratified.none[[i]] <- lidaRtRee::ABAmodel(plots[subsample, variable], metrics.points[subsample,], transform="none", xy = plots[subsample,c("X", "Y")]) } } # combine list of models into single object diff --git a/R/area-based.3.mapping.and.inference.Rmd b/R/area-based.3.mapping.and.inference.Rmd index a9f559471d734af7f6fc2b3336f906f6ad1fa587..727f40ca4baba611f686a0587faf00cd726d64b6 100644 --- a/R/area-based.3.mapping.and.inference.Rmd +++ b/R/area-based.3.mapping.and.inference.Rmd @@ -68,25 +68,35 @@ sf::st_crs(public) <- sf::st_crs(forest) <- 2154 ## ALS metrics mapping -First the ALS metrics used as independent variables in the ABA prediction models have to be mapped on the whole area. For this the SAME functions used for computing the metrics from the plot-extent cloud metrics should be used. They are used inside the function `lidR::grid_metrics` which performs the computation the area of interest divided in square pixels. Pixel surface should be of the same size as calibration plots. Using round values is advisable in case other remote sensing data are considered for use. +First the ALS metrics used as independent variables in the ABA prediction models have to be mapped on the whole area. For this the SAME functions used for computing the metrics from the plot-extent cloud metrics should be used. Metrics are computed in each cell of the area of interest divided in square pixels. Pixel surface should be of the same size as calibration plots. Using round values is advisable in case other remote sensing data are considered for use. ```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6} # resolution for metrics map resolution <- 25 ``` -### Example on one tile +The following paragraph show how to compute metrics for a 100 x 100 m^2^ subsample of the first ALS tile. -The following code computes point metrics for a 100 x 100 m^2^ subsample of the first ALS tile. - -```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6} +```{r loadTile, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 6} # load ALS tile -point.cloud.height <- lidR::clip_rectangle(cata.height, 900100, 6447100, 900400, 6447400) -point.cloud.ground <- lidR::filter_ground(lidR::clip_rectangle(cata.altitude, 900100, 6447100, 900400, 6447400)) +point.cloud.height <- lidR::clip_rectangle(cata.height, 899600, 6447100, 899900, 6447400) +point.cloud.ground <- lidR::filter_ground(lidR::clip_rectangle(cata.altitude, 899600, 6447100, 899900, 6447400)) +``` + +### Point cloud metrics + +For computation of point cloud metrics, the same function used in `lidaRtRee::cloudMetrics` is now supplied to `lidR::grid_metrics`. Some metrics are displayed hereafter. + +```{r pointMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} # compute point metrics -point.metrics <- lidR::grid_metrics(point.cloud.height, ~lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2), res=resolution) -raster::plot(point.metrics[c("")]) +metrics.points <- lidR::grid_metrics(point.cloud.height, aba.pointMetricsFUN, res=resolution) +par(mfrow=c(1,3)) +for (i in c("zmean", "zsd", "imean")) +{ + raster::plot(metrics.points[[i]], main = i) +} ``` +### Tree metrics For tree metrics, three steps are required: @@ -94,7 +104,9 @@ For tree metrics, three steps are required: * summary statistics of tree features by target pixels (values calculated based on trees which apices are inside the pixel); * additional statistics about the percentage of cover and the mean canopy height (values calculated based on segmented tree crowns which are inside the target pixel). -```{r treeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6} +It is very important that tree segmentation parameters are the same as in the calibration steps, and that the function for aggregating tree attributes into raster metrics is the also the same as in the calibration step (except for the area which might differ slightly between the field plots and target cells). + +```{r treeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 6, fig.height = 6} # compute chm chm <- lidaRtRee::points2DSM(point.cloud.height, res = aba.resCHM) # replace NA, low and high values @@ -129,77 +141,239 @@ dummy <- lidaRtRee::rasterMetrics(r.treechm, names(dummy) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot") metrics.trees <- raster::addLayer(metrics.trees, dummy) # -metrics.trees +par(mfrow=c(1,1)) +raster::plot(chm, main = "CHM and detected trees") +sp::plot(trees, cex = trees$h/60, add = TRUE) +# +``` + +```{r plottreeMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +par(mfrow=c(1,3)) +for (i in c("Tree.meanH", "Tree.sdH", "Tree.density")) +{ + raster::plot(metrics.trees[[i]], main = i) +} ``` -Terrain statistics can be computed from the altitude values of ground points. +### Terrain metrics + +Terrain metrics can be computed from the altitude values of ground points. A wrapping function is required to pass `lidaRtRee::points2terrainStats` to `lidR::grid_metrics`. -```{r terrainMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 4, fig.height = 6} +```{r terrainMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} # compute terrain metrics f <- function(x, y, z) { as.list(lidaRtRee::points2terrainStats(data.frame(x, y, z))) } -terrain.metrics <- lidR::grid_metrics(point.cloud.ground, ~f(X, Y, Z), res=resolution) +metrics.terrain <- lidR::grid_metrics(point.cloud.ground, ~f(X, Y, Z), res=resolution) +par(mfrow=c(1,3)) +#for (i in c("altitude", "azimut.gr", "slope.gr")) +#{ +# raster::plot(metrics.terrain[[i]], main = i) +#} +raster::plot(metrics.terrain[["altitude"]], main = "altitude") +raster::plot(metrics.terrain$azimut.gr, main = "azimut (gr)", breaks = seq(from = 0, to = 400, length.out = 40), col=rainbow(40)) +raster::plot(metrics.terrain$slope.gr, main = "Slope (gr)") ``` -```{r old, include = FALSE} -# workflow associated with package lidaRtRee -# Copyright Irstea -# Author(s): Jean-Matthieu Monnet -# Licence: CeCill 2.1 -#######################" -# workflow for computation of ALS metrics on LAZ dataset -# application of prediction model to obtain forest parameter map -# inference function to compute estimates from whole map -# -# resolution for metrics map -resolution <- 25 -taille.dalle <- 1000 -################ -# metrics computation -# -# compute metrics if not already done -if (1)#!file.exists("./data/map.metrics.rda")) -{ - # create catalog of normalized LAZ files : should be rectangular tiles, with edges aligned on multiple values of the resolution - cata <- lidR::catalog("/media/reseau/services/LESSEM/ProjetsCommuns/Lidar/data/73-Bauges/norm.laz.pnr") - # cata <- lidR::catalog("/media/data/las/4montagnes/laznorm/") - library(foreach) - doFuture::registerDoFuture() - future::plan("multisession", workers = 2L) - # - map.metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar% +### Batch processing of tiles + +For the batch processing of tiles, parallel processing with packages `foreach` and `doFuture` is used. The processing capabilities of `lidR` are not yet used. + +```{r setupCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +library(foreach) +# create parallel frontend, specify to use two parallel sessions +doFuture::registerDoFuture() +future::plan("multisession", workers = 2L) +``` + +A buffering procedure is designed to handle the border effects when detecting trees at tile edges. A 10 m buffer is enough for tree metrics. One can also specify points classes to use or apply a threshold to high points. Some parameters specified in the previous paragraphs are also required for the batch processing (output map resolution, chm resolution for tree segmentation, functions for metrics computation). + +```{r paramCluster, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +# processing parameters +b.size <- 10 +# height threshold (m) for high points removal (points above this threshold are considered +# as outliers) +h.points <- 60 +# points classes to retain for analysis (vegetation + ground) +# ground should be 2 +class.points <- c(0, 1, 2, 3, 4, 5) +``` + +This first chunk of code computes tree and point metrics from the normalized ALS tiles. + +```{r batchProcessing, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +# processing by tile +metrics <- foreach::foreach (i=1:nrow(cata.height), .errorhandling="remove") %dopar% { - centroid.x <- mean(as.numeric(cata@data[i,c("Min.X", "Max.X")])) - centroid.y <- mean(as.numeric(cata@data[i,c("Min.Y", "Max.Y")])) - origine.x <- floor(centroid.x/taille.dalle)*taille.dalle - origine.y <- floor(centroid.y/taille.dalle)*taille.dalle - a <- lidR::readLAS(cata$filename[i]) - a <- lidR::lasfilter(a, X>=origine.x & X<origine.x+taille.dalle & Y>=origine.y & Y<origine.y+taille.dalle) + # tile extent + b.box <- as.numeric(cata.height@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")]) + # load tile extent plus buffer + a <- try(lidR::clip_rectangle(cata.height, b.box[1]-b.size, b.box[2]-b.size, b.box[3]+b.size, + b.box[4]+b.size)) + # + # check if points are successfully loaded + if (class(a)=="try-error") return(NULL) + # add 'buffer' flag to points in buffer with TRUE value in this new attribute + a <- lidR::add_attribute(a, + a$X < b.box[1] | a$Y < b.box[2] | a$X >= b.box[3] | a$Y >= b.box[4], + "buffer") + # remove unwanted point classes, and points higher than height threshold + a <- lidR::filter_poi(a, is.element(Classification, class.points) & Z<=h.points) + # check number of remaining points + if (length(a)==0) return(NULL) # set negative heights to 0 a@data$Z[a@data$Z<0] <- 0 - dummy <- lidR::grid_metrics(a, lidaRtRee::ABAmodelMetrics(Z, Intensity, ReturnNumber, Classification, 2), res=resolution) - - return(dummy) + # + # compute chm + chm <- lidaRtRee::points2DSM(a, res=aba.resCHM) + # replace NA, low and high values by 0 + chm[is.na(chm) | chm<0 | chm>h.points] <- 0 + # + # compute tree metrics + # tree top detection (default parameters) + segms <- lidaRtRee::treeSegmentation(chm) + # extraction to data.frame + trees <- lidaRtRee::treeExtraction(segms$filled.dem, segms$local.maxima, segms$segments.id) + # remove trees outside of tile + trees <- trees[trees$x>=b.box[1] & trees$x<b.box[3] & trees$y>=b.box[2] & trees$y<b.box[4],] + # compute raster metrics + metrics.trees <- lidaRtRee::rasterMetrics(trees[,-1], res=resolution, + fun=function(x){ + lidaRtRee::stdTreeMetrics(x, resolution^2/10000) + }, + output="raster") + # compute canopy cover in trees and canopy mean height in trees + # in region of interest, because it is not in previous step. + r.treechm <- segms$filled.dem + # set chm to NA in non segment area + r.treechm[segms$segments.id==0] <- NA + # compute raster metrics + metrics.trees2 <- lidaRtRee::rasterMetrics( + raster::crop(r.treechm, raster::extent(b.box[1], b.box[3],b.box[2], b.box[4])), + res=resolution, + fun=function(x){ + c(sum(!is.na(x$filled.dem))/(resolution/aba.resCHM)^2, + mean(x$filled.dem,na.rm=T)) + } + , output="raster") + names(metrics.trees2) <- c("TreeCanopy.coverInPlot", "TreeCanopy.meanHeightInPlot") + # + # compute 1D height metrics + # remove buffer points + a <- lidR::filter_poi(a, buffer==0) + # + if (length(a)==0) return(NULL) + # all points metrics + metrics.points <- lidR::grid_metrics(a, aba.pointMetricsFUN, res = resolution) + # + # extend / crop to match metrics.points + metrics.trees <- raster::extend(metrics.trees, metrics.points, values = NA) + metrics.trees2 <- raster::extend(metrics.trees2, metrics.points, values = NA) + metrics.trees <- raster::crop(metrics.trees, metrics.points) + metrics.trees2 <- raster::crop(metrics.trees2, metrics.points) + # merge rasterstacks + metrics <- metrics.points + metrics <- raster::addLayer(metrics, metrics.trees) + metrics <- raster::addLayer(metrics, metrics.trees2) + return(metrics) } - # rbind data.frame with option "fill" to handle data.frames with 0 rows and less columns - names.metrics <- names(map.metrics[[1]]) - test <- Reduce(raster::merge, map.metrics) - map.metrics <- do.call(raster::merge, map.metrics) - names(map.metrics) <- names.metrics +# mosaic rasters +names.metrics <- names(metrics[[1]]) +metrics.map <- do.call(raster::merge, metrics) +names(metrics.map) <- names.metrics +``` -} else {load("./data/map.metrics.rda")} + +The second chunk of code computes terrain metrics from the ALS tiles with altitude values. + +```{r batchProcessingTerrain, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +# function to compute terrain metrics to input to grid_metrics +f <- function(x, y, z) +{ + as.list(lidaRtRee::points2terrainStats(data.frame(x, y, z))) +} # -# ADD TERRAIN METRICS <- REQUIRES NON NORMALIZED POINT CLOUD -# ADD TREE METRICS <- REQUIRES ADDITIONNAL CODE TO COMPUTE THEM ON GRID -################# -# rasterize strata -public <- rgdal::readOGR("./data/", "Public4Montagnes") -map.metrics@crs <- public@proj4string -map.metrics$strata <- raster::rasterize(public, map.metrics) -map.metrics$strata[is.na(map.metrics$strata)] <- 0 -levels(map.metrics$strata) <- data.frame(ID=c(0,1), propriete=c("private", "public")) +metrics.terrain <- foreach::foreach (i=1:nrow(cata.altitude), .errorhandling="remove") %dopar% + { + # tile extent + b.box <- as.numeric(cata.altitude@data[i, c("Min.X", "Min.Y", "Max.X", "Max.Y")]) + # load tile extent plus buffer + a <- try(lidR::clip_rectangle(cata.altitude, b.box[1], b.box[2], b.box[3], + b.box[4])) + # check if points are successfully loaded + if (class(a)=="try-error") return(NULL) + # keep only ground points + a <- lidR::filter_ground(a) + lidR::grid_metrics(a, ~f(X, Y, Z), res=resolution) + } +# mosaic rasters +names.metrics <- names(metrics.terrain[[1]]) +metrics.terrain <- do.call(raster::merge, metrics.terrain) +names(metrics.terrain) <- names.metrics +``` + +Finally all metrics are combined in a single object. + +```{r allMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +metrics.terrain <- raster::extend(metrics.terrain, metrics.map, values = NA) +metrics.terrain <- raster::crop(metrics.terrain, metrics.map) +metrics.terrain <- raster::dropLayer(metrics.terrain, "adjR2.plane") +# merge rasterstacks +metrics.map <- raster::addLayer(metrics.map, metrics.terrain) +``` + +If the metrics maps are to be displayed in external software, the tif format can be used to produce one single with all bands. Meanwhile band names are not retained, so it is useful to save them in a separated R archive. + +```{r saveMetrics, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +metrics.names <- names(metrics.map) +save(metrics.names, file = "../data/aba.model/output/metrics.names.rda") +raster::writeRaster(metrics.map, file = "../data/aba.model/output/metrics.map.tif", overwrite = TRUE) +save(metrics.map, file = "../data/aba.model/output/metrics.map.rda") +``` + +## Forest parameters mapping + +### Single (not stratified) model + +Forest parameters are then mapped using the mapped metrics and the previously calibrated model with the function `lidaRtRee::ABApredict`. + +```{r mapForest, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +prediction.map <- lidaRtRee::ABApredict(model.ABA, metrics.map) +``` + +### Stratified models + +In case strata-specific models have been calibrated, it is necessary to add a layer containing the strata information in the metrics raster. For this, the existing vector layer of public forests has to be rasterized and added to the metrics raster. The levels in the layer metadata have to be filled, both in the ID (code, numeric) and stratum (label, character) fields. + +```{r rasterizeStrat, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +metrics.map$stratum <- raster::rasterize(public, metrics.map, field = 1) +# set to 0 in non-public forests areas +metrics.map$stratum[is.na(metrics.map$stratum)] <- 0 +# label levels in the raster metadata +levels(metrics.map$stratum) <- data.frame(ID=c(0,1), stratum=c("private", "public")) +# plot raster and vector +raster::plot(metrics.map$stratum, legend=FALSE, col=c("green", "blue"),xaxt='n', yaxt='n') +legend("bottom", legend=c("private","public"), fill=c("green","blue")) +# crop public vector +public.cropped <- sf::st_crop(public, raster::extent(metrics.map)) +plot(public.cropped, col = NA, add = TRUE) +``` + +```{r rasterizeStrat, include=TRUE, message=FALSE, warning=FALSE, fig.width = 12, fig.height = 4} +prediction.map.mixed <- lidaRtRee::ABApredict(model.ABA.stratified.mixed, metrics.map, stratum="stratum") +# verification par comparaison avec model prive / public +model.private <- list(model = model.ABA.stratified.mixed$model$private, stats = model.ABA.stratified.mixed$stats["private",]) +prediction.map.private <- lidaRtRee::ABApredict(model.private, metrics.map) +model.public <- list(model = model.ABA.stratified.mixed$model$public, stats = model.ABA.stratified.mixed$stats["public",]) +prediction.map.public <- lidaRtRee::ABApredict(model.public, metrics.map) + +``` + +```{r old, include = FALSE} + +} else {load("./data/map.metrics.rda")} + # ################ # requires previous calibration of models (script workflow.prediction.model.calibration.R) diff --git a/R/forest.structure.metrics.Rmd b/R/forest.structure.metrics.Rmd index 45788b3664d68066f329473725c225b23f83cbb4..f8d542f490ee012d77972e32f2ab7f554b2d5a34 100644 --- a/R/forest.structure.metrics.Rmd +++ b/R/forest.structure.metrics.Rmd @@ -467,7 +467,7 @@ metrics <- foreach::foreach (i=1:nrow(cata), .errorhandling="remove") %dopar% a@data$Z[a@data$Z<0] <- 0 # # compute chm - chm <- lidaRtRee::points2DSM(a,res=resCHM) + chm <- lidaRtRee::points2DSM(a, res=resCHM) # replace NA, low and high values chm[is.na(chm) | chm<0 | chm>h.points] <- 0 # ########################## diff --git a/data/aba.model/output/terrain.stats.rda b/data/aba.model/output/terrain.stats.rda deleted file mode 100644 index cf7d4ac75ae5aba08ce9a7cd16f39f7825bbc66d..0000000000000000000000000000000000000000 Binary files a/data/aba.model/output/terrain.stats.rda and /dev/null differ