catalog_make_CHM_Tiles.R

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.