library(terra)
terra 1.7.78
Attaching package: ‘terra’
The following objects are masked from ‘package:lidR’:
area, crs, crs<-, is.empty
# this is my workplace
setwd("~/Library/CloudStorage/OneDrive-UniversityofArizona/Documents/sandbox/Reproducibility/remote_sensing_fp/3_dsm_ortho/1_dsm")
Warning: The working directory was changed to /Users/nkwiflores/Library/CloudStorage/OneDrive-UniversityofArizona/Documents/sandbox/Reproducibility/remote_sensing_fp/3_dsm_ortho/1_dsm inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.
# Bring point cloud to the environmental
exp_viron_point = readLAS("Final_Project_dsm_100cm.las")
Warning: Invalid data: 142596 points with a return number equal to 0 found.Warning: Invalid data: 142596 points with a number of returns equal to 0 found.
# plot the point cloud for the number of point.
plot(exp_viron_point, color = "RGB")
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
Error in dyn.load(dynlib <- getDynlib(dir)) :
unable to load shared object '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rgl/libs/rgl.so':
dlopen(/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rgl/libs/rgl.so, 0x0006): Library not loaded: /opt/X11/lib/libGLU.1.dylib
Referenced from: <C3FB2FAF-8517-3EBF-A366-EAF5C1C74ED7> /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rgl/libs/rgl.so
Reason: tried: '/opt/X11/lib/libGLU.1.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/opt/X11/lib/libGLU.1.dylib' (no such file), '/opt/X11/lib/libGLU.1.dylib' (no such file), '/Library/Frameworks/R.framework/Resources/lib/libGLU.1.dylib' (no such file), '/Library/Java/JavaVirtualMachines/jdk-11.0.18+10/Contents/Home/lib/server/libGLU.1.dylib' (no such file), '/var/folders/1c/1g5v5wx51hsghdscvwj0ymk80000gn/T/rstudio-fallback-library-path-1588176732/libGLU.1.dylib' (no such file)
Warning: Loading rgl's DLL failed.
This build of rgl depends on XQuartz, which failed to load.
See the discussion in https://stackoverflow.com/a/66127391/2554330Warning: Trying without OpenGL...Error: .onLoad failed in loadNamespace() for 'rgl', details:
call: rgl.init(initValue, onlyNULL)
error: OpenGL is not available in this build
Ground Filter It uses the lidT package to separate ground from
vegetation points called cloth simulation filter
(csf) Ref. Zhang et al. 2016
Ground points and canopy points rendered in RGB
Point clouds are different than the previous when the class_threshold
parameter was set to 0.1.
p_expl_experiment
[1] 525419.1 3570857.2
Creating Digital Terrain Model

p_explearning_dtm <- plot(explearning_dtm, xlab = "UTM_X", ylab = "UTM_Y", main = "Expt. Learning DTM")
mtext("Elevation (meter)", side = 4)
p_explearning_dtm
NULL
Export files out of R.
# export bigelow_DTM.tif file format, which is a Geotif format
writeRaster(explearning_dtm, filename = "explearning_dtm.tif", format = "GTiff", datatype = 'FLT4S', overwrite = TRUE)
Calculating vegetation height above ground.
Find the difference between canopy point clound and the digital
terrain model (dtm) to define the vegetation height.
#Normalization plot of the heights
explearning_agl = normalize_height(explearning_canopy_points, explearning_dtm, copy = TRUE)
Warning: 2036 points do not belong in the raster. Nearest neighbor was used to assign a value.
p_explearning_agl <- plot(explearning_agl)
# Range of height values uncleaned
range(explearning_agl$Z)
[1] -5.493 36.192
#Clean points for mistakes created by the photogrammetry program
explearning_agl_clean = filter_poi(explearning_agl, Z>0)
#range of height values cleaned
range(explearning_agl_clean$Z)
[1] 0.001 36.192
#plot the cleaned set
plot(explearning_agl_clean)
Voxelize Point Cloud The voxel is a 3D volume pixel for establising a
3D grid in the point cloud space.
compare canopy voxel (device 2) and explearning canoppy
points (device 1)
length(explearning_canopy_voxel)
[1] 1
number of volexs in the voxelized point cloud
print(explearning_num_voxels)
[1] 81781
total vegetation volume
print(total_veg_vol)
[1] 81.781
ID Individual Trees
#Segmented trees individually
treeseg4 = segment_trees(explearning_canopy_voxel, li2012(dt1 = 1.5, dt2 = 2, R= 2, Zu = 10, hmin = 1.5, speed_up = 5), attribute = "treeID")
p_treeseg4 <- plot(treeseg4, color = "treeID")
Unique number of trees there are 95 unique trees
# the number of unique individual trees
unique_trees = unique(treeseg4$treeID)
length(unique_trees)
[1] 4576
LS0tCnRpdGxlOiAiRXhwLiBMZWFybmluZyBFbnZpcm9ubWVudCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCmBgYHtyIHBhY2thZ2VzIGFuZCBsaWJyYXJpZXN9Cmluc3RhbGwucGFja2FnZXMoInJhc3RlciIpCmluc3RhbGwucGFja2FnZXMoInRlcnJhIikKaW5zdGFsbC5wYWNrYWdlcygibGlkUiIpCmluc3RhbGwucGFja2FnZXMoIlJDU0YiKQppbnN0YWxsLnBhY2thZ2VzKCJzcCIpCmluc3RhbGwucGFja2FnZXMoInJnbCIsIGRlcGVuZGVuY2llcyA9IFRSVUUpCmluc3RhbGwucGFja2FnZXMoInhxdWFydHoiKQpsaWJyYXJ5KHNwKQpsaWJyYXJ5KGxpZFIpCmxpYnJhcnkocmFzdGVyKQpsaWJyYXJ5KHRlcnJhKQpsaWJyYXJ5KFJDU0YpCmBgYAoKYGBge3IgbG9hZGluZyBmaWxlc30KIyB0aGlzIGlzIG15IHdvcmtwbGFjZQpzZXR3ZCgifi9MaWJyYXJ5L0Nsb3VkU3RvcmFnZS9PbmVEcml2ZS1Vbml2ZXJzaXR5b2ZBcml6b25hL0RvY3VtZW50cy9zYW5kYm94L1JlcHJvZHVjaWJpbGl0eS9yZW1vdGVfc2Vuc2luZ19mcC8zX2RzbV9vcnRoby8xX2RzbSIpCgojIEJyaW5nIHBvaW50IGNsb3VkIHRvIHRoZSBlbnZpcm9ubWVudGFsCmV4cF92aXJvbl9wb2ludCA9IHJlYWRMQVMoIkZpbmFsX1Byb2plY3RfZHNtXzEwMGNtLmxhcyIpCgojIHBsb3QgdGhlIHBvaW50IGNsb3VkIGZvciB0aGUgbnVtYmVyIG9mIHBvaW50LgogcGxvdChleHBfdmlyb25fcG9pbnQsIGNvbG9yID0gIlJHQiIpCmBgYApHcm91bmQgRmlsdGVyCkl0IHVzZXMgdGhlIGxpZFQgcGFja2FnZSB0byBzZXBhcmF0ZSBncm91bmQgZnJvbSB2ZWdldGF0aW9uIHBvaW50cyBjYWxsZWQgKioqY2xvdGggc2ltdWxhdGlvbiBmaWx0ZXIgKGNzZikqKiogUmVmLiBaaGFuZyBldCBhbC4gMjAxNgoKYGBge3IgR3JvdW5kIEZpbHRlcn0KIyBDbGFzc2lmeSBncm91bmQgYW5kIG5vbi1ncm91bmQgcG9pbnRzIHVzaW5nIHRoZSBjbG90aCBzaW11bGF0aW9uCmV4cF9lbnZfY2xhc3NpZmllZCA9IGNsYXNzaWZ5X2dyb3VuZChleHBfdmlyb25fcG9pbnQsIGFsZ29yaXRobSA9IGNzZihzbG9vcF9zbW9vdGggPSBGQUxTRSwgY2xhc3NfdGhyZXNob2xkID0gMC4xLCBjbG90aF9yZXNvbHV0aW9uID0gMC41LCByaWdpZG5lc3MgPSAyKSkKCiMgUGxvdCB0aGUgY2xhc3NpZnkgZ3JvdW5kIGFuZCBub24tZ3JvdW5kIHBvaW50cwpwbG90KGV4cF9lbnZfY2xhc3NpZmllZCwgY29sb3IgPSAiQ2xhc3NpZmljYXRpb24iKQpgYGAKCgpHcm91bmQgcG9pbnRzIGFuZCBjYW5vcHkgcG9pbnRzIHJlbmRlcmVkIGluIFJHQgpgYGB7ciBSR0IgcmVuZGVyfQojIGNyZWF0ZXMgZXhwZXJpbWVudGFsIGxlYXJuaW5nIGVudmlyb25tZW50IGdyb3VuZCBwb2ludHMKCiMgU2VwYXJhdGUgdGhlIHBvaW50IGNsb3VkIGludG86CiMgMS4gZ3JvdW5kIHBvaW50cyAoY2xhc3NpZmljYXRpb24gPT0gMikKZXhwX2Vudl9ncm91bmRfcG9pbnRzID0gZmlsdGVyX3BvaShleHBfZW52X2NsYXNzaWZpZWQsIENsYXNzaWZpY2F0aW9uID09IDIpCnBsb3RfZ3JvdW5kX3BvaW50IDwtIHBsb3QoZXhwX2Vudl9ncm91bmRfcG9pbnRzLCBjb2xvcj0iUkdCIikKcGxvdF9ncm91bmRfcG9pbnQKCiMyLiB2ZWdhdGF0aW9uIHBvaW50cyAoY2xhc3NpZmljYXRpb24gPT0gMCBvciAxKQpleHBsZWFybmluZ19jYW5vcHlfcG9pbnRzID0gZmlsdGVyX3BvaShleHBfZW52X2NsYXNzaWZpZWQsIENsYXNzaWZpY2F0aW9uID09IDApCnBfY2Fub3B5X3BvaW50cyA8LSBwbG90KGV4cGxlYXJuaW5nX2Nhbm9weV9wb2ludHMsIGNvbG9yID0gIlJHQiIpCgpgYGAKClBvaW50IGNsb3VkcyBhcmUgZGlmZmVyZW50IHRoYW4gdGhlIHByZXZpb3VzIHdoZW4gdGhlIGNsYXNzX3RocmVzaG9sZCBwYXJhbWV0ZXIgd2FzIHNldCB0byAwLjEuCgpgYGB7ciBjbGFzc190cmVzaG9sZH0KI3N0ZXAgMjAsIGNoYWdpbmcgdGhlIHRocmVzaG9sZApleHBsZWFybmluZ19jbGFzc2lmaWVkX2V4cGVyaW1lbnQgPSBjbGFzc2lmeV9ncm91bmQoZXhwX3Zpcm9uX3BvaW50LCBhbGdvcml0aG0gPSBjc2Yoc2xvb3Bfc21vb3RoID0gRkFMU0UsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBjbGFzc190aHJlc2hvbGQgPSAxLCBjbG90aF9yZXNvbHV0aW9uID0gMC41LCByaWdpZG5lc3MgPSAyKSkKCnBfZXhwbF9leHBlcmltZW50IDwtIHBsb3QoZXhwbGVhcm5pbmdfY2xhc3NpZmllZF9leHBlcmltZW50LCBjb2xvciA9ICJDbGFzc2lmaWNhdGlvbiIpCnBuZyhmaWxlID0gInBfZXhwbGVhcm5pbmdfZXhwZXJpbWVudC5wbmciLCB3aWR0aCA9IDYwMCwgaGVpZ2h0ID0gNDUwKQpwX2V4cGxfZXhwZXJpbWVudApkZXYub2ZmKCkKYGBgCgoKQ3JlYXRpbmcgRGlnaXRhbCBUZXJyYWluIE1vZGVsCgpgYGB7cn0KIyBTdGVwIDIyLiBVc2UgdGhlIGJpZ2Vsb3cgZ3JvdW5kIHBvaW50cyBmcm9tIHN0ZXAgMTggdG8gZ2VuZXJhdGUgZ3JpZGRlZCBkaWdpdGFsIG1vZGVsIChEVE0pCgojIFNwYXRpYWwgaW50ZXJwb2xhdGlvbjogay1uZWFyZXN0IG5laWdoYm9yIChLTk4pCiMgaW52ZXJzZS1kaXN0YW5jZSB3ZWlnaHRpbmcgKElEVykKCmV4cGxlYXJuaW5nX2R0bSA9IGdyaWRfdGVycmFpbihleHBfZW52X2dyb3VuZF9wb2ludHMsIHJlcyA9IDAuMSwgYWxnb3JpdGhtID0ga25uaWR3KGsgPSAxMCwgcCA9IDIpKQoKI1Bsb3QgdGhlIGdyaWRkZWQgZGlnaXRhbCBtb2RlbCAoZHRtKQpwX2V4cGxlYXJuaW5nX2R0bSA8LSBwbG90KGV4cGxlYXJuaW5nX2R0bSwgeGxhYiA9ICJVVE1fWCIsIHlsYWIgPSAiVVRNX1kiLCBtYWluID0gIkV4cHQuIExlYXJuaW5nIERUTSIpCm10ZXh0KCJFbGV2YXRpb24gKG1ldGVyKSIsIHNpZGUgPSA0KQpwX2V4cGxlYXJuaW5nX2R0bQoKYGBgCgpFeHBvcnQgZmlsZXMgb3V0IG9mIFIuCgpgYGB7cn0KIyBleHBvcnQgYmlnZWxvd19EVE0udGlmIGZpbGUgZm9ybWF0LCB3aGljaCBpcyBhIEdlb3RpZiBmb3JtYXQKd3JpdGVSYXN0ZXIoZXhwbGVhcm5pbmdfZHRtLCBmaWxlbmFtZSA9ICJleHBsZWFybmluZ19kdG0udGlmIiwgZm9ybWF0ID0gIkdUaWZmIiwgZGF0YXR5cGUgPSAnRkxUNFMnLCBvdmVyd3JpdGUgPSBUUlVFKQoKYGBgCgpDYWxjdWxhdGluZyB2ZWdldGF0aW9uIGhlaWdodCBhYm92ZSBncm91bmQuCi0tLQpGaW5kIHRoZSBkaWZmZXJlbmNlIGJldHdlZW4gY2Fub3B5IHBvaW50IGNsb3VuZCBhbmQgdGhlIGRpZ2l0YWwgdGVycmFpbiBtb2RlbCAoZHRtKSB0byBkZWZpbmUgdGhlIHZlZ2V0YXRpb24gaGVpZ2h0LgotLS0KYGBge3J9CiNOb3JtYWxpemF0aW9uIHBsb3Qgb2YgdGhlIGhlaWdodHMKZXhwbGVhcm5pbmdfYWdsID0gbm9ybWFsaXplX2hlaWdodChleHBsZWFybmluZ19jYW5vcHlfcG9pbnRzLCBleHBsZWFybmluZ19kdG0sIGNvcHkgPSBUUlVFKQpwX2V4cGxlYXJuaW5nX2FnbCA8LSBwbG90KGV4cGxlYXJuaW5nX2FnbCkKCiMgUmFuZ2Ugb2YgaGVpZ2h0IHZhbHVlcyB1bmNsZWFuZWQKcmFuZ2UoZXhwbGVhcm5pbmdfYWdsJFopCgojQ2xlYW4gcG9pbnRzIGZvciBtaXN0YWtlcyBjcmVhdGVkIGJ5IHRoZSBwaG90b2dyYW1tZXRyeSBwcm9ncmFtCmV4cGxlYXJuaW5nX2FnbF9jbGVhbiA9IGZpbHRlcl9wb2koZXhwbGVhcm5pbmdfYWdsLCBaPjApCiNyYW5nZSBvZiBoZWlnaHQgdmFsdWVzIGNsZWFuZWQgCnJhbmdlKGV4cGxlYXJuaW5nX2FnbF9jbGVhbiRaKQojcGxvdCB0aGUgY2xlYW5lZCBzZXQKcGxvdChleHBsZWFybmluZ19hZ2xfY2xlYW4pCmBgYAoKClZveGVsaXplIFBvaW50IENsb3VkClRoZSB2b3hlbCBpcyBhIDNEIHZvbHVtZSBwaXhlbCBmb3IgZXN0YWJsaXNpbmcgYSAzRCBncmlkIGluIHRoZSBwb2ludCBjbG91ZCBzcGFjZS4KCioqKmNvbXBhcmUgY2Fub3B5IHZveGVsIChkZXZpY2UgMikgYW5kIGV4cGxlYXJuaW5nIGNhbm9wcHkgcG9pbnRzIChkZXZpY2UgMSkgKioqCmBgYHtyfQpleHBsZWFybmluZ19jYW5vcHlfdm94ZWwgPSB2b3hlbGl6ZV9wb2ludHMoZXhwbGVhcm5pbmdfYWdsX2NsZWFuLCByZXMgPSAwLjEpCmxlbmd0aChleHBsZWFybmluZ19jYW5vcHlfdm94ZWwpCnBfY2Fub3B5X3ZveGVsIDwtIHBsb3QoZXhwbGVhcm5pbmdfY2Fub3B5X3ZveGVsKQpgYGAKCm51bWJlciBvZiB2b2xleHMgaW4gdGhlIHZveGVsaXplZCBwb2ludCBjbG91ZApgYGB7cn0KIyB2b3hlbCBwb2ludCBudW1iZXIgdXNpbmcgbWF0cml4CmV4cGxlYXJuaW5nX2Nhbm9weV92b3hlbF9udW0gPC0gdm94ZWxfbWV0cmljcyhleHBsZWFybmluZ19jYW5vcHlfcG9pbnRzLCB+bGVuZ3RoKFopLCByZXMgPSAwLjEpCiNjb3VudCB2b3hlbHMgYnkgY29udmVydGluZyB0byBkYXRhIGZyYW0gdG8gY291bnQgZWFyY2ggcm93IHJlcHJlc2VudGluZyB0aGUgdm94ZWwKZXhwbGVhcm5pbmdfbnVtX3ZveGVscyA9IG5yb3coYXMuZGF0YS5mcmFtZShleHBsZWFybmluZ19jYW5vcHlfdm94ZWxfbnVtKSkKcHJpbnQoZXhwbGVhcm5pbmdfbnVtX3ZveGVscykKYGBgCgoKdG90YWwgdmVnZXRhdGlvbiB2b2x1bWUgCmBgYHtyfQojIHZvbHVtZSBvZiBzaW5nbGUgdm9sZXgsIHdoaWNoIGlzIDEwY20geCAxMGNtIHggMTBjbSA9IGN1YmUKIyB2b2x1bWUgcGYgdm94ZWwgPSAoMC4xbSleMyA9IDAuMDAxIG1eMwp2b3hlbF92b2wgPSAwLjAwMQp0b3RhbF92ZWdfdm9sID0gZXhwbGVhcm5pbmdfbnVtX3ZveGVscyAqIHZveGVsX3ZvbApwcmludCh0b3RhbF92ZWdfdm9sKQpgYGAKCklEIEluZGl2aWR1YWwgVHJlZXMKCmBgYHtyfQojU2VnbWVudGVkIHRyZWVzIGluZGl2aWR1YWxseQp0cmVlc2VnNCA9IHNlZ21lbnRfdHJlZXMoZXhwbGVhcm5pbmdfY2Fub3B5X3ZveGVsLCBsaTIwMTIoZHQxID0gMS41LCBkdDIgPSAyLCBSPSAyLCBadSA9IDEwLCBobWluID0gMS41LCBzcGVlZF91cCA9IDUpLCBhdHRyaWJ1dGUgPSAidHJlZUlEIikKcF90cmVlc2VnNCA8LSBwbG90KHRyZWVzZWc0LCBjb2xvciA9ICJ0cmVlSUQiKQpgYGAKCgpVbmlxdWUgbnVtYmVyIG9mIHRyZWVzCnRoZXJlIGFyZSA5NSB1bmlxdWUgdHJlZXMgCmBgYHtyfQojIHRoZSBudW1iZXIgb2YgdW5pcXVlIGluZGl2aWR1YWwgdHJlZXMKdW5pcXVlX3RyZWVzID0gdW5pcXVlKHRyZWVzZWc0JHRyZWVJRCkKbGVuZ3RoKHVuaXF1ZV90cmVlcykKCmBgYAoK