-
-
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 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
@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.
maybe we should check this right now in the session.
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