Make interactive map by spatially averaging high resolution observations

map_spatial_average.r

#!/usr/bin/env Rscript

options(dplyr.summarise.inform = F)

library(data.table)
library(htmlwidgets)
library(leaflet)
library(tidyverse)


get_observations <- function() {
  # Returns time joined observations from 3_primary_keys_nulled table
  dirname <- '/uufs/chpc.utah.edu/common/home/lin-group9/measurements/data/gsv/tables'
  filename <- '3_primary_keys_nulled.csv'
  fread(file.path(dirname, filename), data.table = F) %>%
    mutate(time = as.POSIXct(time, tz = 'UTC', format = '%Y-%m-%d %H:%M:%S')) %>%
    select(c(-starts_with('qc_'), -polygon)) %>%
    filter(!is.na(longitude), !is.na(latitude)) %>%
    arrange(system_id, time)
}


get_grid_uid <- function(x, y) {
  # Returns a unique id for each sequentially distinct x,y pair, e.g.
  #   x  y  uid
  #   1  1  1
  #   1  2  2
  #   1  2  2
  #   2  2  3
  #   1  2  4
  grid_id <- paste(x, y, sep = ',')
  run <- rle(grid_id)
  rep.int(1:length(run$lengths), run$lengths)
}


grid_observations <- function(df) {
  df %>%
    # Set x, y position to 0.002deg cell centers
    mutate(longitude = round(trunc(longitude * 1e3 / 2) * 2 / 1e3 - 1e-3, 8),
           latitude = round(trunc(latitude * 1e3 / 2) * 2 / 1e3 + 1e-3, 8)) %>%
    # Grid data separately for each vehicle
    group_by(system_id) %>%
    group_modify(function(rows, key) {
      rows %>%
        arrange(time) %>%
        # Identify sequence of grid cells traveled
        group_by(grid_uid = get_grid_uid(longitude, latitude)) %>%
        # Remove receptors if duration out of range
        # mutate(duration_seconds = as.numeric(max(time) - min(time),
        #                                      units = 'secs')) %>%
        # ungroup() %>%
        # filter(duration_seconds > 5, duration_seconds < 1800) %>%
        # select(-duration_seconds) %>%
        # Recompute sequence of grid cells after filtering
        # mutate(grid_uid = get_grid_uid(longitude, latitude)) %>%
        # group_by(grid_uid) %>%
        summarize(
          across(everything(), median, na.rm = T),
          across(where(is.numeric) & !ends_with(c('_ex', '_base')),
                 n = function(x) sum(!is.na(x)))
        ) %>%
        ungroup() %>%
        select(-grid_uid)
    }) %>%
    ungroup() %>%
    # Reorder columns and rows in output
    select(order(colnames(.))) %>%
    relocate(system_id, time) %>%
    arrange(system_id, time)
}


observations <- get_observations()
grid <- grid_observations(observations)
saveRDS(grid, 'data/gridded.rds')


# Add STILT receptor conventions
# gridded %>%
#   mutate(run_time = as.POSIXct(trunc(time, units = 'mins')),
#          lati = round(latitude, 8),
#          long = round(longitude, 8),
#          zagl = 3,
#          job_id = paste(strftime(run_time, '%Y%m%d%H%M', 'UTC'),
#                         long,
#                         lati,
#                         zagl,
#                         sep = '_')) %>%
#   saveRDS('data/receptors.rds')


grid_average <- grid %>%
  group_by(longitude, latitude) %>%
  summarize(across(where(is.numeric), mean, na.rm = T))


tracer_name <- 'co2d_ppm_ex'
color_range <- range(grid_average[[tracer_name]], na.rm = T)
colors <- colorNumeric('Spectral', rev = T, na.color = 'transparent',
                       domain = color_range)
map <- leaflet() %>%
  addProviderTiles('CartoDB.Positron') %>%
  # Replace CartoDB.Positron (white basemap) with Esri.WorldImagery (satellite)
  # addProviderTiles('Esri.WorldImagery') %>%
  addCircles(lng = grid_average$longitude,
             lat = grid_average$latitude,
             radius = 20,
             fillOpacity = 0.5,
             stroke = F,
             color = colors,
             popup = paste0(tracer_name, ': ', grid_average[[tracer_name]])) %>%
  addLegend(position = 'bottomright', 
            pal = colors,
            values = color_range)

# Use different pandoc to avoid file not found errors from CHPC's system install
Sys.setenv('RSTUDIO_PANDOC' = '/uufs/chpc.utah.edu/sys/installdir/pandoc/2.7.1/bin/')
saveWidget(map, paste0(tracer_name, '.html'))