R Spatial Final Project, Alex Jacobson

Setup

knitr::opts_chunk$set(echo = TRUE)

# Making sure all needed packages are installed

# If running this project for the first time: install.packages(c("tidyverse", "magrittr", "janitor", "lidR", "terra", "sf", "tmap")

#Making sure all needed packages are loaded

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
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(magrittr);

## 
## Attaching package: 'magrittr'
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(janitor);

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(lidR);
library(terra);

## terra 1.8.42
## 
## Attaching package: 'terra'
## 
## The following objects are masked from 'package:lidR':
## 
##     area, crs, crs<-, is.empty, watershed
## 
## The following object is masked from 'package:janitor':
## 
##     crosstab
## 
## The following objects are masked from 'package:magrittr':
## 
##     extract, inset
## 
## The following object is masked from 'package:tidyr':
## 
##     extract

library(sf);

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
## 
## Attaching package: 'sf'
## 
## The following object is masked from 'package:lidR':
## 
##     st_concave_hull

library(tmap);
library(maptiles)

# Loading in LiDAR data using lidR readLAS function

# First, the four tiles for the 2017 data: 2017 data uses smaller tiles than the 2014 data, so it was necessary to have four 2017 tiles to approximate one 2014 tile

a2017 <- readLAS("~/Desktop/Data Viz with R/R Spatial Project Alex/2017/30165.las")

## Warning: There are 3702152 points flagged 'withheld'.

b2017 <- readLAS("~/Desktop/Data Viz with R/R Spatial Project Alex/2017/30167.las")

## Warning: There are 3377152 points flagged 'withheld'.

c2017 <- readLAS("~/Desktop/Data Viz with R/R Spatial Project Alex/2017/32165.las")

## Warning: There are 2028861 points flagged 'withheld'.

d2017 <- readLAS("~/Desktop/Data Viz with R/R Spatial Project Alex/2017/32167.las")

## Warning: There are 3185123 points flagged 'withheld'.

tile2014 <- readLAS("~/Desktop/Data Viz with R/R Spatial Project Alex//18TWK985970.las")

# First experimental plot using 2014 data, copied from lidR documentation to begin learning how to use plots

# Plot object by scan angle, 
# make the background white, 
# display XYZ axis and  scale colors
lidR::plot(tile2014, color = "ScanAngleRank", bg = "white", axis = TRUE, legend = TRUE)

Data Cleaning

# Validating data, following steps in lidR documentation
las_check(tile2014)

## 
##  Checking the data
##   - Checking coordinates... ✓
##   - Checking coordinates type... ✓
##   - Checking coordinates range... ✓
##   - Checking coordinates quantization... ✓
##   - Checking attributes type... ✓
##   - Checking ReturnNumber validity... ✓
##   - Checking NumberOfReturns validity... ✓
##   - Checking ReturnNumber vs. NumberOfReturns... ✓
##   - Checking RGB validity... ✓
##   - Checking absence of NAs... ✓
##   - Checking duplicated points...
##     ⚠ 2404 points are duplicated and share XYZ coordinates with other points
##   - Checking degenerated ground points...
##     ⚠ There were 1 degenerated ground points. Some X Y Z coordinates were repeated
##   - Checking attribute population... ✓
##   - Checking gpstime incoherances
##     ✗ 4864066 pulses (points with the same gpstime) have points with identical ReturnNumber
##   - Checking flag attributes... ✓
##   - Checking user data attribute... ✓
##  Checking the header
##   - Checking header completeness... ✓
##   - Checking scale factor validity... ✓
##   - Checking point data format ID validity... ✓
##   - Checking extra bytes attributes validity... ✓
##   - Checking the bounding box validity... ✓
##   - Checking coordinate reference system... ✓
##  Checking header vs data adequacy
##   - Checking attributes vs. point format... ✓
##   - Checking header bbox vs. actual content... ✓
##   - Checking header number of points vs. actual content... ✓
##   - Checking header return number vs. actual content... ✓
##  Checking coordinate reference system...
##   - Checking if the CRS was understood by R...
##     ✗ EPSG code 32767 unknown
##  Checking preprocessing already done 
##   - Checking ground classification... yes
##   - Checking normalization... no
##   - Checking negative outliers...
##     ⚠ 1927436 points below 0
##   - Checking flightline classification... yes
##  Checking compression
##   - Checking attribute compression...
##    -  Synthetic_flag is compressed
##    -  Keypoint_flag is compressed
##    -  Withheld_flag is compressed
##    -  UserData is compressed

# Cleaning up 2014 data based on las_check()

# Step 1: Fix the issue with the CRS being unknown

st_crs(tile2014) <- 6347


# Step 2: Getting rid of the 2404 duplicate points
tile2014 <- filter_duplicates(tile2014)

# Repeating las_check to see what's fixed and what's still a problem

las_check(tile2014)

## 
##  Checking the data
##   - Checking coordinates... ✓
##   - Checking coordinates type... ✓
##   - Checking coordinates range... ✓
##   - Checking coordinates quantization... ✓
##   - Checking attributes type... ✓
##   - Checking ReturnNumber validity... ✓
##   - Checking NumberOfReturns validity... ✓
##   - Checking ReturnNumber vs. NumberOfReturns... ✓
##   - Checking RGB validity... ✓
##   - Checking absence of NAs... ✓
##   - Checking duplicated points... ✓
##   - Checking degenerated ground points... ✓
##   - Checking attribute population... ✓
##   - Checking gpstime incoherances
##     ✗ 4862387 pulses (points with the same gpstime) have points with identical ReturnNumber
##   - Checking flag attributes... ✓
##   - Checking user data attribute... ✓
##  Checking the header
##   - Checking header completeness... ✓
##   - Checking scale factor validity... ✓
##   - Checking point data format ID validity... ✓
##   - Checking extra bytes attributes validity... ✓
##   - Checking the bounding box validity... ✓
##   - Checking coordinate reference system... ✓
##  Checking header vs data adequacy
##   - Checking attributes vs. point format... ✓
##   - Checking header bbox vs. actual content... ✓
##   - Checking header number of points vs. actual content... ✓
##   - Checking header return number vs. actual content... ✓
##  Checking coordinate reference system...
##   - Checking if the CRS was understood by R... ✓
##  Checking preprocessing already done 
##   - Checking ground classification... yes
##   - Checking normalization... no
##   - Checking negative outliers...
##     ⚠ 1927213 points below 0
##   - Checking flightline classification... yes
##  Checking compression
##   - Checking attribute compression...
##    -  Synthetic_flag is compressed
##    -  Keypoint_flag is compressed
##    -  Withheld_flag is compressed
##    -  UserData is compressed

print(tile2014)

## class        : LAS (v1.2 format 1)
## memory       : 782.9 Mb 
## extent       : 598500, 6e+05, 4497000, 4498500 (xmin, xmax, ymin, ymax)
## coord. ref.  : NAD83(2011) / UTM zone 18N 
## area         : 2.24 km²
## points       : 12.83 million points
## density      : 5.71 points/m²
## density      : 5.31 pulses/m²

print(tile2014)

## class        : LAS (v1.2 format 1)
## memory       : 782.9 Mb 
## extent       : 598500, 6e+05, 4497000, 4498500 (xmin, xmax, ymin, ymax)
## coord. ref.  : NAD83(2011) / UTM zone 18N 
## area         : 2.24 km²
## points       : 12.83 million points
## density      : 5.71 points/m²
## density      : 5.31 pulses/m²

# Negative elevation points are probably fine because we're looking at shoreline

# Normalization...no
# In the documentation, it mentions the goal of normalization is to get flat terrain, which we don't need cause we're focusing on ground points/bare earth


# Looking at classifications
unique(tile2014$Classification)

## [1]  1  9  7 10  2

# Filtering out all but last returns, as per project instructions 

tile2014 <- filter_poi(tile2014, Classification == 2)

print(tile2014)

## class        : LAS (v1.2 format 1)
## memory       : 145 Mb 
## extent       : 598500, 6e+05, 4497000, 4498500 (xmin, xmax, ymin, ymax)
## coord. ref.  : NAD83(2011) / UTM zone 18N 
## area         : 1.59 km²
## points       : 2.38 million points
## density      : 1.5 points/m²
## density      : 1.32 pulses/m²

# Merging 2017 tiles into one tile for easier processing later on

# Using LAScatalog processing engine for better results, based on what it recommends in the documentation 

#removing smaller datasets to free up space, since am going with this new approach

#rm(a2017)
#rm(b2017)
#rm(c2017)
#rm(d2017)

ctg2017 <- readLAScatalog("~/Desktop/Data Viz with R/R Spatial Project Alex/2017")

las2017 <- readLAS(ctg2017)

## Warning: There are 12293288 points flagged 'withheld'.

las_check(las2017)

## 
##  Checking the data
##   - Checking coordinates... ✓
##   - Checking coordinates type... ✓
##   - Checking coordinates range... ✓
##   - Checking coordinates quantization... ✓
##   - Checking attributes type... ✓
##   - Checking ReturnNumber validity... ✓
##   - Checking NumberOfReturns validity... ✓
##   - Checking ReturnNumber vs. NumberOfReturns... ✓
##   - Checking RGB validity... ✓
##   - Checking absence of NAs... ✓
##   - Checking duplicated points...
##     ⚠ 19267 points are duplicated and share XYZ coordinates with other points
##   - Checking degenerated ground points...
##     ⚠ There were 1 degenerated ground points. Some X Y Z coordinates were repeated
##     ⚠ There were 4 degenerated ground points. Some X Y coordinates were repeated but with different Z coordinates
##   - Checking attribute population... ✓
##   - Checking gpstime incoherances
##     ✗ 6336543 pulses (points with the same gpstime) have points with identical ReturnNumber
##   - Checking flag attributes...
##     🛈 12293288 points flagged 'withheld'
##   - Checking user data attribute...
##     🛈 36203865 points have a non 0 UserData attribute. This probably has a meaning
##  Checking the header
##   - Checking header completeness... ✓
##   - Checking scale factor validity... ✓
##   - Checking point data format ID validity... ✓
##   - Checking extra bytes attributes validity... ✓
##   - Checking the bounding box validity... ✓
##   - Checking coordinate reference system... ✓
##  Checking header vs data adequacy
##   - Checking attributes vs. point format... ✓
##   - Checking header bbox vs. actual content... ✓
##   - Checking header number of points vs. actual content... ✓
##   - Checking header return number vs. actual content... ✓
##  Checking coordinate reference system...
##   - Checking if the CRS was understood by R... ✓
##  Checking preprocessing already done 
##   - Checking ground classification... yes
##   - Checking normalization... no
##   - Checking negative outliers...
##     ⚠ 24165921 points below 0
##   - Checking flightline classification... yes
##  Checking compression
##   - Checking attribute compression...
##    -  Synthetic_flag is compressed
##    -  Keypoint_flag is compressed

# Checking CRS as advised

print(las2017)

## class        : LAS (v1.4 format 6)
## memory       : 3.6 Gb 
## extent       : 1030000, 1035000, 165000, 170000 (xmin, xmax, ymin, ymax)
## coord. ref.  : NAD83 / New York Long Island (ftUS) + NAVD88 height (ftUS) 
## area         : 21.62 kus-ft²
## points       : 46.29 million points
## density      : 2.14 points/us-ft²
## density      : 1.61 pulses/us-ft²

#Based on an online search, EPSG is EPSG:2263

las2017 <- filter_duplicates(las2017)

#Checking to see what the classification value for withheld points is, so I can remove them
unique(las2017$Classification)

## [1]  7  2  1 10  9 45 40 41

# Getting rid of withheld points - based on online resources, the code for withheld points should be 7

las2017 <- filter_poi(las2017, !Withheld_flag)

las_check(las2017)

## 
##  Checking the data
##   - Checking coordinates... ✓
##   - Checking coordinates type... ✓
##   - Checking coordinates range... ✓
##   - Checking coordinates quantization... ✓
##   - Checking attributes type... ✓
##   - Checking ReturnNumber validity... ✓
##   - Checking NumberOfReturns validity... ✓
##   - Checking ReturnNumber vs. NumberOfReturns... ✓
##   - Checking RGB validity... ✓
##   - Checking absence of NAs... ✓
##   - Checking duplicated points... ✓
##   - Checking degenerated ground points...
##     ⚠ There were 4 degenerated ground points. Some X Y coordinates were repeated but with different Z coordinates
##   - Checking attribute population... ✓
##   - Checking gpstime incoherances
##     ✗ 4290368 pulses (points with the same gpstime) have points with identical ReturnNumber
##   - Checking flag attributes... ✓
##   - Checking user data attribute...
##     🛈 26546448 points have a non 0 UserData attribute. This probably has a meaning
##  Checking the header
##   - Checking header completeness... ✓
##   - Checking scale factor validity... ✓
##   - Checking point data format ID validity... ✓
##   - Checking extra bytes attributes validity... ✓
##   - Checking the bounding box validity... ✓
##   - Checking coordinate reference system... ✓
##  Checking header vs data adequacy
##   - Checking attributes vs. point format... ✓
##   - Checking header bbox vs. actual content... ✓
##   - Checking header number of points vs. actual content... ✓
##   - Checking header return number vs. actual content... ✓
##  Checking coordinate reference system...
##   - Checking if the CRS was understood by R... ✓
##  Checking preprocessing already done 
##   - Checking ground classification... yes
##   - Checking normalization... no
##   - Checking negative outliers...
##     ⚠ 17514707 points below 0
##   - Checking flightline classification... yes
##  Checking compression
##   - Checking attribute compression...
##    -  Synthetic_flag is compressed
##    -  Keypoint_flag is compressed

# Seems all good for now

# Filtering by last returns

las2017 <- filter_poi(las2017, Classification == 2)

unique(las2017$Classification)

## [1] 2

plot(las2017)

# Based on the NY state website with the original datasets, the 2014 data and the combined 2017 data do not fully overlap. 
 
# I will need to clip one to the other to make sure I'm analyzing the same area. 

# This means I need to reproject one of the two to the other's projection first

# For 2014, current projection is EPSG:6347, for 2017, current projection is EPSG:2263

las2017_transf <- st_transform(las2017, 6347)

# Converting Z values to meters manually to allow for better comparisons between tiles. 

las2017_transf@data$Z <- las2017_transf@data$Z * 0.3048

print(las2017_transf)

## class        : LAS (v1.4 format 6)
## memory       : 491.7 Mb 
## extent       : 598507.8, 600046.9, 4497165, 4498702 (xmin, xmax, ymin, ymax)
## coord. ref.  : NAD83(2011) / UTM zone 18N 
## area         : 1.29 km²
## points       : 5.6 million points
## density      : 4.36 points/m²
## density      : 3.65 pulses/m²

# Decided to use a TIN initially to interpolate the LiDAR data for efficiency, but considering other methods and weighing pros and cons of different options
# https://r-lidar.github.io/lidRbook/dtm.html#sec-tin

DEM2014 <- rasterize_terrain(tile2014, res = 1, algorithm = tin())

## 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)

## Warning: Interpolation of 514 points failed because they are too far from
## ground points. Nearest neighbour was used but interpolation is weak for those
## points

Basic Plotting

plot(DEM2014)

# Writing cleaned and transformed datasets 

writeLAS(las2017_transf, "~/Desktop/Data Viz with R/R Spatial Project Alex/las2017_transf.las")

writeLAS(tile2014, "~/Desktop/Data Viz with R/R Spatial Project Alex/las2014.las")

DEM2017 <- rasterize_terrain(las2017_transf, res = 1, algorithm = tin())

## 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)Inverse distance weighting: [=================================-----------------] 67% (1 threads)Inverse distance weighting: [==================================----------------] 68% (1 threads)Inverse distance weighting: [==================================----------------] 69% (1 threads)Inverse distance weighting: [===================================---------------] 70% (1 threads)Inverse distance weighting: [===================================---------------] 71% (1 threads)Inverse distance weighting: [====================================--------------] 72% (1 threads)Inverse distance weighting: [====================================--------------] 73% (1 threads)Inverse distance weighting: [=====================================-------------] 74% (1 threads)Inverse distance weighting: [=====================================-------------] 75% (1 threads)Inverse distance weighting: [======================================------------] 76% (1 threads)Inverse distance weighting: [======================================------------] 77% (1 threads)Inverse distance weighting: [=======================================-----------] 78% (1 threads)Inverse distance weighting: [=======================================-----------] 79% (1 threads)Inverse distance weighting: [========================================----------] 80% (1 threads)Inverse distance weighting: [========================================----------] 81% (1 threads)Inverse distance weighting: [=========================================---------] 82% (1 threads)Inverse distance weighting: [=========================================---------] 83% (1 threads)Inverse distance weighting: [==========================================--------] 84% (1 threads)Inverse distance weighting: [==========================================--------] 85% (1 threads)Inverse distance weighting: [===========================================-------] 86% (1 threads)Inverse distance weighting: [===========================================-------] 87% (1 threads)Inverse distance weighting: [============================================------] 88% (1 threads)Inverse distance weighting: [============================================------] 89% (1 threads)Inverse distance weighting: [=============================================-----] 90% (1 threads)Inverse distance weighting: [=============================================-----] 91% (1 threads)Inverse distance weighting: [==============================================----] 92% (1 threads)Inverse distance weighting: [==============================================----] 93% (1 threads)Inverse distance weighting: [===============================================---] 94% (1 threads)Inverse distance weighting: [===============================================---] 95% (1 threads)Inverse distance weighting: [================================================--] 96% (1 threads)Inverse distance weighting: [================================================--] 97% (1 threads)Inverse distance weighting: [=================================================-] 98% (1 threads)Inverse distance weighting: [=================================================-] 99% (1 threads)Inverse distance weighting: [==================================================] 100% (1 threads)

## Warning: Interpolation of 2077 points failed because they are too far from
## ground points. Nearest neighbour was used but interpolation is weak for those
## points

plot(DEM2017)

# Write rasters to files to avoid having to remake 
#terra::writeRaster(DEM2017, "~/Desktop/Data Viz with R/R Spatial Project Alex//DEM2017.tif")

#terra::writeRaster(DEM2014, "~/Desktop/Data Viz with R/R Spatial Project Alex/DEM2014.tif")

DEM2017_2 <- rast("~/Desktop/Data Viz with R/R Spatial Project Alex/DEM2017.tif")
DEM2014_2 <- rast("~/Desktop/Data Viz with R/R Spatial Project Alex/DEM2014.tif")

DEM2014_aligned <- resample(DEM2014_2, DEM2017_2, method="bilinear")

DEM2017_aligned <- resample(DEM2017_2, DEM2014_2, method="bilinear")

plot(DEM2017_aligned)

plot(DEM2014_aligned)

# Finding overlap extent of the two 

dem_overlap <- intersect(ext(DEM2014_aligned), ext(DEM2017_aligned))

DEM2014_clip <- crop(DEM2014_aligned, dem_overlap)
DEM2017_clip <- crop(DEM2017_aligned, dem_overlap)
plot(DEM2017_clip)

plot(DEM2014_clip)

# Simple difference grid by just subtracting 2014 from 2017

dem_change <- DEM2017_clip - DEM2014_clip

plot(dem_change, bg='white')

# Looks off, the triangles from the TINs seem to be producing strange results for the 2017 data 

Troubleshooting

# The documentation mentions that edge artifacts are an issue with TINs- it certainly seems to be an issue with this data. Attempting to clip to valid data with buffers to cut out areas that are not ground from the DEMs before repeating the DEM change calculation

# Getting the coordinates from my final filtered LAS files

GroundCoords2017 <- st_as_sf(las2017_transf, coords=c('X', 'Y', 'Z'), crs=st_crs(las2017_transf))
GroundCoords2014 <- st_as_sf(tile2014, coords=c('X', 'Y', 'Z'), crs=st_crs(tile2014))

# Creating convex hulls around each as I would do in QGIS or ArcGIS Pro, since the coastline does not easily fit within bounding boxes

hull2017 <- st_convex_hull(GroundCoords2017)

hull2014 <- st_convex_hull(GroundCoords2014)

class(hull2017)

## [1] "sf"         "data.table" "data.frame"

# Convert the sf object to SpatVector
hull2017_v <- vect(hull2017)

hull2014_v <- vect(hull2014)

buffer2017_15<- terra::buffer(hull2017_v, width=15)

buffer2014_15 <- terra::buffer(hull2014_v, width=15)

#plot(buffer2014_15)
#plot(buffer2017_15)

DEM2017_final <- mask(DEM2017_clip, buffer2017_15)
DEM2014_final <- mask(DEM2014_clip, buffer2014_15)

plot(DEM2017_final)

plot(DEM2014_final)

dem_overlap_fin <- intersect(ext(DEM2017_final), ext(DEM2014_final))

DEM2017finOverlap <- mask(DEM2017_final, DEM2014_final)

DEM2014finOverlap <- mask(DEM2014_final, DEM2017_final)

plot(DEM2017finOverlap)

plot(DEM2014finOverlap)

# Rewrite rasters with final versions

#terra::writeRaster(DEM2017finOverlap, "~/Desktop/Data Viz with R/R Spatial Project Alex//DEM2017finOverlap_15.tif")
#terra::writeRaster(DEM2014finOverlap, "~/Desktop/Data Viz with R/R Spatial Project Alex/DEM2014finOverlap_15.tif")

# Load rasters back in to avoid rerunning earlier code

DEM2017finOverlap <- rast("~/Desktop/Data Viz with R/R Spatial Project Alex//DEM2017finOverlap_15.tif")
DEM2014finOverlap <- rast("~/Desktop/Data Viz with R/R Spatial Project Alex//DEM2014finOverlap_15.tif")

plot(DEM2014finOverlap)

plot(DEM2017finOverlap)

dem_change_f <- DEM2017finOverlap - DEM2014finOverlap

plot(dem_change_f, bg='white')

Final Calculations and Final Plot Development

# Using Wang et al. 2023's equation to calculate organic matter inventory (OMI)

# y = 0.15x + 0.12
# x = organic matter, y = accretion rate
# y - 0.12 = x
# (y - 0.12)/0.15 = x
# DEM2017 - DEM2014 should yield accretion rate

OMI <- (dem_change_f - 0.12)/0.15
plot(OMI)

# Making map scales more consistent for the DEMs before making final plots

min(DEM2017finOverlap$Z)

## class       : SpatRaster 
## dimensions  : 1335, 1493, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 598507, 6e+05, 4497165, 4498500  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83(2011) / UTM zone 18N (EPSG:6347) 
## source(s)   : memory
## varname     : DEM2017finOverlap_15 
## name        :   min 
## min value   : -0.86 
## max value   :  7.31

max(DEM2017finOverlap$Z)

## class       : SpatRaster 
## dimensions  : 1335, 1493, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 598507, 6e+05, 4497165, 4498500  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83(2011) / UTM zone 18N (EPSG:6347) 
## source(s)   : memory
## varname     : DEM2017finOverlap_15 
## name        :   max 
## min value   : -0.86 
## max value   :  7.31

min(DEM2014finOverlap$Z)

## class       : SpatRaster 
## dimensions  : 1335, 1493, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 598507, 6e+05, 4497165, 4498500  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83(2011) / UTM zone 18N (EPSG:6347) 
## source(s)   : memory
## varname     : DEM2014finOverlap_15 
## name        :   min 
## min value   : -1.14 
## max value   :  7.12

max(DEM2014finOverlap$Z)

## class       : SpatRaster 
## dimensions  : 1335, 1493, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 598507, 6e+05, 4497165, 4498500  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83(2011) / UTM zone 18N (EPSG:6347) 
## source(s)   : memory
## varname     : DEM2014finOverlap_15 
## name        :   max 
## min value   : -1.14 
## max value   :  7.12

# min from 2014 is lower, min will be -1.14
# max from 2017 is higher, max will be 7.31

min_z <- -1.14
max_z <- 7.31

breaks <- seq(min_z, max_z, length.out = 12)

# Testing with basic terra plots

terra::plot(DEM2014finOverlap, breaks=breaks, main = "2014 Ground Elevation")

terra::plot(DEM2017finOverlap, breaks=breaks, main="2017 Ground Elevation")

# 2014 Static Map

marsh_palette = c("#006994", "#2E8B57", "yellowgreen", "tan", "grey", "lightgrey")


tm_shape(DEM2014finOverlap) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Ground Elevation (Meters)"), col.scale = tm_scale(breaks = breaks, values = marsh_palette)) + 
  tm_title("Rulers Bar Hassock: Digital Elevation Model,2014")

tm_shape(DEM2017finOverlap) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Ground Elevation (Meters)"), col.scale = tm_scale(breaks = breaks, values = marsh_palette)) + tm_title("Rulers Bar Hassock: Digital Elevation Model, 2017")

tm_shape(dem_change_f) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Change in Elevation, 2014-2017"), col.scale = tm_scale(values = "terrain")) 

min(dem_change_f)

## class       : SpatRaster 
## dimensions  : 1335, 1493, 1  (nrow, ncol, nlyr)
## resolution  : 1, 1  (x, y)
## extent      : 598507, 6e+05, 4497165, 4498500  (xmin, xmax, ymin, ymax)
## coord. ref. : NAD83(2011) / UTM zone 18N (EPSG:6347) 
## source(s)   : memory
## varname     : DEM2017finOverlap_15 
## name        :   min 
## min value   : -1.75 
## max value   :  2.08

# Vast majority of values are between negative 1 and 1, changing min and max to better represent that range

min_change = -1.75
max_change = 2.1

change_breaks <- seq(min_change, max_change, length.out = 20)

tm_shape(dem_change_f) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Change in Elevation, 2014-2017"), col.scale = tm_scale(breaks= change_breaks, values = c("darkblue", "tan", "brown"))) + tm_title("Rulers Bar Hassock: Change in Elevation, 2014-17")

## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

tm_shape(OMI) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Est. Change in OMI, 2014-2017"), col.scale = tm_scale(values = c("brown", "darkgrey","yellowgreen", "darkgreen" ))) + tm_title("Rulers Bar Hassock: Organic Mattery Inventory (OMI), 2014-17")

## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

Final Plots

tm_shape(DEM2014finOverlap) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Ground Elevation (Meters)"), col.scale = tm_scale(breaks = breaks, values = marsh_palette)) + 
  tm_title("Rulers Bar Hassock: Digital Elevation Model,2014")

tm_shape(DEM2017finOverlap) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Ground Elevation (Meters)"), col.scale = tm_scale(breaks = breaks, values = marsh_palette)) + tm_title("Rulers Bar Hassock: Digital Elevation Model, 2017")

tm_shape(dem_change_f) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Change in Elevation, 2014-2017"), col.scale = tm_scale(breaks= change_breaks, values = c("darkblue", "tan", "brown"))) + tm_title("Rulers Bar Hassock: Change in Elevation, 2014-17")

## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

tm_shape(OMI) + 
  tm_raster(col = "Z", col.legend = tm_legend(title="Est. Change in OMI, 2014-2017"), col.scale = tm_scale(values = c("brown", "darkgrey","yellowgreen", "darkgreen" ))) + tm_title("Rulers Bar Hassock: Organic Mattery Inventory, 2014-17")

## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.