2014 Tile Cleaning

# Load 2014 tile
las_2014 <- readLAS("18TWK970955.las")
# Problems to fix: 

# EPSG code 32767 unknown
# There are 1,111 duplicated points
# Remove duplicate points
las_2014 <- filter_duplicates(las_2014)
# Assign CRS to header, which can be found using summary() function
st_crs(las_2014) <- 6347
# Change meters to feet for Z factor, for later analysis
las_2014@data$Z <- las_2014@data$Z * 3.28084
# Ensure point cloud is valid after changing Z factor
las_quantize(las_2014)
# Update object and header
las_2014 <- las_update(las_2014)

2017 Tile Cleaning

# Load 2017 tile
las_2017 <- readLAS("27160.las")
# Problems to fix:

# 2,841 points are duplicated and share XYZ coordinates with other points
# 4,397,012 points flagged 'withheld'
# Remove duplicate points
las_2017 <- filter_duplicates(las_2017)
# Remove withheld points (based on classmate Alex Jacobson's solution)
las_2017 <- filter_poi(las_2017, !Withheld_flag)

Last Returns

# Filter ground classification (2)
ground_2014 = filter_poi(las_2014, Classification == 2)

ground_2017 = filter_poi(las_2017, Classification == 2)
plot(ground_2014)
plot(ground_2017)

TIN Interpolation

# Reproject to 2017's CRS
tin_2014_aligned <- project(tin_2014, tin_2017)
# 2014 DEM
wgs_2014 <- project(tin_2014_aligned, "EPSG:4326")
terra::plot(wgs_2014,
            main = "Yellowbar DEM, 2014",
            plg = list(title = "Elevation (ft)")
)

# 2017 DEM
wgs_2017 <- project(tin_2017, "EPSG:4326")
terra::plot(wgs_2017,
            main = "Yellowbar DEM, 2017",
            plg = list(title = "Elevation (ft)")
)

Calculating Elevation Change!

elevation_change <- tin_2017 - tin_2014_aligned
wgs_elev <- project(elevation_change, "EPSG:4326")
terra::plot(wgs_elev,
            main = "Yellowbar Elevation Change",
            plg = list(title = "Feet")
)

Wang et al. (2023) Organic Matter Equation

# x = organic matter
# y = accretion rate

# y = 0.15x + 0.12

# So to solve for organic matter inventory,

# x = (y - 0.12)/0.15

organic <- (elevation_change - 0.12) / 0.15
wgs_organic <- project(organic, "EPSG:4326")
terra::plot(wgs_organic,
            main = "Organic Matter Gain/Loss",
            plg = list(title = "Feet")
)