Skip to content

Instantly share code, notes, and snippets.

@gisma
Last active January 15, 2022 18:49
Show Gist options
  • Save gisma/4172ef049b116abb1454233c8950d587 to your computer and use it in GitHub Desktop.
Save gisma/4172ef049b116abb1454233c8950d587 to your computer and use it in GitHub Desktop.
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"))
@GeoKL
Copy link

GeoKL commented Dec 5, 2021

Good evening,
trying to go thru this script I got into some problems:

In line 16, when running projection(mof100_ctg) <- epsg_number, I get the Error # Error: object 'epsg_number' not found.

Later on, in line 34, I have trouble running
dtm_knnidw_1m <- grid_terrain(mof100_ctg, res=1, algorithm = knnidw(k = 6L, p = 2))

The Output is:
Processing [--------] 0% (0/2) eta: ?sWarning: There were 1 degenerated ground points. Some X Y Z coordinates were repeated. They were removed.
Warning: There were 45 degenerated ground points. Some X Y coordinates were repeated but with different Z coordinates. min Z were retained.
Processing [========] 100% (2/2) eta: 0s
An error occurred when processing the chunk 2. Try to load this chunk with:
chunk <- readRDS("C:\Users\Kathi\AppData\Local\Temp\Rtmp6PGiVm/chunk2.rds")
las <- readLAS(chunk)
cannot allocate vector of size 1.8 Gb.

I'm sorry, I can't phrase actual questions about this as I am still working on understanding all of the scripts.
Best regards

@gisma
Copy link
Author

gisma commented Dec 7, 2021

@GeoKL sorry for delay (I got no notification you can activate this by using @gisma) as discussed you need to add the content of the other snippets as well as the typicall setup sequence to this script-snippet. Due to the missing EPSG I assume that this is not done correctly.
The warnings are ok.
However the line
chunk <- readRDS("C:\Users\Kathi\AppData\Local\Temp\Rtmp6PGiVm/chunk2.rds") is not correct.

Even if you are doing so (loading something called chunk2.rds) whichis not part of the script, you are additionally not following the setup rules. The addressed temp folder is a system generated irregular temporary temp folder (have a look at ?tempfile() to learn more), which is probably during your call not available any more (what results in the error).

Please follow exactly the suggested setup procedure add the missings parts that you find on the web page and try again

@GeoKL
Copy link

GeoKL commented Dec 7, 2021

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

@gisma
Copy link
Author

gisma commented Dec 7, 2021

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