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