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 habitat type has experienced the most damage by both area and percent of total area (Meixler, 34).

Part 1: Study Area Map

Ruffle Bar Marsh is an uninhabited salt marsh located in Jamaica Bay, Queens, New York. It is a part of the Gateway National Recreation Area on the eastern portion of the bay. Ruffle Bar’s habitat consists of the North American Low Salt Marsh habitat made up of mostly cordgrass or Spartina alterniflora. As one of the bay’s interior marsh islands, Ruffle Bar has undergone many changes in edge erosion and marsh loss in recent decades that we can observe more closely in this project during namely Hurricane Sandy in 2012.

A bbox_marsh is transformed from shape coordinates into an sf object dataframe, this data is then displayed using leaflet’s setView().

# Define bbox and convert to WGS84 for leaflet
bbox_marsh <- st_bbox(c(
  xmin = -73.880,
  xmax = -73.835,
  ymin = 40.585,
  ymax = 40.615
), crs = st_crs(4326)) |>
  st_as_sfc()

bbox_marsh_sf <- st_sf(geometry = bbox_marsh)

leaflet() |> addProviderTiles("Esri.WorldImagery") |>
  setView(lng = -73.857715, lat = 40.599519, zoom = 16) |>
  addPolygons(data = bbox_marsh_sf,
              fill = FALSE,
              label = "Ruffle Bar Study Area")

Part 2: Source and download LiDAR data

The most readily available and thorough LiDAR .laz point cloud data available for the study area of Jamaica Bay were from the 2014 and 2017 LiDAR surveys.

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

Both data sources for 2014 and 2017 are sourced from NOAA. To reduce memory constraints, readLAScatalog() method is used instead of readLAS() to index the tile metadata as a catalog as opposed to loading all points into memory.

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

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

2017 Plot

This checks to see if dem_tiles_2017 already exists, if it doesn’t then it creates the directory and populates it with the rastarized .tif tiles.

opt_output_files method docs this method saves tiff files into a folder where tiles are merged using terra library mosaic() into a single raster.

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 = "Ruffle Bar, Jamaica Bay DEM 2017")

2014 Plot

This does the same as above but for dem_tiles_2014 rastarized .tif tiles.

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 = "Ruffle Bar, Jamaica Bay DEM 2014")

Part 3: Clip to Salt Marsh

A bbox_marsh is the bounding box set to WGS84 and reprojected to the catalog’s ftUS 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().

Ground points (class 2) represent the bare marsh surface elevation and bathymetric points (class 40) capture the underwater bay bottom to detect above and below water elevation change.

# Reproject bbox to catalog CRS (ftUS)
bbox_marsh_proj <- bbox_marsh |>
  st_transform(crs = st_crs(ctg_2017))

marsh_ext <- st_bbox(bbox_marsh_proj)

# 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 4: Marsh DEMs Side by Side

The 2014 and 2017 LiDAR datasets are from two survey missions with different tile boundaries, because of this after clipping to the same Ruffle Bar area the two rasters had slightly different extents and grids and one map appeared larger than the other.

To align the two datasets I created boundary variables for both LAS files to first be clipped to the same boundary using clip_rectangle() and then the DEMs were then reprojected to EPSG:32618, and aligned using resample() followed by ext() from terra’s library.

The methodsknnidw() is used over tin() as the lidR book (section 5.2) recommends it for small clipped areas without buffers. Using resample() followed by ext() is similar to the R equivalent of ArcGIS’s snap raster, where each pixel is forced into grid alignment before plotting side by side, after printing both have the same extents.

In Meixler’s article, the author identifies North American Atlantic Low Salt Marsh as the most impacted habitat type in Jamaica Bay with 18.13% of its total area damaged by Hurricane Sandy but there isn’t much of a notable a difference in these maps (Meixler,p.34). It could be compelling to see a 2005 LiDAR between these two years as a before and after comparison.

xmin <- 1022133
xmax <- 1026046
ymin <- 156406
ymax <- 158922

las_2017 <- clip_rectangle(las_2017, xmin, ymin, xmax, ymax)
las_2014 <- clip_rectangle(las_2014, xmin, ymin, xmax, ymax)

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))
## Inverse distance weighting: [==============================================----] 93% (6 threads)Inverse distance weighting: [===============================================---] 94% (6 threads)Inverse distance weighting: [===============================================---] 95% (6 threads)Inverse distance weighting: [================================================--] 96% (6 threads)Inverse distance weighting: [================================================--] 97% (6 threads)Inverse distance weighting: [=================================================-] 98% (6 threads)Inverse distance weighting: [=================================================-] 99% (6 threads)Inverse distance weighting: [==================================================] 100% (6 threads)

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

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

print(ext(dem_2017))
## SpatExtent : 596146.531345339, 597350.145551851, 4494519.03365111, 4495300.6211042 (xmin, xmax, ymin, ymax)
print(ext(dem_2014))
## SpatExtent : 596146.531345339, 597350.145551851, 4494519.03365111, 4495300.6211042 (xmin, xmax, ymin, ymax)

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

par(mfrow = c(1, 1))

Part 5: Elevation Change Statistics Results

A new raster or DEM of Difference is created by subtracting the 2014 DEM is from the 2017 DEM, cell by cell where each pixel shows how much elevation has changed as a list of numbers dod_values. The average avg_elevation_change, the maximum erosion value pixel max_erosion and the maximum gained in elevation or max_accretion.

The USGS’s measured the long term rate of accretion for Ruffle Bar of 4.9 mm per year (Wang H et al.,2017,Table 4). Using this number as a baseline we compare it to the dod_values observed elevation change and create a new var dod_rate_mmyr calculated by converting each pixel from meters to millimeters and dividing by 3 years to account for 2014 to 2017.

The marsh lost elevation at -15.71 mm per year between 2014 and 2017, was more than four times faster than the natural accretion rate of USGS’s rate of 4.9 mm per year (Wang H. et al., 2017, Table 4). This leaves a deficit of -20.61 mm per year, which suggests Ruffle Bar had not fully recovered from Hurricane Sandy by 2017, although the interior marsh platform showed signs of recovery, the outer edges continued to erode.

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

avg_elevation_change <- round(mean(dod_values), 3)
max_erosion <- round(min(dod_values), 3)
max_accretion <- round(max(dod_values), 3)

cat("Average elevation change:", avg_elevation_change, "m \n")
## Average elevation change: -0.047 m
cat("Max erosion:",  max_erosion, "m \n")
## Max erosion: -4.129 m
cat("Max accretion:", max_accretion, "m \n")
## Max accretion: 4.856 m

usgs_accretion_mmyr <- 4.9  # mm per year from Wang H Table 4
dod_rate_mmyr <- (mean(dod_values) * 1000) / 3 # mm per year over 2014-2017

cat("USGS's accretion rate:", usgs_accretion_mmyr, "mm per year\n")
## USGS's accretion rate: 4.9 mm per year
cat("Observed net change rate:", round(dod_rate_mmyr, 2), "mm per year\n")
## Observed net change rate: -15.71 mm per year
cat("Deficit:", round(dod_rate_mmyr - usgs_accretion_mmyr, 2), "mm per year\n")
## Deficit: -20.61 mm per year

Part 6: Elevation Change Map

The elevation change map maps the the DEM of Difference or DoD raster, where each pixel shows how much the ground surface change between 2014 and 2017. The color scale shows the change where red indicates erosion and blue indicates accretion using the RdBu palette from RColorBrewer.

This new map shows that the interior of the marsh has slight accretion seen in blue between 2014 and 2017, but the outer edges show consistent erosion in red, mostly by the southern and western edges. This pattern is consistent with the findings in the USGS Wang H. et al. 2017 report of edge erosion at Ruffle Bar as well as supports Meixler’s five year prediction that marsh interiors were recovering while the edges or boundary loss continued post storm.

pal_dod <- colorRampPalette(brewer.pal(11, "RdBu"))(100)
plot(dod, col = pal_dod, range = c(-2, 2), main  = "Ruffle Bar, Jamaica Bay\nElevation Change from 2014 to 2017 (m)\nRed = Erosion and Blue = Accretion")

Part 7: Carbon Loss and Gain Estimate Map

This carbon change map confirms what we found in the Difference of DEMs map and aligns with the results of edge erosion at Ruffle Bar between 2014 and 2017. These results show a measurable carbon loss concentrated along the southern and western marsh boundaries and where the interior marsh remains relatively carbon stable.

Carbon change is estimated using the biomass-height regression forumla from the Wang C et al.(2023) paper in Figure 8d where biomass (g/m sq) = 205.693 + 492.731 × Marsh Height(m). The elevation change from the DoD is multiplied by the regression slope of 492.731 to estimate biomass change delta_biomass (Wang et al., 2023, Figure 8d). This is then converted to carbon delta_carbon using a conversion factor of 0.45 for tidal salt marsh grasses, from the Coastal Blue Carbon Manual, in this case the dominant grass is cordgrass or Spartina alterniflora species in Jamaica Bay (Howard et al., 2014, p. 151).

Since the biomass change is calculated as a density (g/m sq), multiplying by cell area of 25m sq at 5m resolution and dividing by 1e6 or 1,000,000 converts it to kilograms per cell so that the final sum() represents actual total carbon lost or gained across the marsh rather than an average rate.

cell_area <- res(dod)[1] * res(dod)[2]# m sq per cell 25 in total

# biomass change per cell using slope of Wang regression figure 8d
delta_biomass <- 492.731 * dod # g/m sq change per cell
delta_carbon <- delta_biomass * 0.45 # carbon change
total_carbon <- delta_carbon * cell_area / 1e6  # converts to kg per cell

carbon_values <- values(total_carbon, na.rm = TRUE) # total carbon
carbon_lost <- carbon_values[carbon_values < 0] # carbon loss
carbon_gained <- carbon_values[carbon_values > 0] # carbon gain

cat("Total carbon lost (kg):", round(sum(carbon_lost), 2),"\n")
## Total carbon lost (kg): -22.39
cat("Total carbon gained (kg):", round(sum(carbon_gained), 2),"\n")
## Total carbon gained (kg): 14.93
cat("Net carbon change (kg):", round(sum(carbon_values), 2),"\n")
## Net carbon change (kg): -7.45

# carbon change map
plot(total_carbon, col= colorRampPalette(brewer.pal(11, "RdBu"))(100), main = "Ruffle Bar Carbon Change 2014-2017 (kg/cell)\n Red = Loss and Blue = Gain")

Part 8: Interactive Map with Layers

This interactive map has all of our maps in one mapview from from week 8’s assignment project as it handles multiple raster layer toggling than one leaflet map. Each layer can be toggled on or off in the map viewer. DEMs and DoD are reprojected to WGS84 for display.

dem_2014_wgs <- project(dem_2014, "EPSG:4326")
dem_2017_wgs <- project(dem_2017, "EPSG:4326")
dod_wgs <- project(dod, "EPSG:4326")

mapview(dem_2014_wgs,layer.name = "DEM 2014",col.regions = rev(brewer.pal(9, "Spectral"))) +
mapview(dem_2017_wgs,layer.name = "DEM 2017",col.regions = rev(brewer.pal(9, "Spectral"))) +
mapview(dod_wgs,layer.name = "Elevation Change 2014-2017",col.regions = brewer.pal(11, "RdBu"))

Future additions to this project:

  • Future work on this project could include adding a 2005 LiDAR dataset to show pre-Sandy conditions for a more complete before and after comparison.
  • Expanding the carbon analysis to running the same DoD analysis on other Jamaica Bay marsh islands like Four Sparrow and Big Egg could show whether Ruffle Bar’s edge erosion pattern is unique or occurs bay wide.
  • I could also overlay the accretion deficit against NOAA’s projected sea level rise data to estimate how long Ruffle Bar will remain a viable above water salt marsh.

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

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

Wang, H., Chen, Q., Hu, K., Snedden, G. A., Hartig, E. K., Couvillion, B. R., Johnson, C. L., & Orton, P. M. (2017). Numerical modeling of the effects of Hurricane Sandy and potential future hurricanes on spatial patterns of salt marsh morphology in Jamaica Bay, New York City (Open-File Report 2017–1016). U.S. Geological Survey. https://doi.org/10.3133/ofr20171016

Howard, J., Hoyt, S., Isensee, K., Pidgeon, E., & Telszewski, M. (Eds.). (2014). Coastal Blue Carbon: Methods for assessing carbon stocks and emissions factors in mangroves, tidal salt marshes, and seagrass meadows. Conservation International, Intergovernmental Oceanographic Commission of UNESCO, International Union for Conservation of Nature. https://www.thebluecarboninitiative.org/manual