-
-
Save gisma/4172ef049b116abb1454233c8950d587 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
set.seed(1000) | |
knitr::knit_global() | |
require(envimaR) | |
# MANDANTORY: defining the root folder DO NOT change this line | |
root_folder = "~/edu/agis" | |
#-- Further customization of the setup by the user this section | |
#-- can be freely customized only the definition of additional packages | |
#-- and directory paths MUST be done using the two variables | |
#-- appendpackagesToLoad and appendProjectDirList | |
#-- feel free to remove this lines if you do not need them | |
# define additional packages comment if not needed | |
appendpackagesToLoad = c("lidR","future","lwgeom") | |
# define additional subfolders comment if not needed | |
appendProjectDirList = c("data/lidar/","data/lidar/l_raster/","data/lidar/l_raw/","data/lidar/l_norm/") | |
## define current projection (It is not magic you need to check the meta data or ask your instructor) | |
## ETRS89 / UTM zone 32N | |
proj4 = "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0" | |
epsg_number = 25832 | |
lidR.progress = FALSE | |
# MANDANTORY: calling the setup script also DO NOT change this line | |
source(file.path(envimaR::alternativeEnvi(root_folder = root_folder),"src/000-rspatial-setup.R"),local = knitr::knit_global()) | |
# 1 - start script | |
#----------------------------- | |
#--- create path link to the original las file | |
utils::download.file(url="https://github.com/gisma/gismaData/raw/master/uavRst/data/lidR_data.zip", | |
destfile=paste0(envrmt$path_tmp,"/chm.zip")) | |
unzip(paste0(envrmt$path_tmp,"/chm.zip"), | |
exdir = envrmt$path_tmp, | |
overwrite = TRUE) | |
las_files = list.files(envrmt$path_tmp, | |
pattern = glob2rx("*.las"), | |
full.names = TRUE) | |
#--- create CHM with lidR catalog | |
#- source: https://github.com/Jean-Romain/lidR/wiki/Rasterizing-perfect-canopy-height-models | |
future::plan(future::multisession) | |
ctg <- lidR::readLAScatalog(las_files[[1]]) | |
lidR::projection(ctg) <- 25832 | |
lidR::opt_chunk_size(ctg) = 100 | |
lidR::opt_chunk_buffer(ctg) <- 5 | |
lidR::opt_progress(ctg) <- FALSE | |
lidR::opt_laz_compression(ctg) <- TRUE | |
ctg@output_options$drivers$Raster$param$overwrite <- TRUE | |
#--- height normalisation within the point cloud | |
lidR::opt_output_files(ctg) <- paste0(envrmt$path_l_norm,"/{ID}","_norm_height") | |
norm_ctg <- lidR::normalize_height(ctg,lidR::tin()) | |
# check new ctg | |
las_check(norm_ctg) | |
# calculate chm | |
chm <- grid_canopy(norm_ctg, res = 1.0, lidR::pitfree(thresholds = c(0,2,5,10,15), c(0,1.5))) | |
chm = get_vrt_img("chm",envrmt$path_l_norm,"_norm_height") | |
# check new ctg | |
las_check(norm_ctg) | |
#- mapview plot ctg | |
mapview::mapview(norm_ctg, zcol="Max.Z") | |
#- mapview plot | |
mapview::mapview(chm[[1]],layer.name = "pitfree chm 1 m² cells height [m]",col.regions=lidR::height.colors(20)) | |
#- classic plot | |
plot( chm, | |
col = lidR::height.colors(20), | |
main = "pitfree chm 1 m² cells", | |
cex.main = 1) | |
#- tmap plot | |
tm_shape(chm) + | |
tm_raster( title = "pitfree chm 1 m² cells", | |
palette = lidR::height.colors(20)) + | |
tm_grid()+ | |
tm_layout(legend.outside = TRUE) | |
#- ggplot plot with stars | |
ggplot() + | |
geom_stars(data = stars::st_as_stars(chm)) + | |
scale_fill_gradientn(colours=lidR::height.colors(20)) + | |
coord_equal()+ | |
guides(fill=guide_legend(title="pitfree chm 1 m² cells")) | |
maybe we should check this right now in the session.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
@gisma no problem.
So on the web page (https://gisma-courses.github.io/gi-modules/post/2021-11-20-forest-information-from-lidar-data/) it said "Take the setup section of the above script and exchange the code part as found under #— start script with the following code chunk.", so I took the setup from above (on the web page) and exchanged everything from "#— start script" on with this script here...
Since "#— start script" comes right after the setup and there are no more snipplets on the web page, I don't exactly understand what other missing parts you mean?
One thing I changed is that I added ".libPaths("E:/edu/Packages")" on the very top of the script (what I do in every script I use to save packages on my external drive).#
Also I added some more packages to appendpackagesToLoad and changed the appendProjectDirList to appendProjectDirList = c("data/lidar/","data/lidar/l_raster","data/l_raw","data/l_norm") since I wanted to have the l_raw and l_norm folder in my data folder. This didn't change the evrmt link though.