diff --git a/R/ALS_data_preprocessing.Rmd b/R/ALS_data_preprocessing.Rmd index 7d985f953e4b341d8d0f5556c224ecba2c1bd1fd..36a770470830394b249dff9a1769a65b730c1c2a 100755 --- a/R/ALS_data_preprocessing.Rmd +++ b/R/ALS_data_preprocessing.Rmd @@ -81,11 +81,12 @@ ALS data are usually stored in files respecting the [ASPRS LAS specifications](h A single LAS file can be loaded in R with the `lidR::readLAS` function, which returns an object of class LAS. Some checks are automatically performed when the file is read. Function options allow to skip loading certain attributes, which might help saving time and memory. The object contains several slots, including: * `header`: metadata read from the file, -* `data`: point cloud attributes, -* `bbox`: the bounding box of the data, -* `proj4string`: projection information. +* `data`: point cloud attributes. +* `crs`: projection information. -It is advisable to fill the projection information if missing, as spatial objects computed from the point cloud will inherit the information. A warning is issued when loading the file because unexpected values are present in the scan angle field. +The bounding box of the data cona be accessed with the function `sf::st_bbox`. + +It is advisable to fill the projection information if missing, as spatial objects computed from the point cloud will inherit the information. A warning is issued when loading this file because unexpected values are present in the scan angle field. ```{r ALS.load, include = TRUE} # read file @@ -105,7 +106,7 @@ point_cloud@header ``` ### Basic attributes -The `data` slot contains point attributes: +The `data` slot contains point attributes. Some attributes may be present or not or have different names depending on the format version of the LAS file. The main attributes (coordinates, intensity, classification, return number) are always present. **`X`, `Y`, `Z`**: coordinates. The point cloud can be displayed with `lidR::plot` ```{r ALS.coordinates, eval=FALSE} @@ -262,24 +263,24 @@ f <- nstrip = length(unique(pid))) # number of unique values of strips } # apply the function to the catalog, indicating which attributes to use, and output resolution -metrics <- lidR::grid_metrics(cata, ~f(ReturnNumber, Classification, PointSourceID), res = res) +metrics <- lidR::pixel_metrics(cata, ~f(ReturnNumber, Classification, PointSourceID), res = res) # convert number to density in first three layers for (i in c("npoints", "npulses", "nground")) { metrics[[i]] <- metrics[[i]] / res ^ 2 } -raster::plot(metrics) +terra::plot(metrics) # percentage of 5 m cells with no ground points -p_no_ground <- round(sum(raster::values(metrics$nground)==0) / length(raster::values(metrics$nground)) * 100, 1) +p_no_ground <- round(sum(terra::values(metrics$nground)==0) / length(terra::values(metrics$nground)) * 100, 1) # percentage of 5 m cells with less than 5 pulses / m2 -p_pulse_5 <- round(sum(raster::values(metrics$npulses)<10) / length(raster::values(metrics$npulses)) * 100, 1) +p_pulse_5 <- round(sum(terra::values(metrics$npulses)<10) / length(terra::values(metrics$npulses)) * 100, 1) ``` From those maps summary statistics can be derived, e.g. the percentage of cells with no ground points (`r p_no_ground`). Checking that the spatial distribution of pulses is homogeneous can be informative. ```{r catalog.maps, fig.width = 9, fig.height = 4} par(mfrow = c(1,2)) -hist(raster::values(metrics$npulses), main = "Histogram of pulse density", xlab = "Pulse density /m2", ylab = "Number of 5m cells") -raster::plot(metrics$npulses < 5, main = "Areas with pulse density < 5 /m2") +hist(terra::values(metrics$npulses), main = "Histogram of pulse density", xlab = "Pulse density /m2", ylab = "Number of 5m cells") +terra::plot(metrics$npulses < 5, main = "Areas with pulse density < 5 /m2") ``` @@ -313,26 +314,26 @@ It represents the altitude of the bare earth surface. It is computed using point ```{r dem.terrain, include = TRUE, message = FALSE} # create dtm at 0.5 m resolution with TIN algorithm -dtm <- lidR::grid_terrain(cata, res = 0.5, lidR::tin()) +dtm <- lidR::rasterize_terrain(cata, res = 0.5, lidR::tin()) dtm ``` -Functions from the `raster` package are available to derive topographical rasters (slope, aspect, hillshade, contour lines) from the DTM. +Functions from the `terra` package are available to derive topographical rasters (slope, aspect, hillshade, contour lines) from the DTM. ```{r dem.terrain.derived, include = TRUE, message = FALSE, fig.dim = c(6, 4)} -dtm_slope <- raster::terrain(dtm) -dtm_aspect <- raster::terrain(dtm, opt = "aspect") -dtm_hillshade <- raster::hillShade(dtm_slope, dtm_aspect) -raster::plot(dtm_hillshade, col = gray(seq(from = 0, to = 1, by = 0.01)), main = "Hillshade") +dtm_slope <- terra::terrain(dtm, unit = "radians") +dtm_aspect <- terra::terrain(dtm, v = "aspect", unit = "radians") +dtm_hillshade <- terra::shade(dtm_slope, dtm_aspect) +terra::plot(dtm_hillshade, col = gray(seq(from = 0, to = 1, by = 0.01)), main = "Hillshade") ``` ### Digital surface model (DSM) -It represents the altitude of the earth surface, including natural and man-made objects. It is computed using all points. The function `lidR::grid_canopy` proposes different algorithms for that purpose, and works either with catalogs or LAS objects. A straightforward way to compute the DSM is to retain the value of the highest point present in each pixel. In this case, choosing an adequate resolution with respect to the point density should limit the proportion of empty cells. +It represents the altitude of the earth surface, including natural and man-made objects. It is computed using all points. The function `lidR::rasterize_canopy` proposes different algorithms for that purpose, and works either with catalogs or LAS objects. A straightforward way to compute the DSM is to retain the value of the highest point present in each pixel. In this case, choosing an adequate resolution with respect to the point density should limit the proportion of empty cells. ```{r dem.surface, include = TRUE, message = FALSE} # create dsm at 0.5 m resolution with highest point per pixel algorithm -dsm <- lidR::grid_canopy(cata, res = 0.5, lidR::p2r()) +dsm <- lidR::rasterize_canopy(cata, res = 0.5, lidR::p2r()) dsm ``` @@ -348,9 +349,9 @@ chm cata_norm <- lidR::catalog("../data/aba.model/ALS/tiles.norm.laz/") # set projection lidR::projection(cata_norm) <- 2154 -chm_norm <- lidR::grid_canopy(cata_norm, res = 0.5, lidR::p2r()) +chm_norm <- lidR::rasterize_canopy(cata_norm, res = 0.5, lidR::p2r()) chm_norm -# boxplot(raster::values(chm - chm_norm)) +# boxplot(terra::values(chm - chm_norm)) ``` If high points are present in the point cloud, a threshold can be applied to the CHM values to remove outliers. @@ -358,7 +359,7 @@ If high points are present in the point cloud, a threshold can be applied to the ```{r dem.chm.threshold, include = TRUE, message = FALSE, fig.dim = c(6, 4)} threshold <- 40 chm[chm > threshold] <- threshold -raster::plot(chm, main = "Canopy Height Model (m)") +terra::plot(chm, main = "Canopy Height Model (m)") ``` ## Batch processing of files @@ -380,7 +381,7 @@ lidR::opt_chunk_buffer(cata) <- 10 # DTM output file template lidR::opt_output_files(cata) <- "dtm_{ORIGINALFILENAME}" # create DTM -dtm <- lidR::grid_terrain(cata, 0.5, lidR::tin()) +dtm <- lidR::rasterize_terrain(cata, 0.5, lidR::tin()) # # laz compression lidR::opt_laz_compression(cata) <- TRUE @@ -394,7 +395,7 @@ lidR::opt_chunk_buffer(cata) <- 0 # DSM output file template lidR::opt_output_files(cata) <- "dsm_{ORIGINALFILENAME}" # create DSM -dsm <- lidR::grid_canopy(cata, 0.5, lidR::p2r()) +dsm <- lidR::rasterize_canopy(cata, 0.5, lidR::p2r()) # # create CHM from normalized files # process by file @@ -404,7 +405,7 @@ lidR::opt_chunk_buffer(cata_norm) <- 0 # CHM output file template lidR::opt_output_files(cata_norm) <- "chm_{ORIGINALFILENAME}" # create CHM -chm <- lidR::grid_canopy(cata_norm, 0.5, lidR::p2r()) +chm <- lidR::rasterize_canopy(cata_norm, 0.5, lidR::p2r()) ``` ## References \ No newline at end of file diff --git a/R/area-based.1.data.preparation.Rmd b/R/area-based.1.data.preparation.Rmd index f0cc4b84b6d8e7c8866f4a82e40435557d8df807..f8369857bcd2b1c0c79bc31bd3c2a4647b3c05c2 100755 --- a/R/area-based.1.data.preparation.Rmd +++ b/R/area-based.1.data.preparation.Rmd @@ -325,7 +325,7 @@ plot(slope_gr ~ azimuth_gr, ) ``` -The following lines create a spatial polygon corresponding to one plot extent. +The following lines create a `sf` spatial polygon corresponding to one plot extent. ```{r plotExtent, include = TRUE} # extract example plot @@ -440,9 +440,9 @@ The tree positions can be plotted with the canopy height model computed with the ```{r plotRS, include = TRUE, out.width = '70%', fig.dim=c(6.5, 4.5), warning=FALSE} # compute background # compute canopy height model -chm <- lidaRtRee::points2DSM(llas_height[[plot_test]], res = 0.5) +chm <- lidR::rasterize_canopy(llas_height[[plot_test]], res = 0.5, algorithm = lidR::p2r(), pkg = "terra") # display CHM -raster::plot(chm, +terra::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "Canopy Height Model and tree positions" ) @@ -461,7 +461,7 @@ Similar output with `ggplot2`. ```{r plotRS.ggplot2, include = TRUE, out.width = '70%', fig.dim=c(6.5, 4.5), warning=FALSE} # convert raster to data.frame for display with ggplot -chm.df <- raster::as.data.frame(chm, xy = TRUE) +chm.df <- terra::as.data.frame(chm, xy = TRUE) # rename columns names(chm.df) <- c("x", "y", "CHM") # create graph diff --git a/R/area-based.2.model.calibration.Rmd b/R/area-based.2.model.calibration.Rmd index 3982cb1989d5fe258b0b1a8058ae6d7be9656256..816eb5789b6a905a5f44bdd045324dcd82622a63 100755 --- a/R/area-based.2.model.calibration.Rmd +++ b/R/area-based.2.model.calibration.Rmd @@ -178,9 +178,9 @@ lidaRtRee::aba_plot(model_aba, 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")], +plot(model_aba$values$predicted,# plots[subsample, c("G_m2_ha")], model_aba$values$residual, - ylab = "Prediction errors", xlab = "Field values", + ylab = "Prediction errors", xlab = "Predicted in LOOCV", col = ifelse(plots$stratum == "public", "green", "blue") ) abline(h = 0, lty = 2) diff --git a/R/tree.detection.Rmd b/R/tree.detection.Rmd index 97765bda3f0f0a09883a92e16979dc48805e5afc..9c21b499ea5a6c90253a581e899b9d1d6d1bb5aa 100755 --- a/R/tree.detection.Rmd +++ b/R/tree.detection.Rmd @@ -112,6 +112,8 @@ In this tutorial, ALS data available in the `lidaRtRee` package is used. # load data in package lidaRtRee (default) LASfile <- system.file("extdata", "las_chablais3.laz", package="lidaRtRee") las_chablais3 <- lidR::readLAS(LASfile) +# set projection +lidR::projection(las_chablais3) <- 2154 ``` Otherwise, ALS data is loaded with functions of package `lidR`. First a catalog of files containing ALS data is built. Then points located inside our ROI are loaded. @@ -143,37 +145,31 @@ From the ALS point cloud, digital elevation models are computed [@Monnet11c, pp. ```{r computeDEMs, include=TRUE, warning = FALSE} # define extent and resolution of raster -output_raster <- raster::raster(resolution = 0.5, - xmn = ROI_range$x[1] - ROI_buffer, - xmx = ROI_range$x[2] + ROI_buffer, - ymn = ROI_range$y[1] - ROI_buffer, - ymx = ROI_range$y[2] + ROI_buffer - ) +output_raster <- terra::rast(resolution = 0.5, + xmin = ROI_range$x[1] - ROI_buffer, + xmax = ROI_range$x[2] + ROI_buffer, + ymin = ROI_range$y[1] - ROI_buffer, + ymax = ROI_range$y[2] + ROI_buffer, + crs = sf::st_crs(las_chablais3)$wkt +) # terrain model computed from points classified as ground -dtm <- lidR::grid_terrain(las_chablais3, output_raster, lidR::tin()) -# dtm_lidaRtRee <- lidaRtRee::points2DTM(lidR::filter_ground(las_chablais3), -# res = 0.5, ROI_range$x[1] - ROI_buffer, -# ROI_range$x[2] + ROI_buffer, ROI_range$y[1] - ROI_buffer, -# ROI_range$y[2] + ROI_buffer -# ) +dtm <- lidR::rasterize_terrain(las_chablais3, output_raster, lidR::tin()) # surface model -dsm <- lidR::grid_canopy(lidR::clip_roi(las_chablais3, raster::extent(output_raster)), output_raster, lidR::p2r()) -# dsm_lidaRtRee <- lidaRtRee::points2DSM(las_chablais3, -# res = 0.5, dtm@extent@xmin, -# dtm@extent@xmax, dtm@extent@ymin, dtm@extent@ymax -# ) +dsm <- lidR::rasterize_canopy(las_chablais3, output_raster, lidR::p2r()) # canopy height model chm <- dsm - dtm +# save for later use +# chm_chablais3 <- terra::wrap(chm); save(chm_chablais3, file = "~/R/lidaRtRee/data/chm_chablais3.rda") ``` -```{r plotDEMs, echo=FALSE, fig.width = 12, fig.height = 4.5, out.width='100%'} +```{r plotDEMs, echo=FALSE, fig.width = 12, fig.height = 3.5, out.width='100%'} par(mfrow = c(1, 3)) # display DTM -raster::plot(dtm, main = "DTM") +terra::plot(dtm, main = "DTM") # display DSM -raster::plot(dsm, main = "DSM") +terra::plot(dsm, main = "DSM") # display CHM -raster::plot(chm, main = "CHM") +terra::plot(chm, main = "CHM") ``` ### Visual comparison of field inventory and ALS data @@ -191,14 +187,16 @@ mask_tree <- lidaRtRee::raster_xy_mask( ) # union of convexHull and tree buffers mask_plot <- max(mask_chull, mask_tree) +# set zeros values to NA +mask_plot[mask_plot==0] <- NA # vectorize plot mask -mask_plot_v <- raster::rasterToPolygons(mask_plot, function(x) (x == 1), dissolve = T) +mask_plot_v <- terra::as.polygons(mask_plot) ``` Displaying inventoried trees on the CHM shows a pretty good agreement of crowns visible in the CHM with trunk locations and sizes. ```{r plotPlot, include = TRUE, out.width = '70%', fig.dim=c(6.5, 4.5), warnings=FALSE} # display CHM -raster::plot(chm, +terra::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "Canopy Height Model and tree positions" ) @@ -208,7 +206,7 @@ lidaRtRee::plot_tree_inventory(tree_inventory_chablais3[, c("x", "y")], species = as.character(tree_inventory_chablais3$s), add = TRUE ) # display plot mask -raster::plot(mask_plot_v, border = "red", add = TRUE) +terra::plot(mask_plot_v, border = "red", add = TRUE) ``` ## Tree delineation @@ -220,19 +218,19 @@ Tree segmentation is performed on the Canopy Height Model by using a general fun + image segmentation with a watershed algorithm for crown delineation. The first two steps are documented in @Monnet11c, pp. 47-52. -```{r plotSegm, include=TRUE, fig.width = 12, fig.height = 4.5, out.width='100%'} +```{r plotSegm, include=TRUE, fig.width = 12, fig.height = 3.5, out.width='100%'} # tree detection (default settings), applied on canopy height model segms <- lidaRtRee::tree_segmentation(chm) # par(mfrow = c(1, 3)) # display pre-processed chm -raster::plot(segms$smoothed_dem, main = "Pre-processed CHM") +terra::plot(segms$smoothed_dem, main = "Pre-processed CHM") # display selected local maxima -raster::plot(segms$local_maxima, main = "Selected local maxima") +terra::plot(segms$local_maxima, main = "Selected local maxima") # display segments, except ground segment dummy <- segms$segments_id dummy[dummy == 0] <- NA -raster::plot(dummy %% 8, main = "Segments (random colors)", col = rainbow(8), legend = FALSE) +terra::plot(dummy %% 8, main = "Segments (random colors)", col = rainbow(8)) ``` ### Extraction of apices positions and attributes @@ -258,16 +256,16 @@ apices <- lidaRtRee::tree_extraction( ) head(apices, n = 3L) # convert segments from raster to polygons -segments_v <- raster::rasterToPolygons(segms$segments_id, dissolve = T) +segments_v <- terra::as.polygons(segms$segments_id) # # display initial image raster::plot(chm, col = gray(seq(0, 1, 1 / 255)), main = "CHM and detected positions") # display segments border -sp::plot(segments_v, border = "white", add = T) +terra::plot(segments_v, border = "white", add = T) # display plot mask -sp::plot(mask_plot_v, border = "red", add = T) +terra::plot(mask_plot_v, border = "red", add = T) # display detected apices -graphics::points(apices$x, apices$y, col = "blue", cex = apices$h / 20, pch = 2) +plot(apices["h"], col = "blue", cex = apices$h / 20, pch = 2, add = TRUE) ``` ## Detection evaluation @@ -279,12 +277,12 @@ To assess detection accuracy, reference (field) trees should be linked to detect # match detected apices with field trees based on relative distance of apices matched <- lidaRtRee::tree_matching( tree_inventory_chablais3[, c("x", "y", "h")], - cbind(apices@coords, apices$h) + cbind(sf::st_coordinates(apices), apices$h) ) # display matching results lidaRtRee::plot_matched( tree_inventory_chablais3[, c("x", "y", "h")], - cbind(apices@coords, apices$h), matched, chm, mask_plot_v + cbind(sf::st_coordinates(apices), apices$h), matched, chm, mask_plot_v ) ``` @@ -294,7 +292,7 @@ lidaRtRee::plot_matched( # height histogram of detections detection_stats <- lidaRtRee::hist_detection( tree_inventory_chablais3[, c("x", "y", "h")], - cbind(apices@coords, apices$h), matched + cbind(sf::st_coordinates(apices), apices$h), matched ) ``` @@ -304,7 +302,7 @@ Detection accuracy is evaluated thanks to the number of correct detections (`r d # height histogram of detections detection_stats <- lidaRtRee::hist_detection( tree_inventory_chablais3[, c("x", "y", "h")], - cbind(apices@coords, apices$h), matched + cbind(sf::st_coordinates(apices), apices$h), matched ) ``` @@ -312,7 +310,7 @@ detection_stats <- lidaRtRee::hist_detection( ```{r heightRegression, include=FALSE} height_reg <- lidaRtRee::height_regression(tree_inventory_chablais3[, c("x", "y", "h")], - cbind(apices@coords, apices$h), + cbind(sf::st_coordinates(apices), apices$h), matched, species = tree_inventory_chablais3$s ) @@ -323,7 +321,7 @@ For true detections, estimated heights can be compared to field-measured heights ```{r plotRegression, include=TRUE, out.width = '40%', fig.dim=c(4.5, 4.5)} # linear regression between reference height and estimated height height_reg <- lidaRtRee::height_regression(tree_inventory_chablais3[, c("x", "y", "h")], - cbind(apices@coords, apices$h), + cbind(sf::st_coordinates(apices), apices$h), matched, species = tree_inventory_chablais3$s ) @@ -374,10 +372,11 @@ head(metrics, n = 3L) Computed metrics are merged with reference trees and detected apices, thanks to the segment id. ```{r mergeMetrics, include=TRUE, eval=TRUE} # associate each reference tree with the segment that contains its trunk. -tree_inventory_chablais3$seg_id <- raster::extract( +dummy <- terra::extract( segms$segments_id, tree_inventory_chablais3[, c("x", "y")] ) +tree_inventory_chablais3$seg_id <- dummy$segments_id # create new data.frame by merging metrics and inventoried trees based on segment id tree_metrics <- base::merge(tree_inventory_chablais3, metrics) # remove non-tree segment @@ -532,23 +531,25 @@ To display the classification results, an image of the classified segments is cr # create image of predicted species species <- segms$segments_id # replace segment id by predicted species in image -raster::values(species) <- +terra::values(species) <- metrics$predicted_s[match(raster::values(segms$segments_id), metrics$seg_id)] # remove ground segment species[segms$segments_id == 0] <- NA # convert to factor raster -species <- raster::as.factor(species) +species <- terra::as.factor(species) # build raster attribute table rat -rat <- raster::levels(species)[[1]] -rat$Species <- levels(metrics$predicted_s) -levels(species)[[1]] <- rat +rat <- data.frame(id = terra::levels(species)[[1]], + Species = levels(metrics$predicted_s)) + # retrieve reference colors -target_col <- lidaRtRee::species_color()[rat$Species, "col"] +rat$col <- lidaRtRee::species_color()[rat$Species, "col"] # set NA color to green -target_col[is.na(target_col)] <- "green" +rat$col[is.na(rat$col)] <- "green" +# +levels(species)[[1]] <- rat # display results -raster::plot(species, col = target_col, legend = FALSE) -legend("bottomright", legend=rat$Species, fill=target_col) +terra::plot(species, col = rat$col, legend = FALSE) +legend("bottomright", legend=rat$Species, fill=rat$col) sp::plot(segments_v, add = TRUE, border = "white") lidaRtRee::plot_tree_inventory(tree_metrics_h [, c("x", "y")], @@ -561,12 +562,14 @@ lidaRtRee::plot_tree_inventory(tree_metrics_h ```{r saveGIS, include=FALSE, warning=FALSE, message=FALSE} # # extract all apices # apices_all <- lidaRtRee::tree_extraction(segms$filled_dem, segms$local_maxima, segms$segments_id) +# # round height values to cm +# apices_all$h <- round(apices_all$h, 2) # # save outputs -# raster::writeRaster(chm, file = "./data/output/chm.tif", overwrite = TRUE) -# raster::writeRaster(species, file = "./data/output/r_species.tif", overwrite = TRUE) -# raster::writeRaster(segms$segments_id, file = "./data/output/r_segments.tif", overwrite = TRUE) -# rgdal::writeOGR(segments_v, "./data/output/", "v_segments", driver = "ESRI Shapefile", overwrite = TRUE) -# write.csv(apices_all, "./data/output/apices.csv", row.names = FALSE) +# terra::writeRaster(chm, file = "./data/output/chm.tif", overwrite = TRUE) +# terra::writeRaster(species, file = "./data/output/r_species.tif", overwrite = TRUE) +# terra::writeRaster(segms$segments_id, file = "./data/output/r_segments.tif", overwrite = TRUE) +# terra::writeVector(segments_v, "./data/output/v_segments.gpkg", overwrite = TRUE) +# sf::st_write(apices_all, "./data/output/apices.gpkg", delete_layer = TRUE) # write.csv(tree_metrics_h, "./data/output/correct_detections.csv", row.names = FALSE) # write.csv(tree_inventory_chablais3, "./data/output/tree_inventory_chablais.csv", row.names = FALSE) ``` @@ -662,6 +665,7 @@ for (i in 1:nrow(tree_inventory_chablais3)) ```{r plotTreeModel2, include=FALSE, eval=FALSE} # virtual apices from detection +coord <- data.frame(sf::st_coordinates(apices)) # shape different based on mean lidar intensity value (threshold 55) library(rLiDAR) rgl::rgl.open() @@ -676,8 +680,8 @@ for (i in 1:nrow(apices)) HCB = 0.4 * apices$h[i], dbh = apices$h[i] / 50, resolution = "high", - X = apices$x[i], - Y = apices$y[i], + X = coord$X[i], + Y = coord$Y[i], mesh = F ) } else { @@ -688,8 +692,8 @@ for (i in 1:nrow(apices)) HCB = 0.5 * apices$h[i], dbh = apices$h[i] / 50, resolution = "high", - X = apices$x[i], - Y = apices$y[i], + X = coord$X[i], + Y = coord$Y[i], mesh = F, stemcolor = "burlywood4", crowncolor = "darkgreen" @@ -766,17 +770,18 @@ future::plan("multisession", workers = 2L) options(future.rng.onMisuse = "ignore") # # COORDINATES OF TILES +bbox <- sf::st_bbox(cata) # low left corner x <- seq( - from = floor(cata@bbox["x", "min"] / tile_size) * tile_size, + from = floor(bbox$xmin / tile_size) * tile_size, by = tile_size, - to = ceiling(cata@bbox["x", "max"] / tile_size - 1) * tile_size + to = ceiling(bbox$xmax / tile_size - 1) * tile_size ) # [1:2] length_x <- length(x) y <- seq( - from = floor(cata@bbox["y", "min"] / tile_size) * tile_size, + from = floor(bbox$ymin / tile_size) * tile_size, by = tile_size, - to = ceiling(cata@bbox["y", "max"] / tile_size - 1) * tile_size + to = ceiling(bbox$ymax / tile_size - 1) * tile_size ) # [1:2] length_y <- length(y) # repeat coordinates for tiling while area @@ -800,26 +805,26 @@ resultats <- future.apply::future_mapply( j + tile_size + buffer_size ) # check if points are present - if (nrow(las_tile@data) > 0) { + if (nrow(las_tile) > 0) { # normalization if required # las_tile <- lidR::normalize_height(las_tile, lidR::tin()) # in this example LAS tiles are already normalized # compute canopy height model - chm <- lidaRtRee::points2DSM(las_tile) + chm <- lidR::rasterize_canopy(las_tile, resolution, algorithm = lidR::p2r(), pkg = "terra") # check all chm is not NA - if (!all(is.na(raster::values(chm)))) { + if (!all(is.na(terra::values(chm)))) { # output chm if asked for if (out_chm) { - output$chm <- raster::crop(chm, raster::extent( + output$chm <- terra::wrap(terra::crop(chm, terra::ext( i, i + tile_size, j, j + tile_size - )) + ))) } # save on disk if (out_save_chm) { - raster::writeRaster(raster::crop( + terra::writeRaster(terra::crop( chm, - raster::extent( + terra::ext( i, i + tile_size, j, j + tile_size ) @@ -838,21 +843,27 @@ resultats <- future.apply::future_mapply( segms$segments_id ) # remove apices in buffer area - apices <- apices[apices$x >= i & apices$x < i + tile_size & - apices$y >= j & apices$y < j + tile_size, ] + # copy coordinates + apices <- cbind(apices, sf::st_coordinates(apices)) + apices <- apices[apices$X >= i & apices$X < i + tile_size & + apices$Y >= j & apices$Y < j + tile_size, ] # add tile id to apices to avoid duplicates in final file apices$tile <- paste0(i, "_", j) # convert to vectors if option is TRUE if (out_vectorize_apices) { # vectorize - apices_v <- raster::rasterToPolygons(segms$segments_id, dissolve = T) + apices_v <- terra::as.polygons(segms$segments_id) # remove polygons which treetop is in buffer apices_v <- apices_v[is.element(apices_v$segments_id, apices$id), ] # names(apices_v) <- "id" # add attributes # errors when using sp::merge so using sp::over even if it is probably slower + # convert to sf + apices_v <- sf::st_as_sf(apices_v) + names(apices_v)[1] <- "id" + apices_v <- merge(apices_v, sf::st_drop_geometry(apices)) # merge(apices_v@data, apices, all.x = TRUE) - apices_v@data <- cbind(apices_v@data, sp::over(apices_v, apices)) + #apices_v@data <- cbind(apices_v@data, sp::over(apices_v, apices)) # save in list output$apices_v <- apices_v } @@ -877,10 +888,10 @@ apices <- do.call(rbind, apices) # # chm if (out_chm) { - chm <- lapply(resultats, function(x) x[["chm"]]) + chm <- lapply(resultats, function(x) terra::rast(x[["chm"]])) # merge chm # no names in list otherwise do.call returns an error - chm_all <- do.call(raster::merge, chm) + chm_all <- do.call(terra::merge, chm) names(chm) <- id } # apices_v @@ -901,24 +912,24 @@ The following image displays the results for the whole area. chm_all[chm_all > 40] <- 40 chm_all[chm_all < 0] <- 0 # display chm -raster::plot(chm_all, +terra::plot(chm_all, main = "Canopy Height Model and segments" ) # display segments border -sp::plot(apices_v, border = "white", add = T) +plot(sf::st_geometry(apices_v), border = "white", add = T) # add apices -sp::plot(apices, cex = apices$h / 40, add = TRUE, pch = 2) +plot(sf::st_geometry(apices), cex = apices$h / 40, add = TRUE, pch = 2) ``` The following lines save outputs to files. ```{r batch.export, eval = FALSE} # merged chm -raster::writeRaster(chm_all, file = "chm.tif", overwrite = TRUE) +terra::writeRaster(chm_all, file = "chm.tif", overwrite = TRUE) # apices -sf::st_write(sf::st_as_sf(apices), 'apices_points.gpkg', 'apices_points', delete_dsn = TRUE) +sf::st_write(apices, "apices_points.gpkg", delete_layer = TRUE) # vectorized apices -if (out_vectorize_apices) sf::st_write(sf::st_as_sf(apices_v), 'v_apices_points.gpkg', 'v_apices_points', delete_dsn = T) +if (out_vectorize_apices) sf::st_write(apices_v, "v_apices_points.gpkg", delete_layer = TRUE) ``` ## References \ No newline at end of file diff --git a/data/aba.model/output/llas_height.rda b/data/aba.model/output/llas_height.rda index d73258906897c47ad80c009aca0f22169cb17e2d..48ca9e4c4ef7559a3e2803ee8734b4f2deb965f7 100755 Binary files a/data/aba.model/output/llas_height.rda and b/data/aba.model/output/llas_height.rda differ diff --git a/data/aba.model/output/metrics_terrain.rda b/data/aba.model/output/metrics_terrain.rda index 00b3dd382e01246c099a03d33eb0a95cc06e8b41..932e5985e9d7a160076960df5cdc69f041ad32cc 100755 Binary files a/data/aba.model/output/metrics_terrain.rda and b/data/aba.model/output/metrics_terrain.rda differ diff --git a/data/aba.model/output/models.rda b/data/aba.model/output/models.rda index 251da88737575c2f098eee1c33158f67489315ac..38983e1e9275261d34f658e9f77b0c8983ae8b09 100755 Binary files a/data/aba.model/output/models.rda and b/data/aba.model/output/models.rda differ diff --git a/data/aba.model/output/plots.rda b/data/aba.model/output/plots.rda index c25408484b357010b3c48ceec9950b9562782591..d45f374f1df3835fcdfcb9a531fa8f83804b8930 100755 Binary files a/data/aba.model/output/plots.rda and b/data/aba.model/output/plots.rda differ