Lab_NC_Analysis_R

Load every session — paste these at the TOP of your script

library(terra) # raster and vector analysis 

terra 1.9.11

library(sf) # simple features vector data 

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

library(ggplot2) # publication-quality plots 
library(tidyterra) # ggplot2 geoms for terra rasters 

Attaching package: 'tidyterra'

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

    filter

library(patchwork) # combine multiple ggplot panels 

Attaching package: 'patchwork'

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

    area

library(leaflet) # interactive web maps 
library(dplyr) # data wrangling 

Attaching package: 'dplyr'

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

    intersect, union

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

    filter, lag

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

    intersect, setdiff, setequal, union

library(classInt) # classification breaks 
library(viridis) # color palettes 

Loading required package: viridisLite

library(ggnewscale) # multiple fill scales in one ggplot 
library(htmlwidgets) # save leaflet as standalone HTML

── CONSTANTS ───────────────────────────────────────────────────

ft_to_m <- 0.3048006096

── LOAD DEMs FROM VRT ──────────────────────────────────────────

dem_2020_raw <- rast("data/dem_2020_noaa.vrt") 
dem_2001_raw <- rast("data/dem_2001_ncfmp.vrt")

── STUDY AREA ──────────────────────────────────────────────────

aoi_wgs <- st_as_sfc(st_bbox(c(xmin = -76.1, xmax = -75.4,
                               ymin =  35.0, ymax =  35.6),
                             crs = 4326))
aoi_nc_stateplane <- st_transform(aoi_wgs, 6543)
aoi_utm           <- st_transform(aoi_wgs, 32618)

── AGGREGATE, CROP, CONVERT UNITS ──────────────────────────────

dem_2020_agg     <- aggregate(dem_2020_raw, fact = 5, fun = "mean")

|---------|---------|---------|---------|
=========================================
                                          

dem_2020_cropped <- crop(dem_2020_agg, vect(aoi_nc_stateplane))

|---------|---------|---------|---------|
=========================================
                                          

dem_2020_m       <- dem_2020_cropped * ft_to_m

|---------|---------|---------|---------|
=========================================
                                          

names(dem_2020_m) <- "elev_2020_m"
dem_2001_agg     <- aggregate(dem_2001_raw, fact = 5, fun = "mean")

|---------|---------|---------|---------|
=========================================
                                          

dem_2001_cropped <- crop(dem_2001_agg, vect(aoi_nc_stateplane))
dem_2001_m       <- dem_2001_cropped * ft_to_m
names(dem_2001_m) <- "elev_2001_m"

── REPROJECT BOTH TO UTM ZONE 18N ──────────────────────────────

dem_2020 <- project(dem_2020_m, "EPSG:32618")

|---------|---------|---------|---------|
=========================================
                                          

dem_2020 <- crop(dem_2020, vect(aoi_utm))

|---------|---------|---------|---------|
=========================================
                                          

dem_2001 <- project(dem_2001_m, "EPSG:32618")
dem_2001 <- crop(dem_2001, vect(aoi_utm))

cat("2020 DEM resolution (m):", res(dem_2020), "\n")

2020 DEM resolution (m): 4.999724 4.999724 

cat("2001 DEM resolution (m):", res(dem_2001), "\n")

2001 DEM resolution (m): 18.28417 18.28417 

# Resample 2001 UP to match the 2020 grid (5m is the better resolution to keep)
dem_2001_matched <- resample(dem_2001, dem_2020, method = "bilinear")

|---------|---------|---------|---------|
=========================================
                                          

DEM difference: positive = gain, negative = loss

elev_change <- dem_2020 - dem_2001_matched

|---------|---------|---------|---------|
=========================================
                                          

names(elev_change) <- "elev_change_m"

Statistical summary of change surface

cat("Elevation change summary (m, 2001 → 2020):\n")

Elevation change summary (m, 2001 → 2020):

print(global(elev_change, fun = c("min", "max", "mean", "sd"), na.rm = TRUE))

                    min      max     mean       sd
elev_change_m -18.12684 6.767681 -1.16616 1.896224

Classify change into five categories

elev_change_cat <- classify(
  elev_change,
  rcl = matrix(
    c(-Inf, -2,   1,
      -2,   -0.5, 2,
      -0.5,  0.5, 3,
      0.5,  2,   4,
      2,   Inf,  5),
    ncol = 3, byrow = TRUE
  )
)

|---------|---------|---------|---------|
=========================================
                                          

names(elev_change_cat) <- "change_class"

Land mask: keep only cells where 2020 elevation is above 0.2 m

land_mask        <- ifel(dem_2020 > 0.2, 1, NA)

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

dem_2020_matched <- mask(dem_2020, land_mask)

|---------|---------|---------|---------|
=========================================
                                          

elev_change      <- mask(elev_change, land_mask)

|---------|---------|---------|---------|
=========================================
                                          

slr_scenarios <- c(0.5, 1.0, 2.0, 3.0)  # meters (IPCC AR6)

lapply loops over all four scenarios in one expression

dem_2020_matched is already masked to land and in meters

inundation_stack <- lapply(slr_scenarios, function(s) {
  mask_layer <- ifel(dem_2020_matched <= s, 1, NA)
  names(mask_layer) <- paste0("slr_", s, "m")
  mask_layer
})

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

Combine into a single multi-layer raster

inundation_stack <- rast(inundation_stack)
inundation_stack <- mask(inundation_stack, land_mask)

|---------|---------|---------|---------|
=========================================
                                          

Calculate inundated area (km²) for each scenario

areas <- sapply(seq_along(slr_scenarios), function(i) {
  n_cells <- global(inundation_stack[[i]], "sum", na.rm = TRUE)[[1]]
  n_cells * prod(res(inundation_stack[[i]])) / 1e6
})

data.frame(scenario_m = slr_scenarios, area_km2 = round(areas, 1))

  scenario_m area_km2
1        0.5     93.5
2        1.0    143.4
3        2.0    174.8
4        3.0    188.8

Calculate slope from the 2020 DEM (in degrees)

slope_r <- terrain(dem_2020_matched, v = "slope", unit = "degrees")

|---------|---------|---------|---------|
=========================================
                                          

Normalization function: scales any raster to the 0–1 range

normalize <- function(r) {
  mn <- global(r, "min", na.rm = TRUE)[[1]]
  mx <- global(r, "max", na.rm = TRUE)[[1]]
  (r - mn) / (mx - mn)
}

Normalize each layer

elev_norm   <- 1 - normalize(dem_2020_matched)  # inverted: low elev = high risk

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

slope_norm  <- normalize(slope_r)

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

change_norm <- normalize(abs(elev_change))       # magnitude only; direction handled by elev_norm

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

Weighted index (weights must sum to 1.0)

w_elev <- 0.50; w_change <- 0.30; w_slope <- 0.20

vulnerability <- (w_elev * elev_norm) + (w_change * change_norm) + (w_slope * slope_norm)

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

names(vulnerability) <- "vulnerability"
vulnerability <- mask(vulnerability, land_mask)

|---------|---------|---------|---------|
=========================================
                                          

cat("Vulnerability index range:", round(minmax(vulnerability), 3), "\n")

Vulnerability index range: 0.035 0.797 

Compute slope and aspect in radians (required by shade())

aspect_r  <- terrain(dem_2020_matched, v = "aspect", unit = "radians")

|---------|---------|---------|---------|
=========================================
                                          

slope_rad <- terrain(dem_2020_matched, v = "slope",  unit = "radians")

|---------|---------|---------|---------|
=========================================
                                          

angle = solar elevation above horizon; direction = compass azimuth of light

hillshade <- shade(slope_rad, aspect_r, angle = 40, direction = 330)

|---------|---------|---------|---------|
=========================================
                                          

hillshade <- clamp(hillshade, 0, 1, values = TRUE)

|---------|---------|---------|---------|
=========================================
                                          

hillshade <- mask(hillshade, land_mask)

|---------|---------|---------|---------|
=========================================
                                          

names(hillshade) <- "hillshade"

map_theme <- theme_void() +
  theme(
    plot.title      = element_text(size = 10, face = "bold", hjust = 0.5,
                                   margin = margin(b = 4)),
    plot.subtitle   = element_text(size = 7, hjust = 0.5, color = "grey40",
                                   margin = margin(b = 2)),
    legend.position   = "bottom",
    legend.title      = element_text(size = 7, face = "bold"),
    legend.text       = element_text(size = 6),
    legend.key.height = unit(0.3, "cm"),
    legend.key.width  = unit(0.8, "cm"),
    plot.margin       = margin(4, 4, 4, 4)
  )

p_elev <- ggplot() +
  # First fill scale: hillshade in grayscale
  geom_spatraster(data = hillshade) +
  scale_fill_gradientn(
    colors   = gray.colors(100, start = 0.1, end = 0.9),
    na.value = "white", guide = "none"
  ) +
  # Open a second fill scale for the elevation overlay
  new_scale_fill() +
  geom_spatraster(data = dem_2020_matched, alpha = 0.6) +
  scale_fill_viridis_c(
    option   = "mako", name = "Elevation (m)",
    na.value = "transparent", limits = c(0, 10), oob = scales::squish
  ) +
  labs(title    = "A — Elevation (2019–2020)",
       subtitle = "UTM Zone 18N / EPSG:32618 — meters") +
  map_theme

<SpatRaster> resampled to 500832 cells.
<SpatRaster> resampled to 500832 cells.

p_change <- ggplot() +
  geom_spatraster(data = elev_change) +
  scale_fill_distiller(
    palette  = "RdBu", direction = 1,
    name     = "Δ Elevation (m)", na.value = "white",
    limits   = c(-2, 2), oob = scales::squish
  ) +
  labs(title    = "B — Elevation Change (2001–2020)",
       subtitle = "Blue = gain  |  Red = loss  |  ±2 m range") +
  map_theme

<SpatRaster> resampled to 500832 cells.

Collapse four layers: first scenario at which each cell floods

first_flooded <- app(inundation_stack, fun = function(x) {
  if (all(is.na(x))) return(NA_real_)
  idx <- which(!is.na(x))[1]
  if (is.na(idx)) NA_real_ else slr_scenarios[idx]
})

|---------|---------|---------|---------|
=========================================
                                          

names(first_flooded) <- "first_flooded_m"

p_slr <- ggplot() +
  geom_spatraster(data = hillshade) +
  scale_fill_gradientn(
    colors = gray.colors(100, 0.1, 0.9), na.value = "white", guide = "none"
  ) +
  new_scale_fill() +
  geom_spatraster(data = first_flooded) +
  scale_fill_stepsn(
    colors = c("#08519c", "#3182bd", "#9ecae1", "#deebf7"),
    breaks = c(0.5, 1.0, 2.0, 3.0),
    name   = "First inundated\nat SLR scenario (m)",
    na.value = "transparent",
    guide  = guide_colorsteps(even.steps = FALSE)
  ) +
  labs(title    = "C — Inundation Scenarios",
       subtitle = "+0.5 / +1.0 / +2.0 / +3.0 m SLR (IPCC AR6)") +
  map_theme

<SpatRaster> resampled to 500832 cells.
<SpatRaster> resampled to 500832 cells.

p_vuln <- ggplot() +
  geom_spatraster(data = hillshade) +
  scale_fill_gradientn(
    colors = gray.colors(100, 0.1, 0.9), na.value = "white", guide = "none"
  ) +
  new_scale_fill() +
  geom_spatraster(data = vulnerability, alpha = 0.75) +
  scale_fill_viridis_c(
    option   = "inferno", name = "Vulnerability\nIndex (0–1)",
    na.value = "transparent", direction = -1
  ) +
  labs(title    = "D — Coastal Vulnerability Index",
       subtitle = "Elevation × Change × Slope (weighted)") +
  map_theme

<SpatRaster> resampled to 500832 cells.
<SpatRaster> resampled to 500832 cells.

figure <- (p_elev | p_change) / (p_slr | p_vuln) +
  plot_annotation(
    title   = "Outer Banks, NC — Coastal Change & Vulnerability Analysis",
    caption = paste0(
      "Data: NOAA NGS Topobathy Lidar 2019–2020 (Dataset 9468) & NCFMP Lidar 2001  |  ",
      "Projection: UTM Zone 18N (EPSG:32618)  |  SLR: IPCC AR6 (2021)  |  R / terra / ggplot2"
    ),
    theme = theme(
      plot.title   = element_text(size = 13, face = "bold", hjust = 0.5),
      plot.caption = element_text(size = 6, color = "grey50", hjust = 0.5)
    )
  )

print(figure)

ggsave(
  "outer_banks_vulnerability.png",
  plot = figure, width = 10, height = 8, dpi = 300, bg = "white"
)
cat("Figure saved to outer_banks_vulnerability.png\n")

Figure saved to outer_banks_vulnerability.png

Extract values — na.rm = FALSE keeps vectors the same length

vuln_vals   <- as.vector(values(vulnerability,    na.rm = FALSE))
elev_vals   <- as.vector(values(dem_2020_matched, na.rm = FALSE))
change_vals <- as.vector(values(elev_change,      na.rm = FALSE))

Keep only rows where all three layers have valid data

complete_idx <- complete.cases(vuln_vals, elev_vals, change_vals)

df <- data.frame(
  vulnerability = vuln_vals[complete_idx],
  elevation     = elev_vals[complete_idx],
  elev_change   = change_vals[complete_idx]
)

Classify into quartiles

df$vuln_quartile <- cut(
  df$vulnerability,
  breaks = quantile(df$vulnerability, probs = seq(0, 1, 0.25), na.rm = TRUE),
  labels = c("Q1 Low", "Q2 Moderate", "Q3 High", "Q4 Very High"),
  include.lowest = TRUE
)

Compute cell area OUTSIDE the pipe to avoid terra/dplyr conflict

cell_area_km2 <- prod(res(vulnerability)) / 1e6

summary_table <- df %>%
  group_by(vuln_quartile) %>%
  summarise(
    n_cells       = n(),
    mean_elev_m   = round(mean(elevation,   na.rm = TRUE), 2),
    mean_change_m = round(mean(elev_change, na.rm = TRUE), 3),
    area_km2      = round(n() * cell_area_km2, 2),
    .groups       = "drop"
  )

print(summary_table)

# A tibble: 4 × 5
  vuln_quartile n_cells mean_elev_m mean_change_m area_km2
  <fct>           <int>       <dbl>         <dbl>    <dbl>
1 Q1 Low        1703620        2.03        -0.206     42.6
2 Q2 Moderate   1703625        0.61        -0.187     42.6
3 Q3 High       1703613        0.68        -0.852     42.6
4 Q4 Very High  1703618        0.83        -4.26      42.6

Reproject to WGS84 — leaflet always uses lat/lon

vuln_wgs <- project(aggregate(vulnerability, fact = 5, fun = "mean"), "EPSG:4326")

|---------|---------|---------|---------|
=========================================
                                          

slr_wgs  <- project(aggregate(first_flooded, fact = 5, fun = "mean"), "EPSG:4326")

|---------|---------|---------|---------|
=========================================
                                          

Color palettes for each layer

pal_vuln <- colorNumeric(
  palette  = "inferno", domain = c(0, 1),
  na.color = "transparent", reverse = FALSE
)
pal_slr <- colorFactor(
  palette  = c("#08519c", "#3182bd", "#9ecae1", "#deebf7"),
  domain   = c(0.5, 1.0, 2.0, 3.0), na.color = "transparent"
)

interactive_map <- leaflet() %>%
  addProviderTiles("CartoDB.Positron",  group = "Basemap (Light)") %>%
  addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
  addRasterImage(x = vuln_wgs, colors = pal_vuln,
                 opacity = 0.7, group = "Vulnerability Index") %>%
  addLegend(position = "bottomleft", pal = pal_vuln, values = c(0, 1),
            title = "Vulnerability<br>Index", group = "Vulnerability Index") %>%
  addRasterImage(x = slr_wgs, colors = pal_slr,
                 opacity = 0.65, group = "SLR Inundation") %>%
  addLegend(position = "bottomright", pal = pal_slr,
            values    = c(0.5, 1.0, 2.0, 3.0),
            title     = "First Flooded<br>at SLR (m)", group = "SLR Inundation",
            labFormat = labelFormat(suffix = " m")) %>%
  addLayersControl(
    baseGroups    = c("Basemap (Light)", "Satellite"),
    overlayGroups = c("Vulnerability Index", "SLR Inundation"),
    options       = layersControlOptions(collapsed = FALSE)
  ) %>%
  addScaleBar(position = "bottomleft") %>%
  addMiniMap(toggleDisplay = TRUE, minimized = TRUE) %>%
  setView(lng = -75.75, lat = 35.25, zoom = 10)

Warning in colors(.): Some values were outside the color scale and will be
treated as NA

print(interactive_map)

Save as a portable, self-contained HTML file

htmlwidgets::saveWidget(
  widget        = interactive_map,
  file          = "outer_banks_interactive.html",
  selfcontained = TRUE
)
cat("Map saved. Open outer_banks_interactive.html in any browser.\n")

Map saved. Open outer_banks_interactive.html in any browser.