Chlorophyll interpolation
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
library(tidyverse)
library(lubridate)
library(sf)
library(gstat)
library(sp)
library(raster)
library(mapview)
library(leaflet)
library(stars)
library(leafem)
data(rswqdat)
data(rsstatloc)
data(segmask)
prj <- "+proj=longlat +datum=WGS84 +no_defs"
mptyps <- c("CartoDB.Positron", "CartoDB.DarkMatter", "OpenStreetMap", "Esri.WorldImagery", "OpenTopoMap")
# non-bay stations
nonbay <- c('BH01', 'P Port 2', 'P Port 3', 'PM Out', '20120409-01', 'PPC41', 'P Port 4', 'PMB01')
# station locations
locs <- rsstatloc %>%
dplyr::select(station)
# current chl data
chldat <- rswqdat %>%
filter(!station %in% nonbay) %>%
# filter(source %in% 'fldep') %>%
filter(var %in% 'chla') %>%
dplyr::select(station, date, var, val) %>%
inner_join(locs, ., by = 'station') %>%
mutate(
date = floor_date(date, unit = 'week')
)
# color palette
colfun <- colorBin(
RColorBrewer::brewer.pal(9, 'Blues'),
domain = c(0, max(chldat$val, na.rm = T)),
bins = 8,
na.color = "#FFFFFF00" ,
alpha = FALSE,
reverse = FALSE
)
# chldat bbox
bbox <- chldat %>%
st_bbox %>%
st_as_sfc() %>%
st_as_sf() %>%
as_Spatial()
# empty grid for interpolation
grd <- as.data.frame(spsample(bbox, "regular", n = 100000))
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE # Create SpatialPixel object
fullgrid(grd) <- TRUE # Create SpatialGrid object
proj4string(grd) <- prj
# interpolate by date
mps <- chldat %>%
group_by(date) %>%
nest %>%
mutate(
mp = purrr::map(data, function(x){
wk <- x %>%
mutate(
col = colfun(val)
)
wkspa <- wk %>%
as_Spatial
# add crs to grid
proj4string(wkspa) <- prj
# interpolate
interp <- gstat::idw(log(val) ~ 1, wkspa, newdata = grd, idp = 2.0) %>%
raster %>%
mask(segmask) %>%
st_as_stars %>%
mutate(
var1.pred = exp(var1.pred)
)
out <- leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addStarsImage(interp, colors = colfun) %>%
addCircleMarkers(
data = wk,
stroke = TRUE,
color = 'black',
fill = TRUE,
fillColor = wk$col,
weight = 0.5,
fillOpacity = 1,
radius= 3,
label = ~val
) %>%
addLegend("topright", title = 'Chla (ug/L)', opacity = 1, values = chldat$val, pal = colfun) %>%
htmltools::tagList() %>%
as.character
return(out)
})
)
dts <- unique(chldat$date) %>% sort
for(i in seq_along(dts)){
mp <- mps %>%
filter(date %in% dts[i]) %>%
pull(mp) %>%
.[[1]]
cat(paste('\n## Week of', dts[i], '\n\n\n\n'))
cat(mp)
}