NEON data processing with ggspatial

neon-ggspatial.R

library(conflicted)
library(hdf5r)
library(raster)
library(sf)
library(fs)
library(dplyr)
library(tidyr)
library(purrr)
library(fst)
library(vroom)
library(ggplot2)
library(ggspatial)
library(ggsflabel)

conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")

aopdir <- "~/projects/misc-data/neon-aop/NEON_spectrometer-orthorectified-surface-directional-reflectance---mosaic"
aop1 <- dir_ls(aopdir, glob = "*.h5")[1]
aop1_hf <- H5File$new(aop1, "r")
h5attributes(aop1_hf[["SERC/Reflectance/Reflectance_Data"]])
epsg <- as.numeric(aop1_hf[["SERC/Reflectance/Metadata/Coordinate_System/EPSG Code"]][])
aop1_hf$close_all()

field_data <- dir_ls("data", glob = "*fieldData*", recurse = TRUE) %>%
  vroom()

# Field data is from other script!
plotmap <- field_data %>%
  filter(siteID == "SERC") %>%
  distinct(plotID, plotType, decimalLatitude, decimalLongitude) %>%
  st_as_sf(coords = c("decimalLongitude", "decimalLatitude"), crs = 4326)
plotmap_t <- st_transform(plotmap, epsg)
pt2 <- plotmap_t %>%
  bind_cols(as_tibble(st_coordinates(plotmap_t))) %>%
  as_tibble() %>%
  mutate(floorX = floor(X / 1000) * 1000,
         floorY = floor(Y / 1000) * 1000)

make_rgb_raster <- function(f) {
  hf1 <- H5File$new(f, "r")
  on.exit(hf1$close_all(), add = TRUE)
  sr <- hf1[["SERC"]][["Reflectance"]]
  meta <- sr[["Metadata"]]
  epsg <- as.numeric(meta[["Coordinate_System"]][["EPSG Code"]][])
  refl <- sr[["Reflectance_Data"]]
  ext <- h5attr(refl, "Spatial_Extent_meters")
  irgb <- c(62, 35, 10)
  rgbmat <- refl[irgb,,]
  rgbmat[rgbmat < 0] <- NA
  rrgb <- brick(aperm(rgbmat, c(3, 2, 1)),
                ext[1], ext[2], ext[3], ext[4],
                crs = CRS(paste0("+init=epsg:", epsg)))
  rrgb
}

aopfiles <- dir_ls(aopdir, glob = "*.h5")
rasters <- unname(lapply(aopfiles, make_rgb_raster))
combined <- do.call(merge, rasters)

myext <- extent(364000, 366000, 4304000, 4307000)
csub <- crop(combined, myext) / 10000
csub[csub > 1] <- NA
csub_small <- aggregate(csub, factor = 20)

plotmap_t_sub <- st_crop(plotmap_t, csub_small)

st_envelope <- function(x, dist) {
  bufs <- st_buffer(x, dist)
  bboxs <- lapply(bufs, st_bbox)
  squares <- lapply(bboxs, st_as_sfc)
  result <- do.call(c, squares)
  st_crs(result) <- st_crs(x)
  result
}

squares <- plotmap_t_sub %>%
  mutate(geometry = st_envelope(geometry, 10))

aop_plot_map <- ggplot() +
  layer_spatial(stretch(csub_small, maxq = 0.995)) +
  geom_sf(data = squares, color = "red", fill = "red") +
  geom_sf_label_repel(aes(label = plotID), plotmap_t_sub, seed = 10) +
  theme_bw() +
  theme(axis.title = element_blank(),
        panel.grid = element_blank())
ggsave("figures/serc-aop-ggplot.png", aop_plot_map,
       width = 5.8, height = 6.65, units = "in", dpi = 300)

zoombox <- st_bbox(c(xmin = -76.565, xmax = -76.557,
                     ymin = 38.894, ymax = 38.899),
                   crs = 4326) %>%
  st_as_sfc() %>%
  st_transform(st_crs(csub_small)) %>%
  st_bbox()
zoomext <- extent(zoombox[c(1, 3, 2, 4)])
csub_zoom <- crop(csub_small, zoomext)
aop_plot_map_zoom <- ggplot() +
  layer_spatial(stretch(csub_zoom, maxq = 0.995)) +
  geom_sf(data = st_crop(squares, zoombox), color = "red", fill = "red") +
  geom_sf_label_repel(aes(label = plotID), st_crop(squares, zoombox),
                      force = 50, seed = 10) +
  theme_bw() +
  theme(axis.title = element_blank(),
        panel.grid = element_blank())
ggsave("figures/serc-aop-map-zoom.png", aop_plot_map_zoom,
       width = 4.9, height = 4.1, units = "in", dpi = 300)

# Just one plot
zoombox <- st_bbox(c(xmin = -76.5617, xmax = -76.5606,
                     ymin = 38.8972, ymax = 38.8979),
                   crs = 4326) %>%
  st_as_sfc() %>%
  st_transform(st_crs(csub_small)) %>%
  st_bbox()
zoomext <- extent(zoombox[c(1, 3, 2, 4)])
csub_zoom <- crop(csub_small, zoomext)
aop_plot_map_zoom <- ggplot() +
  layer_spatial(stretch(csub_zoom, maxq = 0.995)) +
  geom_sf(data = st_crop(squares, zoombox), color = "red", fill = alpha("red", 0.1)) +
  geom_sf_label_repel(aes(label = plotID), st_crop(squares, zoombox),
                      force = 50, seed = 10) +
  theme_bw() +
  theme(panel.grid = element_blank()) +
  labs(x = "UTM 18N Easting (m)", y = "UTM 18N Northing (m)") +
  coord_sf(datum = st_crs(squares))
## aop_plot_map_zoom
ggsave("figures/serc-aop-map-zoommax.png", aop_plot_map_zoom,
       width = 4.9, height = 4.1, units = "in", dpi = 300)

dat <- read_fst("data/extracted-serc-spectra.fst") %>%
  as_tibble() %>%
  filter(plotID == "SERC_001") %>%
  group_by(wavelengths) %>%
  mutate(uid = row_number()) %>%
  ungroup() %>%
  filter(wavelengths >= 400) %>%
  mutate(value = case_when(wavelengths > 1700 & wavelengths < 1950 ~ NA_real_,
                           wavelengths > 1350 & wavelengths < 1400 ~ NA_real_,
                           TRUE ~ value))

serc_onespec_plot <- ggplot(dat) +
  aes(x = wavelengths, y = value / 10000, group = uid) +
  geom_line(color = alpha("gray20", 0.1)) +
  theme_bw() +
  theme(panel.grid = element_blank()) +
  labs(x = "Wavelength (nm)", y = "Reflectance [0,1]")
ggsave("figures/serc-onespec-plot.png", serc_onespec_plot,
       width = 5, height = 3.3, units = "in", dpi = 300)


# NEON plots are 20 x 20 meters.
# We create a 10 meter buffer around the point
plot_buffers <- st_buffer(plotmap_t, 10)

extract_coords <- function(r, x) {
  cell_list <- extract(r, x, cellnumbers = TRUE)
  indices <- lapply(cell_list, function(x) {
    if (is.null(x)) return(x)
    cbind(row = rowFromCell(r, x[, "cell"]),
          col = colFromCell(r, x[, "cell"]))
  })
  indices
}
buffer_ind_list <- lapply(rasters, extract_coords, x = plot_buffers)

get_spectra <- function(indlist, file) {
  hf <- H5File$new(file, "r")
  on.exit(hf$close_all(), add = TRUE)
  sr <- hf[["SERC"]][["Reflectance"]]
  meta <- sr[["Metadata"]]
  wavelengths <- meta[["Spectral_Data"]][["Wavelength"]][]

  hrefl <- sr[["Reflectance_Data"]]
  outlist <- lapply(indlist, function(x) {
    if (is.null(x)) return(NULL)
    out <- matrix(NA_real_, length(wavelengths), nrow(x))
    for (i in seq_len(nrow(x))) {
      # NOTE: Rows and columns are switched here because we used transpose on
      # the raster before.
      out[, i] <- hrefl[, x[i, 2], x[i, 1]]
    }
    out
  })
  outlist
}

spectra_list <- Map(get_spectra, buffer_ind_list, aopfiles)
pt2_spectra <- pt2 %>%
  mutate(spectra = spectra_list)

pt_inorder <- plotmap_t %>%
  bind_cols(as_tibble(st_coordinates(plotmap_t))) %>%
  as_tibble() %>%
  mutate(floorX = floor(X / 1000) * 1000,
         floorY = floor(Y / 1000) * 1000) %>%
  select(-geometry, -X, -Y)

pt_spec <- pt_inorder %>%
  left_join(pt2_spectra)

pt_long <- pt_spec %>%
  unnest(spectra) %>%
  filter(map_lgl(spectra, ~!is.null(.))) %>%
  mutate(spectra = map(spectra, as_tibble, .name_repair = "unique"),
         wavelengths = list(wavelengths)) %>%
  unnest(c(wavelengths, spectra))

pt_long2 <- pt_long %>%
  pivot_longer(starts_with("..."), names_to = "id", values_to = "value")

write_fst(pt_long2, "./data/extracted-serc-spectra.fst")

specplot <- pt_long2 %>%
  mutate(value = value / 10000) %>%
  group_by(plotID, wavelengths) %>%
  summarize(Mean = mean(value, na.rm = TRUE),
            lo = quantile(value, 0.1, na.rm = TRUE),
            hi = quantile(value, 0.9, na.rm = TRUE)) %>%
  ggplot() +
  aes(x = wavelengths, y = Mean,
      ymin = lo, ymax = hi) +
  geom_ribbon(fill = "gray80") +
  geom_line() +
  facet_wrap(vars(plotID)) +
  coord_cartesian(ylim = c(0, 0.75)) +
  labs(x = "Wavelength (nm)", y = "Reflectance [0,1]") +
  theme_bw() +
  theme(panel.grid = element_blank())
ggsave("figures/neon-extracted-spectra.png", specplot,
       width = 6.5, height = 4, units = "in", dpi = 300)