Project

Assessing Salt Marsh Change in Jamaica Bay Using LiDAR Data (2014–2017)

In this report, I will be focusing on working with LiDAR data to analyze whether the salt marsh in Jamaica Bay, Queens goes up or down between 2014 and 2017. The LiDAR data will be obtained from the National Oceanic and Atmospheric Administration (NOAA). This will be done by processing point cloud data into DEM maps and comparing the differences in elevation data from both years. This will help me determine whether the marsh is healthy and growing or slowly eroding over time.

library(sf)

Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(lidR)

Attaching package: 'lidR'

The following object is masked from 'package:sf':

    st_concave_hull

library(RCSF)
library(terra)

terra 1.9.1

Attaching package: 'terra'

The following object is masked from 'package:lidR':

    watershed

library(raster)

Loading required package: sp

Attaching package: 'sp'

The following object is masked from 'package:lidR':

    wkt

Attaching package: 'raster'

The following objects are masked from 'package:lidR':

    projection, projection<-

Processing, Classifying and Extracting Data

I read in the data files and reduced the number of points to improve processing speed. Afterwards, I classified the LiDAR point cloud into ground and non-ground points, isolating ground-level elevation from vegetation and other above-ground features.

# read las file
las_2014 <- readLAS("/Users/Ally/Desktop/LIDAR_Project/25170_2014.las")

Warning: There are 3874701 points flagged 'withheld'.

las_2017 <- readLAS("/Users/Ally/Desktop/LIDAR_Project/25170_2017.las")

# decreasing number of points 
las_2014 <- decimate_points(las_2014, homogenize(1, res = 1))
las_2017 <- decimate_points(las_2017, homogenize(1, res = 1))

# classify LiDAR point clouds into ground and non-ground points
las_2014 <- classify_ground(las_2014, algorithm = csf(time_step = 5))

Original dataset already contains 1132272 ground points. These points were reclassified as 'unclassified' before performing a new ground classification.

las_2017 <- classify_ground(las_2017, algorithm = csf(time_step = 5))

Original dataset already contains 162513 ground points. These points were reclassified as 'unclassified' before performing a new ground classification.

# filtering ground point classification 
ground_2014 <- filter_poi(las_2014, Classification == 2) 
ground_2017 <- filter_poi(las_2017, Classification == 2) 

Interpolating Data

Using the classified ground points, I generated a DEM for each year displaying the salt marsh elevation across Jamaica Bay. I then subtracted the two grids to produce a difference map showing elevation change between 2014 and 2017.

2014 Salt Marsh in Jamaica Bay

# clean and convert scattered LiDAR points into an elevation grid
dem_2014 <- rasterize_terrain(ground_2014, res = 1, algorithm = tin()) 

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)

plot(dem_2014, 
     main = "2014 Salt Marsh Elevation - Jamaica Bay", 
     xlab = "Longitude", 
     ylab = "Latitude")

2017 Salt Marsh in Jamaica Bay

# clean and convert scattered LiDAR points into an elevation grid
dem_2017 <- rasterize_terrain(ground_2017, res = 1, algorithm = tin()) 
crs(dem_2017) <- "EPSG:32618" 

# reproject 2017 to match 2014 coordinate system
dem_2017 <- project(dem_2017, crs(dem_2014))
dem_2017 <- resample(dem_2017, dem_2014)

plot(dem_2017, 
     main = "2017 Salt Marsh Elevation - Jamaica Bay", 
     xlab = "Longitude", 
     ylab = "Latitude")

Difference between Salt Marsh in Jamaica Bay

dem_difference <- dem_2017 - dem_2014
plot(dem_difference, 
     main = "Elevation Change in Jamaica Bay (2014-2017)", 
     xlab = "Longitude", 
     ylab = "Latitude")