#' Get zones geometry from the database or built-in files
#'
#' @param zones_shp_name chr. The name of the zones_shp object in \code{data} to use if
#' the database does not contain the WKT geometry for zones. Possible values are
#' "zones_shp", "zones_shp_swim25".
#' @param db class \code{src}. The connection to the database.
#' @param tbl_name chr. The name of the table in \code{db} containing the geometry
#' Defaults to ALLZONES.
#' @param wkt_col chr. The name of the column in \code{tbl_name} containing the
#' WKT geometry. Defaults to WKTSURFACE
#' @param proj4string chr. The proj4 string for the
#' \code{wkt_col} coordinates. Defaults to Oregon Lambert.
#' @details This function extracts the geography for the zones from the scenario
#' database. It automatically reprojects the data to WGS84 using \code{st_transform}.
#' If the database does not contain geography, then data(zones_shp) is returned instead,
#' unless otherwise specified by the user.
#'
#' @return a \code{SpatialPolygonsDataFrame} in WGS84 lat long coordinates.
#'
#' @export
extract_zones <- function(zones_shp_name = NULL,
db = NULL,
tbl_name='ALLZONES',
wkt_col='WKTSURFACE',
proj4string = "+proj=lcc +lat_1=43 +lat_2=45.5 +lat_0=41.75 +lon_0=-120.5 +x_0=399999.9999999999 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=ft +no_defs"
){
wgs84 <- "+proj=longlat +datum=WGS84 +no_defs"
if ( !is.null(zones_shp_name) ){
get(zones_shp_name)
} else if ( !wkt_col %in% colnames(tbl(db, tbl_name)) ){
warning('Column ', wkt_col, ' not found in ', tbl_name, ' in ', db$con@dbname,
'. Using "zones_shp" instead.')
get('zones_shp')
} else {
zones_tbl <- dplyr::tbl(db, tbl_name) %>%
dplyr::collect()
zones_tbl <- zones_tbl %>%
#dplyr::filter(!is.na(!!quo(wkt_col))) # Doesn't work
dplyr::filter(!is.na(.[[wkt_col]]) & !grepl('EMPTY', .[[wkt_col]])) %>%
dplyr::mutate(id = as.character(row_number() - 1)) %>%
dplyr::left_join(y = regions, by='COUNTY')
# Make sure AZONE is capitalized because other functions assume that it is.
names(zones_tbl)[grep('Azone', names(zones_tbl))] <- 'AZONE'
zones_sf <- sf::st_as_sf(as.data.frame(zones_tbl), wkt=wkt_col, crs=proj4string) %>%
sf::st_transform(crs=wgs84)
zones_shp <- as(zones_sf, "Spatial")
}
return(zones_shp)
}
#' Leaflet plot of change over time in one scenario.
#'
#' @param db The scenario database
#' @param year1 The first year
#' @param year2 The second year
#' @param ... Additional parameters passed to \code{\link{extract_zones}}.
#' Specify the source for the zones here.
#' @details This function creates an interactive Leaflet plot of the alpha zone
#' system, showing annualized implied growth rates of key socioeconomic
#' variables between two years. For instance, the user can select 2010 and
#' 2030 as analysis years, and the plot will show the implied average
#' exponential growth rate between them.
#'
#' As a note, SWIM does not generate alpha-zone level statistics in years when
#' the transport model does not run. If the user asks for year that does not
#' exist, then the plot will use the nearest available year.
#'
#'
#' @export
change_leaflet <- function(db, year1 = 2010, year2 = 2030, ...){
# Get nearest years in db to year1 and year2
#
# The database doesn't have population and employment at
# alpha zone level in years where the transport model doesn't run.
# This simply corrects the value of year1 and year 2 if the user asks for
# a year that doesn't exist.
years <- 1990 + as.numeric(names(table(
dplyr::tbl(db, "AZONE") %>% dplyr::select(TSTEP) %>% dplyr::collect(n=Inf) %>% .$TSTEP)))
year1 <- years[which(abs(years - year1) == min(abs(years - year1)))]
year2 <- years[which(abs(years - year2) == min(abs(years - year2)))]
# Get scenario information and put it onto the shapefile for
# leaflet plotting.
shp <- extract_zones(db=db, ...)
shp@data <- shp@data %>%
dplyr::left_join(extract_zonedata(db, year1, year2))
# Create a palette for
palq <- leaflet::colorFactor(
palette = "RdBu",
domain = cut_diverror(shp@data$pop_rate)
)
zone_leaflet(shp) %>%
leaflet::addPolygons(
group = "Population", fill = TRUE, color = FALSE, stroke=FALSE,
fillColor = ~palq(cut_diverror(pop_rate)),
popup = change_popup(shp, "pop", year1, year2)
) %>%
leaflet::addPolygons(
group = "Employment", fill = TRUE, color = FALSE, stroke=FALSE,
fillColor = ~palq(cut_diverror(emp_rate)),
popup = change_popup(shp, "emp", year1, year2)
) %>%
leaflet::addPolygons(
group = "HH", fill = TRUE, color = FALSE, stroke=FALSE,
fillColor = ~palq(cut_diverror(hh_rate)),
popup = change_popup(shp, "hh", year1, year2)
) %>%
leaflet::addLayersControl(
overlayGroups = c("Population", "Employment", "HH"),
options = leaflet::layersControlOptions(collapsed = FALSE)
)
}
#' Leaflet plot of difference between two scenarios.
#'
#' @param db1 Reference scenario VIZ database connection.
#' @param db2 Current scenario VIZ database connection.
#' @param year Year in which to compare the two scenarios.
#' @param variable Which variable to show in difference plot; currently
#' \code{variable = c("Population", "Employment", "HH")}
#' @param scen_names Names of the scenarios in the comparison. Defaults to
#' "Reference", "Current".
#' @param ... Additional parameters passed to \code{\link{extract_zones}}.
# Specify the source for the zones here.
#' @details Note that if the zones geography is pulled from the database (the default),
#' only zones from \code{db1} are used for the map. In case zones differ
#' between \code{db1} and \code{db2}, note that the difference is based
#' on a join between the zone data based on \code{AZONE}.
#' @export
diff_leaflet <- function(db1, db2, year,
variable = c("Population", "Employment", "HH"),
scen_names = c("Reference", "Current"), ...){
if(variable == "Population"){
crit <- "variable == 'pop'"
} else if(variable == "Employment"){
crit <- "variable == 'emp'"
} else {
crit <- "variable == 'hh'"
}
# Get nearest years in db to year
#
# The database doesn't have population and employment at
# alpha zone level in years where the transport model doesn't run.
# This simply corrects the value of year if the user asks for
# a year that doesn't exist.
years <- 1990 + as.numeric(names(table(
dplyr::tbl(db1, "AZONE") %>% dplyr::select(TSTEP) %>% dplyr::collect(n=Inf) %>% .$TSTEP)))
year <- years[which(abs(years - year) == min(abs(years - year)))]
# get se data from both scenarios and join together ----
se <- extract_zonedata(db1, year) %>%
dplyr::left_join(extract_zonedata(db2, year), by = "AZONE") %>%
tidyr::gather(var, value, -AZONE) %>%
tidyr::separate(var, c("variable", "year", "sc")) %>%
dplyr::select(-year) %>%
tidyr::spread(sc, value) %>%
dplyr::filter_(.dots = crit) %>%
# calculate absolute difference between scenarios
dplyr::mutate(
diff = y - x,
pct = diff / x * 100
)
# Get scenario information and put it onto the shapefile for
# leaflet plotting.
shp <- extract_zones(db=db1, ...)
shp@data <- shp@data %>%
dplyr::left_join(se, by = "AZONE")
# Create a palette for
palq <- leaflet::colorFactor(
palette = "PRGn",
domain = cut_abserror(shp@data$diff)
)
# Fix legend as shown here:
# https://github.com/rstudio/leaflet/issues/477#issuecomment-678542769
zone_leaflet(shp) %>%
leaflet::addPolygons(
group = "Absolute", stroke = FALSE, fillOpacity=0.7,
color = ~palq(cut_abserror(diff)),
popup = diff_popup(shp, variable, scen_names)
) %>%
leaflet::addLegend(
"bottomright",
pal = palq,
values = ~cut_abserror(diff),
group = "Absolute",
className = "info legend Absolute",
title = paste0("Change in ", variable, " per sqmi<br>",
scen_names[1], " - ", scen_names[2] )
) %>%
leaflet::addPolygons(
group = "Percent", stroke = FALSE, fillOpacity=0.7,
color = ~palq(cut_abserror(pct)),
popup = diff_popup(shp, variable, scen_names)
) %>%
leaflet::addLegend(
"bottomright",
pal=palq,
values = ~cut_abserror(pct),
group = "Percent",
className = "info legend Percent",
title = paste0("Pct change in ", variable, " per sqmi<br>",
scen_names[1], " - ", scen_names[2] )
) %>%
leaflet::addLayersControl(
baseGroups = c("Absolute", "Percent"),
options = leaflet::layersControlOptions(collapsed = FALSE)
) %>%
htmlwidgets::onRender("
function(el, x) {
var updateLegend = function () {
var selectedGroup = document.querySelectorAll('input:checked')[0].nextSibling.innerText.substr(1);
document.querySelectorAll('.legend').forEach(a => a.hidden=true);
document.querySelectorAll('.legend').forEach(l => {
if (l.classList.contains(selectedGroup)) l.hidden=false;
});
};
updateLegend();
this.on('baselayerchange', el => updateLegend());
}"
)
}
#' Extract se data for leaflet zone plots
#'
#' @inheritParams change_leaflet
#'
#' @details If \code{year2 = NULL}, then will return only the values from year 1
#' and will not calculate the implied exponential growth rate.
#'
extract_zonedata <- function(db, year1, year2 = NULL){
if(is.null(year2)){
year2 <- year1
}
# pre-construct growth rate function call in order to mix variables
# and field names. See vignettes("nse") for details.
grt_exp <- lazyeval::interp(
~ calc_exprate(p1, p2, year1, year2) * 100,
p1 = as.name(year1), p2 = as.name(year2))
# Get socioeconomic data from database.
se <- dplyr::tbl(db, "AZONE") %>%
dplyr::transmute(AZONE, pop = POPULATION, emp = EMPLOYMENT,
hh = TOTALHHS, year = TSTEP + 1990) %>%
dplyr::filter(year %in% c(year1, year2)) %>%
# return density of variable
dplyr::left_join(
dplyr::tbl(db, "ALLZONES") %>%
dplyr::transmute(AZONE = Azone, sqmi = AREASQFT * 3.587048e-8),
by = "AZONE") %>%
dplyr::collect(n=Inf) %>%
dplyr::mutate_at(vars(pop:hh), funs( ./sqmi)) %>%
dplyr::select(-sqmi) %>%
tidyr::gather(variable, value, -AZONE, -year) %>%
tidyr::spread(year, value)
if(year2 != year1){
se <- se %>%
# calculate implied growth rate
dplyr::mutate_("rate" = grt_exp)
}
# reformat
se %>%
tidyr::gather(period, value, -AZONE, -variable) %>%
unite(var, c(variable, period)) %>%
tidyr::spread(var, value, fill = NA)
}
#' Construct a base leaflet zone plot
#'
#' This function creates a standard leaflet base image that we can add
#' information to.
#'
#' @param shp An object of class \code{SpatialPolygonsDataFrame} with the
#' plotting fields already appended.
#'
#'
zone_leaflet <- function(shp){
leaflet::leaflet(shp) %>%
leaflet::addProviderTiles("CartoDB.Positron") %>%
leaflet::addPolygons(color = "grey", weight = 0.5, fill = FALSE)
}
#' Build change popup tag for leaflet
#'
#' @inheritParams zone_leaflet
#' @inheritParams change_leaflet
#' @param var Variable to create popup tag for (e.g., \code{"pop"})
#'
change_popup <- function(shp, var, year1, year2){
# build popup
zone_info <- paste0("<strong>Alpha Zone: </strong>", shp@data$AZONE)
var_info <- paste0(
"<strong>", year1, " ", var, " density: </strong>",
round(shp@data[, paste0("", var, "_", year1)], digits = 2), "<br>",
"<strong>", year2, " ", var, " density: </strong>",
round(shp@data[, paste0("", var, "_", year2)], digits = 2), "<br>",
"<strong>Growth rate </strong>",
round(shp@data[, paste0(var, "_rate")], digits = 3), "%"
)
paste(zone_info, var_info, sep = "</br>")
}
#' Build diff popup tag for leaflet
#'
#' @inheritParams zone_leaflet
#' @inheritParams diff_leaflet
#' @param var Variable to create popup tag for (e.g., \code{"pop"})
#'
diff_popup <- function(shp, var, scen_names){
# build popup
zone_info <- paste0("<strong>Alpha Zone: </strong>", shp@data$AZONE)
var_info <- paste0(
"<strong>", scen_names[1], " ", var, " per sqmi: </strong>",
round(shp@data[, "x"], digits = 3), "<br>",
"<strong>", scen_names[2], " ", var, " per sqmi: </strong>",
round(shp@data[, "y"], digits = 3), "<br>",
"<strong>Difference </strong>",
round(shp@data[, ("diff")], digits = 3), "<br>",
"<strong>Percent diff </strong>",
round(shp@data[, "pct"], digits = 3)
)
paste(zone_info, var_info, sep = "</br>")
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.