library(lidR)

library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
## 
## Attaching package: 'sf'
## The following object is masked from 'package:lidR':
## 
##     st_concave_hull
library(sp)
## 
## Attaching package: 'sp'
## The following object is masked from 'package:lidR':
## 
##     wkt
library(raster)
## 
## Attaching package: 'raster'
## The following objects are masked from 'package:lidR':
## 
##     projection, projection<-
library(terra)
## terra 1.7.71
## 
## Attaching package: 'terra'
## The following objects are masked from 'package:lidR':
## 
##     area, crs, crs<-, is.empty
library(ggplot2)

library(mapview)

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract(), raster::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ dplyr::select()  masks raster::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Read in tiles

las.data.2014 <- readLAS("/Users/jamesreo/Desktop/Blue Carbon Data/18TWK970955.las")

las.data.2017 <- readLAS("/Users/jamesreo/Desktop/Blue Carbon Data/25162.las")
## Warning: There are 6661036 points flagged 'withheld'.

Filter

filter.las.data.2014 <- filter_last(las.data.2014)

filter.las.data.2017 <- filter_last(las.data.2017)

#clean house to free storage 
rm(las.data.2014)

rm(las.data.2017)
#Original data for 2017 filtered data is too voluminus, classifying as per Andy's method allows for dsm creation

classified.2017 <- classify_poi(filter.las.data.2017, class = as.integer(9), poi = ~Classification %in% c(40, 41, 45))

Rasterize

dsm2014 <- rasterize_terrain(filter.las.data.2014, res=1, algorithm=tin(), pkg= 'raster')
## Delaunay rasterization[================----------------------------------] 32% (1 threads)
Delaunay rasterization[================----------------------------------] 33% (1 threads)
Delaunay rasterization[=================---------------------------------] 34% (1 threads)
Delaunay rasterization[=================---------------------------------] 35% (1 threads)
Delaunay rasterization[==================--------------------------------] 36% (1 threads)
Delaunay rasterization[==================--------------------------------] 37% (1 threads)
Delaunay rasterization[===================-------------------------------] 38% (1 threads)
Delaunay rasterization[===================-------------------------------] 39% (1 threads)
Delaunay rasterization[====================------------------------------] 40% (1 threads)
Delaunay rasterization[====================------------------------------] 41% (1 threads)
Delaunay rasterization[=====================-----------------------------] 42% (1 threads)
Delaunay rasterization[=====================-----------------------------] 43% (1 threads)
Delaunay rasterization[======================----------------------------] 44% (1 threads)
Delaunay rasterization[======================----------------------------] 45% (1 threads)
Delaunay rasterization[=======================---------------------------] 46% (1 threads)
Delaunay rasterization[=======================---------------------------] 47% (1 threads)
Delaunay rasterization[========================--------------------------] 48% (1 threads)
Delaunay rasterization[========================--------------------------] 49% (1 threads)
Delaunay rasterization[=========================-------------------------] 50% (1 threads)
Delaunay rasterization[=========================-------------------------] 51% (1 threads)
Delaunay rasterization[==========================------------------------] 52% (1 threads)
Delaunay rasterization[==========================------------------------] 53% (1 threads)
Delaunay rasterization[===========================-----------------------] 54% (1 threads)
Delaunay rasterization[===========================-----------------------] 55% (1 threads)
Delaunay rasterization[============================----------------------] 56% (1 threads)
Delaunay rasterization[============================----------------------] 57% (1 threads)
Delaunay rasterization[=============================---------------------] 58% (1 threads)
Delaunay rasterization[=============================---------------------] 59% (1 threads)
Delaunay rasterization[==============================--------------------] 60% (1 threads)
Delaunay rasterization[==============================--------------------] 61% (1 threads)
Delaunay rasterization[===============================-------------------] 62% (1 threads)
Delaunay rasterization[===============================-------------------] 63% (1 threads)
Delaunay rasterization[================================------------------] 64% (1 threads)
Delaunay rasterization[================================------------------] 65% (1 threads)
Delaunay rasterization[=================================-----------------] 66% (1 threads)
Delaunay rasterization[=================================-----------------] 67% (1 threads)
Delaunay rasterization[==================================----------------] 68% (1 threads)
Delaunay rasterization[==================================----------------] 69% (1 threads)
Delaunay rasterization[===================================---------------] 70% (1 threads)
Delaunay rasterization[===================================---------------] 71% (1 threads)
Delaunay rasterization[====================================--------------] 72% (1 threads)
Delaunay rasterization[====================================--------------] 73% (1 threads)
Delaunay rasterization[=====================================-------------] 74% (1 threads)
Delaunay rasterization[=====================================-------------] 75% (1 threads)
Delaunay rasterization[======================================------------] 76% (1 threads)
Delaunay rasterization[======================================------------] 77% (1 threads)
Delaunay rasterization[=======================================-----------] 78% (1 threads)
Delaunay rasterization[=======================================-----------] 79% (1 threads)
Delaunay rasterization[========================================----------] 80% (1 threads)
Delaunay rasterization[========================================----------] 81% (1 threads)
Delaunay rasterization[=========================================---------] 82% (1 threads)
Delaunay rasterization[=========================================---------] 83% (1 threads)
Delaunay rasterization[==========================================--------] 84% (1 threads)
Delaunay rasterization[==========================================--------] 85% (1 threads)
Delaunay rasterization[===========================================-------] 86% (1 threads)
Delaunay rasterization[===========================================-------] 87% (1 threads)
Delaunay rasterization[============================================------] 88% (1 threads)
Delaunay rasterization[============================================------] 89% (1 threads)
Delaunay rasterization[=============================================-----] 90% (1 threads)
Delaunay rasterization[=============================================-----] 91% (1 threads)
Delaunay rasterization[==============================================----] 92% (1 threads)
Delaunay rasterization[==============================================----] 93% (1 threads)
Delaunay rasterization[===============================================---] 94% (1 threads)
Delaunay rasterization[===============================================---] 95% (1 threads)
Delaunay rasterization[================================================--] 96% (1 threads)
Delaunay rasterization[================================================--] 97% (1 threads)
Delaunay rasterization[=================================================-] 98% (1 threads)
Delaunay rasterization[=================================================-] 99% (1 threads)
Delaunay rasterization[==================================================] 100% (1 threads)
dsm2017 <- rasterize_terrain(classified.2017, res=1, algorithm=tin(), pkg= 'raster')
## Delaunay rasterization[=-------------------------------------------------] 3% (1 threads)
Delaunay rasterization[==------------------------------------------------] 4% (1 threads)
Delaunay rasterization[==------------------------------------------------] 5% (1 threads)
Delaunay rasterization[===-----------------------------------------------] 6% (1 threads)
Delaunay rasterization[===-----------------------------------------------] 7% (1 threads)
Delaunay rasterization[====----------------------------------------------] 8% (1 threads)
Delaunay rasterization[====----------------------------------------------] 9% (1 threads)
Delaunay rasterization[=====---------------------------------------------] 10% (1 threads)
Delaunay rasterization[=====---------------------------------------------] 11% (1 threads)
Delaunay rasterization[======--------------------------------------------] 12% (1 threads)
Delaunay rasterization[======--------------------------------------------] 13% (1 threads)
Delaunay rasterization[=======-------------------------------------------] 14% (1 threads)
Delaunay rasterization[=======-------------------------------------------] 15% (1 threads)
Delaunay rasterization[========------------------------------------------] 16% (1 threads)
Delaunay rasterization[========------------------------------------------] 17% (1 threads)
Delaunay rasterization[=========-----------------------------------------] 18% (1 threads)
Delaunay rasterization[=========-----------------------------------------] 19% (1 threads)
Delaunay rasterization[==========----------------------------------------] 20% (1 threads)
Delaunay rasterization[==========----------------------------------------] 21% (1 threads)
Delaunay rasterization[===========---------------------------------------] 22% (1 threads)
Delaunay rasterization[===========---------------------------------------] 23% (1 threads)
Delaunay rasterization[============--------------------------------------] 24% (1 threads)
Delaunay rasterization[============--------------------------------------] 25% (1 threads)
Delaunay rasterization[=============-------------------------------------] 26% (1 threads)
Delaunay rasterization[=============-------------------------------------] 27% (1 threads)
Delaunay rasterization[==============------------------------------------] 28% (1 threads)
Delaunay rasterization[==============------------------------------------] 29% (1 threads)
Delaunay rasterization[===============-----------------------------------] 30% (1 threads)
Delaunay rasterization[===============-----------------------------------] 31% (1 threads)
Delaunay rasterization[================----------------------------------] 32% (1 threads)
Delaunay rasterization[================----------------------------------] 33% (1 threads)
Delaunay rasterization[=================---------------------------------] 34% (1 threads)
Delaunay rasterization[=================---------------------------------] 35% (1 threads)
Delaunay rasterization[==================--------------------------------] 36% (1 threads)
Delaunay rasterization[==================--------------------------------] 37% (1 threads)
Delaunay rasterization[===================-------------------------------] 38% (1 threads)
Delaunay rasterization[===================-------------------------------] 39% (1 threads)
Delaunay rasterization[====================------------------------------] 40% (1 threads)
Delaunay rasterization[====================------------------------------] 41% (1 threads)
Delaunay rasterization[=====================-----------------------------] 42% (1 threads)
Delaunay rasterization[=====================-----------------------------] 43% (1 threads)
Delaunay rasterization[======================----------------------------] 44% (1 threads)
Delaunay rasterization[======================----------------------------] 45% (1 threads)
Delaunay rasterization[=======================---------------------------] 46% (1 threads)
Delaunay rasterization[=======================---------------------------] 47% (1 threads)
Delaunay rasterization[========================--------------------------] 48% (1 threads)
Delaunay rasterization[========================--------------------------] 49% (1 threads)
Delaunay rasterization[=========================-------------------------] 50% (1 threads)
Delaunay rasterization[=========================-------------------------] 51% (1 threads)
Delaunay rasterization[==========================------------------------] 52% (1 threads)
Delaunay rasterization[==========================------------------------] 53% (1 threads)
Delaunay rasterization[===========================-----------------------] 54% (1 threads)
Delaunay rasterization[===========================-----------------------] 55% (1 threads)
Delaunay rasterization[============================----------------------] 56% (1 threads)
Delaunay rasterization[============================----------------------] 57% (1 threads)
Delaunay rasterization[=============================---------------------] 58% (1 threads)
Delaunay rasterization[=============================---------------------] 59% (1 threads)
Delaunay rasterization[==============================--------------------] 60% (1 threads)
Delaunay rasterization[==============================--------------------] 61% (1 threads)
Delaunay rasterization[===============================-------------------] 62% (1 threads)
Delaunay rasterization[===============================-------------------] 63% (1 threads)
Delaunay rasterization[================================------------------] 64% (1 threads)
Delaunay rasterization[================================------------------] 65% (1 threads)
Delaunay rasterization[=================================-----------------] 66% (1 threads)
Delaunay rasterization[=================================-----------------] 67% (1 threads)
Delaunay rasterization[==================================----------------] 68% (1 threads)
Delaunay rasterization[==================================----------------] 69% (1 threads)
Delaunay rasterization[===================================---------------] 70% (1 threads)
Delaunay rasterization[===================================---------------] 71% (1 threads)
Delaunay rasterization[====================================--------------] 72% (1 threads)
Delaunay rasterization[====================================--------------] 73% (1 threads)
Delaunay rasterization[=====================================-------------] 74% (1 threads)
Delaunay rasterization[=====================================-------------] 75% (1 threads)
Delaunay rasterization[======================================------------] 76% (1 threads)
Delaunay rasterization[======================================------------] 77% (1 threads)
Delaunay rasterization[=======================================-----------] 78% (1 threads)
Delaunay rasterization[=======================================-----------] 79% (1 threads)
Delaunay rasterization[========================================----------] 80% (1 threads)
Delaunay rasterization[========================================----------] 81% (1 threads)
Delaunay rasterization[=========================================---------] 82% (1 threads)
Delaunay rasterization[=========================================---------] 83% (1 threads)
Delaunay rasterization[==========================================--------] 84% (1 threads)
Delaunay rasterization[==========================================--------] 85% (1 threads)
Delaunay rasterization[===========================================-------] 86% (1 threads)
Delaunay rasterization[===========================================-------] 87% (1 threads)
Delaunay rasterization[============================================------] 88% (1 threads)
Delaunay rasterization[============================================------] 89% (1 threads)
Delaunay rasterization[=============================================-----] 90% (1 threads)
Delaunay rasterization[=============================================-----] 91% (1 threads)
Delaunay rasterization[==============================================----] 92% (1 threads)
Delaunay rasterization[==============================================----] 93% (1 threads)
Delaunay rasterization[===============================================---] 94% (1 threads)
Delaunay rasterization[===============================================---] 95% (1 threads)
Delaunay rasterization[================================================--] 96% (1 threads)
Delaunay rasterization[================================================--] 97% (1 threads)
Delaunay rasterization[=================================================-] 98% (1 threads)
Delaunay rasterization[=================================================-] 99% (1 threads)
Delaunay rasterization[==================================================] 100% (1 threads)
crs(dsm2014) <- CRS("EPSG:6347")

dsm2014.reproject <- projectRaster(dsm2014, crs = crs(dsm2017))

Crop, resample 2014 dsm to match 2017 grid cells,correct unit, and create difference surface.

dsm2014.cropped<- crop(dsm2014.reproject,dsm2017)

dsm2014.resample <- resample(dsm2014.cropped,dsm2017,method="bilinear")

#Correct unit for 2014 to be in feet 
dsm2014.ft <- dsm2014.resample * 3.28

dsm.difference.ft <- dsm2017- dsm2014.ft


#free up storage 
rm(classified.2017)
rm(filter.las.data.2014)
rm(filter.las.data.2017)

#Write Raster
writeRaster(dsm2014.ft, "dsm2014.ft.tif", overwrite=TRUE)
writeRaster(dsm2017,'dsm2017.tif',overwrite=TRUE)
writeRaster(dsm.difference.ft,'difference.ft.tif',overwrite=TRUE)

Plot 2014

plot(dsm2014.ft,main='Elevation 2014(ft)')

Plot 2017

plot(dsm2017,main='Elevation 2017(ft)')

plot(dsm.difference.ft,main='Elevation Change 2014-2017(ft)')

#change to cm 
 dsm.difference.cm <- (dsm.difference.ft / 3) * 30.48

hist(dsm.difference.cm,
     breaks = 6,
     main = 'Annual Accretion Change 2014-2017',
     col ='green',
     xlab='Elevation Change(cm)')

organic.matter <- 6.6667 * dsm.difference.cm - 0.8

mapview(organic.matter)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 6250000 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  6250000 '
hist(organic.matter,
     breaks = 6,
     main = 'Calculated Organic Matter Inventory(g cm^2)',
     col='blue', 
     xlab='OMI(g/cm^2)')