LiDAR DEM Jamaica Bay Project

The goal of this project is to compare ground and habitat erosion damage after Hurricane Sandy in Jamaica Bay,and to analyse the damage particularly in Ruffle Bar a part of the Gateway National Recreation Area. This sample area contains the habitat type of North American Atlantic Low Salt Marshes, this area has experienced the most habitat loss/damage by both area and percent of total area (Meixler, 34).

Part 1: Source and download LiDAR data in an approximate year after Hurricane Sandy:

  • 2014 NOAA NGS Topobathy Lidar: Post-Sandy (SC to NY) (4800)
  • 2017 NYC Topobathy Lidar: New York (9306)

Notes: Both data sources are sourced from NOAA dataviewer portal here. In their dataviewer you’re able to choose the class of points, area of interest and other filters to help lessen the amount of points to download.

Debugging notes: I also downloaded the same data from USGS and I might switch because this data is giving me problems with the grid and extent, when I try different ways to clip the marsh area as a rectangle or polygon it just keeps showing the full map of the Ruffle Bar. Another issue is that when I use par() and plot() to display both map plots side by side, the 2014 map appears larger. Added opt_progress() to remove the dem catalog printing.

ctg_2017 <- readLAScatalog("./2017_nyc_topobathy")
ctg_2014 <- readLAScatalog("./2014_NGS_postSandy_topobathy")

opt_progress(ctg_2017) <- FALSE
opt_progress(ctg_2014) <- FALSE

opt_output_files method docs this method saved the tiff files into a folder where I merged the tiles using terra library mosaic() into a raster.

2017 Plot

if(length(list.files("./dem_tiles_2017", pattern = "\\.tif$")) == 0) {
  opt_output_files(ctg_2017) <- "./dem_tiles_2017/{*}_dem"
  rasterize_terrain(ctg_2017, res = 3, algorithm = tin())
}

tif_files <- list.files("./dem_tiles_2017", pattern = "\\.tif$", full.names = TRUE)
dem_merged_2017 <- mosaic(sprc(lapply(tif_files, rast)))
plot(dem_merged_2017, col = colorRampPalette(rev(brewer.pal(9, "Spectral")))(100), main = "Jamaica Bay DEM 2017")

2014 Plot

if(length(list.files("./dem_tiles_2014", pattern = "\\.tif$")) == 0) {
  opt_output_files(ctg_2014) <- "./dem_tiles_2014/{*}_dem"
  rasterize_terrain(ctg_2014, res = 3, algorithm = tin())
}

tif_files <- list.files("./dem_tiles_2014", pattern = "\\.tif$", full.names = TRUE)
dem_merged_2014 <- mosaic(sprc(lapply(tif_files, rast)))
plot(dem_merged_2014, col = colorRampPalette(rev(brewer.pal(9, "Spectral")))(100), main = "Jamaica Bay DEM 2014")

Part 2: Clip to Salt Marsh and Display side by side plots for 2014 and 2017

Working on different ways to clip the data to the outer edge of land as a rectangle or polygon to zoom into the salt marsh area

Notes: bbox_marsh is the bounding box set to WGS84 and reprojected to the las catalog’s (in US ft) CRS using sf’s st_transform(), then used to clip both LiDAR datasets to the Ruffle Bar marsh area, retaining only ground (class 2) and bathymetric (class 40) points via lidR’s filter_poi().

Debugging Notes: I chose ground points (class 2) for the bare marsh surface elevation and bathymetric points (class 40) to get the underwater bay bottom, I wanted to use these to detect the above and below water elevation change. In the Meixler article, the author points out that Ruffle Bar as a low-lying marsh platform surrounded by bay waters had the most errosion 18.13% of its total area was impacted but I can’t really see much of a difference. There is a 2005 LiDAR dataset that might be compelling to add if there’s a more stark comparison between these two years. Meixler’s study predicted that these marshes would recover in five years and the findings from the plots below suggest that with slightly higher elevations in the 2017 plot, the salt marshes could be recovering.

# Similar to prev lab where we made a boundary extent for the map and set the crs
bbox_marsh <- st_bbox(c(
  xmin = -73.880,
  xmax = -73.835,
  ymin = 40.585,
  ymax = 40.615
), crs = st_crs(4326)) |>
  st_as_sfc() |>
  st_transform(crs = st_crs(ctg_2017))

marsh_ext <- st_bbox(bbox_marsh)

# Clip data to marsh coordinates (didn't work)
# las_2017 <- readLAS(ctg_2017 , filter = "-keep_class 2")
# las_2017 <- filter_poi(las_2017, Classification %in% c(2, 40))
# las_2017 <- clip_rectangle(las_2017, 1018957, 154254, 1028688, 160827) 

# las_2014 <- readLAS(ctg_2014, filter = "-keep_class 2")
# las_2014 <- filter_poi(las_2014, Classification %in% c(2, 40))
# las_2014 <- clip_rectangle(las_2014, 1018957, 154254, 1028688, 160827)

# Another way to clip that doesn't work 
# 2017
las_2017 <- readLAS(ctg_2017@data$filename[1], filter = "-keep_class 2")
las_2017 <- filter_poi(las_2017, Classification %in% c(2, 40))
las_2017 <- clip_rectangle(las_2017, marsh_ext["xmin"], marsh_ext["ymin"], marsh_ext["xmax"], marsh_ext["ymax"])
# 2014
las_2014 <- readLAS(ctg_2014@data$filename[1], filter = "-keep_class 2")
las_2014 <- filter_poi(las_2014, Classification %in% c(2, 40))
las_2014 <- clip_rectangle(las_2014, marsh_ext["xmin"], marsh_ext["ymin"],marsh_ext["xmax"], marsh_ext["ymax"])

print(las_2017)
print(las_2014)

Part 3: Marsh comparison, 2014 and 2017 plots side by side

Notes: I used knnidw() over algorithm = tin() from the lidR doc online book (section 5.2) recommends it over TIN for small clipped areas because it handles edges better where TIN algorithm creates unrealistic triangles at the boundaries when there’s no buffer of extra points around the study area. Both rasters were aligned to one grid using terra library’s resample() and crop() before plotting side by side.

Debugging Notes: The problem is that the 2014 map appears larger than the 2017 map even after using resample() and crop() to align them, even after clipping and reprojecting they don’t overlap and appear to have different boundaries also.

dem_2017 <- rasterize_terrain(las_2017, res=5, algorithm = knnidw(k=10L, p=2))
dem_2014 <- rasterize_terrain(las_2014, res=5, algorithm = knnidw(k=10L, p=2))

dem_2017 <- project(dem_2017, "EPSG:32618")
dem_2014 <- project(dem_2014, "EPSG:32618")

dem_2014 <- resample(dem_2014, dem_2017, method = "bilinear")
shared   <- intersect(ext(dem_2017), ext(dem_2014))
dem_2017 <- crop(dem_2017, shared)
dem_2014 <- crop(dem_2014, shared)

pal <- colorRampPalette(rev(brewer.pal(9, "Spectral")))(100)
par(mfrow = c(1, 2))
plot(dem_2014, col = pal, main ="2014 Post-Sandy Marsh", legend = FALSE)
plot(dem_2017, col = pal, main ="2017 Marsh DEM")

par(mfrow = c(1, 1))

Looking at the data to see elevation change

Notes: The average elevation change -0.049 ft shows the marsh didn’t change much during early post-Sandy recovery, with the extremes is most likely due to it being in an area where sediments can be moved around more frequently.

Debugging Notes: Maybe I should choose a different sample area in Jamaica Bay with marshes at the shoreline or not all surrounded by water. I still need to complete the lost habitat and carbon loss part, maybe there will be different more varied results in Ruffle Bar marsh.

# elevation diff
dod <- dem_2017 - dem_2014
dod_values <- values(dod, na.rm = TRUE)

cat("Avg elevation change:", round(mean(dod_values), 3), "ft \n")
## Avg elevation change: -0.047 ft
cat("Max erosion:", round(min(dod_values), 3), "ft \n")
## Max erosion: -4.114 ft
cat("Max gained or accretion:", round(max(dod_values), 3), "ft \n")
## Max gained or accretion: 4.858 ft

Part 4: Elevation Change Map (Difference of the two DEMs)

pal_dod <- colorRampPalette(brewer.pal(11, "RdBu"))(100)
plot(dod,
     col  = pal_dod,
     main = "Elevation Change 2014-2017 [Red = Erosion and Blue = Accretion]")

Part 5: Amount of organic matter lost

Part 6: Interactive map

Sources:

Meixler, M. S. (2017). Assessment of Hurricane Sandy damage and resulting loss in ecosystem services in a coastal-urban setting. Ecosystem Services, 24, 28–46. https://doi.org/10.1016/j.ecoser.2016.12.009

Campbell, A., Wang, Y., Christiano, M., & Stevens, S. (2017). Salt Marsh Monitoring in Jamaica Bay, New York from 2003 to 2013: A Decade of Change from Restoration to Hurricane Sandy. Remote Sensing (Basel, Switzerland), 9(2), 131–131. https://doi.org/10.3390/rs9020131

Wang, C., Morgan, G. R., & Morris, J. T. (2023). Drone Lidar Deep Learning for Fine-Scale Bare Earth Surface and 3D Marsh Mapping in Intertidal Estuaries. Sustainability, 15(22), 15823. https://doi.org/10.3390/su152215823