Skip to content

Instantly share code, notes, and snippets.

@benfasoli
Created March 26, 2021 15:34
Show Gist options
  • Save benfasoli/b31994ba5084851f1d448feb2166168a to your computer and use it in GitHub Desktop.
Save benfasoli/b31994ba5084851f1d448feb2166168a to your computer and use it in GitHub Desktop.
Make interactive map by spatially averaging high resolution observations
#!/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'))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment