Weighted interpolation of sparse STILT footprints

gaussian_weighted_interpolation.r

#!/usr/bin/env Rscript
# Ben Fasoli
#
# Interpolation utilities for STILT footprints to predict footprints for
# finely gridded receptors using coarse STILT-calculated footprints for
# intermediate receptor locations.
#
# Uses weights derived from 2d Gaussian kernels to average footprints from
# the surrounding spatial area
#
# TODO: expand from 2d (time-integrated) footprints to 3d
# TODO: footprint prediction loss as a function of simulation grid density

library(raster)
library(tidyverse)

FOOTPRINT_PATH <- 'out/footprints'


#' Identifies path to given footprint file, raising an error if not found
#' @param time POSIXct receptor timestamp
#' @param long numeric x position of receptor
#' @param lati numeric y position of receptor
#' @param zagl numeric z position of receptor
#' @param path path to access footprint .nc files
find_footprint_file <- function(time, long, lati, zagl, path) {
  file <- paste(
    strftime(time, tz = 'UTC', format = '%Y%m%d%H'),
    long,
    lati,
    zagl,
    'foot.nc',
    sep = '_'
  )
  file <- file.path(path, file)
  stopifnot(file.exists(file))
  file
}


#' Identify index of footprint in given footprint meta dataframe
#' @param meta data.frame as returned by `make_footprint_meta(...)`
#' @param time POSIXct receptor timestamp
#' @param x numeric x position of receptor
#' @param y numeric y position of receptor
find_meta_index <- function(meta, time, x, y) {
  idx <- with(meta, which(time == time & long == x & lati == y))
  ifelse(length(idx) > 0, idx, 0)
}


#' Read footprint from file
#' @param file path to footprint .nc file
get_footprint_from_file <- function(file) {
  suppressWarnings(readAll(raster(file)))
}


#' Create footprint metadata record for array of file paths
#' @param files path to footprint .nc files
make_footprint_meta <- function(files) {
  base_files <- basename(files)
  meta_split <- strsplit(base_files, '_')
  tibble(
    file = files,
    time = as.POSIXct(sapply(meta_split, function(x) x[1]),
                      tz = 'UTC', format = '%Y%m%d%H'),
    long = as.numeric(sapply(meta_split, function(x) x[2])),
    lati = as.numeric(sapply(meta_split, function(x) x[3])),
    zagl = as.numeric(sapply(meta_split, function(x) x[4])),
    stringsAsFactors = F
  )
}


#' Make object containing shared grid definitions
#' @param r template raster from which to extract grid data
make_grid <- function(r) {
  xy <- coordinates(r)
  xy_res <- res(r)
  g <- list(
    x = sort(unique(xy[, 1])),
    y = sort(unique(xy[, 2])),
    dx = xy_res[1],
    dy = xy_res[2]
  )
  g$nx <- length(g$x)
  g$ny <- length(g$y)
  g$xi <- 1:g$nx
  g$yi <- 1:g$ny
  g
}


#' Make gaussian weight matrix from grid definitions
#' @param grid shared grid definitions as returned by `make_grid(...)`
#' @param r raster used to size output weight matrix
#' @param k kernel bandwidth scaling factor in n_grid_cell units, defaults to 1
make_weights <- function(grid, r, k = 1) {
  w <- list()
  w$fw <- focalWeight(r, k * max(c(grid$dx, grid$dy)), type = 'Gauss')
  w$nx <- ncol(w$fw)
  w$ny <- nrow(w$fw)
  w$center_xi <- (w$nx + 1) / 2
  w$center_yi <- (w$ny + 1) / 2
  w$dx <- 1:w$nx - w$center_xi
  w$dy <- 1:w$ny - w$center_yi
  w
}


# Fetch data -------------------------------------------------------------------
files <- dir(FOOTPRINT_PATH, full.names = T)


template <- get_footprint_from_file(files[1])
grid <- make_grid(template)
wgt <- make_weights(grid, template, k = 1)


# Footprint array small for testing so cache into memory
footprint <- make_footprint_meta(files) %>%
  mutate(foot = lapply(file, get_footprint_from_file),
         foot_interp = NA)

# Sparsely populate dummy footprint grid
spacing <- 2  # sparse receptors placed every `spacing` grid points
x_rem <- footprint$long %% (spacing * grid$dx)
y_rem <- footprint$lati %% (spacing * grid$dy)
x_rem <- signif(x_rem, 2)  # C-level rounding error from source raster
y_rem <- signif(y_rem, 2)  # C-level rounding error from source raster
x_rem_min <- min(x_rem)
y_rem_min <- min(y_rem)
footprint_sparse <- filter(footprint,
                           x_rem == x_rem_min,
                           y_rem == y_rem_min)


# Where are the sparse receptors?
xy <- footprint_sparse[, c('long', 'lati')]
with(footprint, plot(long, lati))
with(footprint_sparse, points(long, lati, col = 'red', pch = 16))


# time and zagl shared across receptors
time <- as.POSIXct('2015-06-18 22:00:00', tz = 'UTC')
zagl <- 5


# xi and yi are grid indices defining the grid location we are predicting
for (xi in 1:grid$nx) {
  for (yi in 1:grid$ny) {
    message('xi: ', xi, '    yi: ', yi)
    
    xiloc <- grid$x[xi]
    yiloc <- grid$y[yi]
    
    sparse_idx <- find_meta_index(footprint_sparse, time, xiloc, yiloc)
    
    if (sparse_idx) {
      idx <- find_meta_index(footprint, time, xiloc, yiloc)
      footprint$foot_interp[idx] <- list(footprint_sparse$foot[[sparse_idx]])
      next
    }
    
    # track total weights applied to scale end weighted average
    fwi_total <- 0
    fiisw_total <- raster(matrix(0), template = template)
    
    # wxi and wyi are indices of the focal weight matrix used to identify
    # footprints within the sparse grid used to calculate a weighted average
    for (wxi in 1:wgt$nx) {
      for (wyi in 1:wgt$ny) {
        dx <- wgt$dx[wxi]
        dy <- wgt$dy[wyi]
        xii <- xi + dx
        yii <- yi + dy
        
        if (xii < 1 | xii > grid$nx | yii < 1 | yii > grid$ny) next
        
        xiiloc <- grid$x[xii]
        yiiloc <- grid$y[yii]
        
        sparse_idx <- find_meta_index(footprint_sparse, time, xiiloc, yiiloc)
        if (!sparse_idx) next
        
        fii <- footprint_sparse$foot[[sparse_idx]]
        fiis <- shift(fii, x = -dx * grid$dx, y = -dy * grid$dy)
        fwi <- wgt$fw[wxi, wyi]
        fiisw <- fiis * fwi
        fwi_total <- fwi_total + fwi
        fiisw_total <- fiisw_total + extend(crop(fiisw, fiisw_total),
                                            fiisw_total, value = 0)
      }
    }
    
    fiisw_wgt_avg <- fiisw_total / fwi_total
    
    idx <- find_meta_index(footprint, time, xiloc, yiloc)
    footprint$foot_interp[idx] <- list(fiisw_wgt_avg)
  }
}