require(tidyverse); require(magrittr);
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.0 ✔ readr 2.1.6
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: magrittr
##
##
## Attaching package: 'magrittr'
##
##
## The following object is masked from 'package:purrr':
##
## set_names
##
##
## The following object is masked from 'package:tidyr':
##
## extract
require(sp);require(sf);require(rgdal);
## Loading required package: sp
## Loading required package: sf
## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
## Loading required package: rgdal
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'rgdal'
require(classInt);require(RColorBrewer);
## Loading required package: classInt
## Loading required package: RColorBrewer
require(ggplot2);require(ggmap);require(leaflet);
## Loading required package: ggmap
## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
## Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service>
## OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
## Attaching package: 'ggmap'
##
## The following object is masked from 'package:magrittr':
##
## inset
##
## Loading required package: leaflet
require(ggrepel); require(ggsn);
## Loading required package: ggrepel
## Loading required package: ggsn
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'ggsn'
library(ggpubr); library(devtools); library(patchwork); library(htmltools);
## Loading required package: usethis
library(prettymapr)
##
## Attaching package: 'prettymapr'
##
## The following objects are masked from 'package:ggmap':
##
## clear_geocode_cache, geocode
require(htmlwidgets); require(mapview)
## Loading required package: htmlwidgets
## Loading required package: mapview
Task#1: Plot at least two high-quality static maps at the zip code
level, with one using the COVID-19 testing data and one using a related
factor such as the number of nursing homes or the elderly population.
You can use either plot method for sf or ggplot method.
library(RColorBrewer)
pal <- brewer.pal(9, "OrRd")
class(pal)
## [1] "character"
acs_pop_by_zip <- st_read("/Users/areeba/Desktop/R/R-Spatial/Week9/acsPopByZip.gpkg")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source
## `/Users/areeba/Desktop/R/R-spatial/Week9/acsPopByZip.gpkg'
## using driver `GPKG'
## Simple feature collection with 180 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 913129 ymin: 120020.9 xmax: 1067113 ymax: 272710.9
## Projected CRS: NAD83 / New York Long Island (ftUS)
names(acs_pop_by_zip)
## [1] "ZIPCODE" "PO_NAME" "POPULAT" "COUNTY" "Positiv" "Total" "totPop"
## [8] "malPctg" "asianPp" "blackPp" "hspncPp" "whitePp" "eldrlyP" "geom"
acs_pop_by_zip <- st_transform(acs_pop_by_zip, 4326)
plot(acs_pop_by_zip["Positiv"],
main = "Positive Covid Cases",
breaks = "quantile", nbreaks = 9,
pal = pal,
graticule = st_crs(4326),
axes=TRUE,
reset=FALSE, )

plot(acs_pop_by_zip["eldrlyP"],
main = "Elderly Population",
breaks = "quantile", nbreaks = 9,
pal = pal,
graticule = st_crs(4326),
axes=TRUE,
reset=FALSE )
addnortharrow(pos = "topright", padin = c(0.1, 0.1), scale = 0.5)
addscalebar(pos = "topleft")
## Autodetect projection: assuming lat/lon (epsg 4326)
## Task#2: Use ggplot2 and other ggplot-compatible packages to create a
multi-map figure illustrating the possible relationship between COVID-19
confirmed cases or rate and another factor (e.g., the number of nursing
homes, neighborhood racial composition, elderly population, etc.).
p_legend <- c(1,2,3,4,5,6)
q_legend <- c(1,2,3,4,5,6)
brks_p <- classIntervals(acs_pop_by_zip$Positiv, n= 6, style = "quantile")
acs_pop_by_zip <- mutate(acs_pop_by_zip, Positive_cat = cut(Positiv, breaks = brks_p$brks, include.lowest = TRUE, labels = p_legend))
brks_e <- classIntervals(acs_pop_by_zip$eldrlyP, n= 6, style = "quantile")
## Warning in classIntervals(acs_pop_by_zip$eldrlyP, n = 6, style = "quantile"):
## var has missing values, omitted in finding classes
acs_pop_by_zip <- mutate(acs_pop_by_zip, eldrly_cat = cut(eldrlyP, breaks = brks_e$brks, include.lowest = TRUE, labels = q_legend))
g3 <- ggplot(acs_pop_by_zip) +
geom_sf(aes(fill=Positive_cat), show.legend = TRUE) +
scale_fill_brewer(palette = "Purples", name='Positive Covid Cases') +
coord_sf(xlim = c(-74.26, -73.70), ylim=c(40.49, 40.92), default_crs = sf::st_crs(4326) ) +
labs(x='Longitude', y='Latitude',
title='Positive Covid Cases',
caption = 'Data Source: [...]')+
theme(plot.title = element_text(size = 14, hjust = 0.5),
axis.text = element_text(size = 9),
legend.text = element_text(size = 8),
legend.position = "right",
legend.key.height = unit(0.6, "cm"),
legend.key.width = unit(0.4, "cm"),
panel.border = element_rect(color = "black", fill = NA))
g4 <- ggplot(acs_pop_by_zip) +
geom_sf(aes(fill=eldrly_cat), show.legend = TRUE) +
scale_fill_brewer(palette = "Blues", name='Elderly Population') +
coord_sf(xlim = c(-74.26, -73.70), ylim=c(40.49, 40.92), default_crs = sf::st_crs(4326) ) +
labs(x='Longitude', y='Latitude',
title='Elderly Population',
caption = 'Data Source: [...]')+
theme(
plot.title = element_text(size = 14, hjust = 0.5),
axis.text = element_text(size = 9),
legend.text = element_text(size = 8),
legend.position = "right",
legend.key.height = unit(0.6, "cm"),
legend.key.width = unit(0.4, "cm"),
panel.border = element_rect(color = "black", fill = NA))
g3 + g4 + plot_layout(ncol = 2, nrow = 1)
## Task#3: Create a web-based interactive map for COIVD-19 data using
tmap, mapview, or leaflet package and save it as a HTML file.
pal_pos <- colorQuantile("Purples", domain = acs_pop_by_zip$Positiv, n = 5)
pal_eld <- colorQuantile("Greens", domain = acs_pop_by_zip$eldrlyP, n = 5)
pal_tot <- colorQuantile("Blues", domain = acs_pop_by_zip$totPop, n = 5)
pos_lab <- paste0("ZIP: ", acs_pop_by_zip$ZIPCODE, "<br>Positive Cases: ", acs_pop_by_zip$Positiv)
eld_lab <- paste0("ZIP: ", acs_pop_by_zip$ZIPCODE, "<br>Elderly Population: ", acs_pop_by_zip$eldrlyP)
tot_lab <- paste0("ZIP: ", acs_pop_by_zip$ZIPCODE, "<br>Total Population: ", acs_pop_by_zip$totPop)
pos_popup <- paste0("<strong>ZIPCODE: </strong>", acs_pop_by_zip$ZIPCODE, "<br/>", "<strong>Positive Cases: </strong>", acs_pop_by_zip$Positiv, "<br/>", "<strong>Total Population: </strong>", acs_pop_by_zip$totPop, "<br/>")
eld_popup <- paste0("<strong>ZIPCODE: </strong>", acs_pop_by_zip$ZIPCODE, "<br/>", "<strong>Elderly Population: </strong>", acs_pop_by_zip$eldrlyP, "<br/>", "<strong>Total Population: </strong>", acs_pop_by_zip$totPop, "<br/>")
tot_popup <- paste0("<strong>ZIPCODE: </strong>", acs_pop_by_zip$ZIPCODE, "<br/>", "<strong>Positive Cases: </strong>", acs_pop_by_zip$Positiv, "<br/>",
"<strong>Elderly Population: </strong>", acs_pop_by_zip$eldrlyP, "<br/>", "<strong>Total Population: </strong>", acs_pop_by_zip$totPop, "<br/>")
acs_pop_by_zip <- st_transform(acs_pop_by_zip, 4326)
covid_elderly_map <- leaflet(acs_pop_by_zip) %>%
addTiles() %>%
addPolygons(
group = "Positive Cases",
stroke = TRUE,
fillColor = ~pal_pos(Positiv),
label = lapply(pos_lab, HTML),
popup = lapply(pos_popup, HTML),
color = "gray40",
weight = 0.5,
fillOpacity = 0.8, smoothFactor = 0.5) %>%
addPolygons(
group = "Elderly Population",
stroke = TRUE,
fillColor = ~pal_eld(eldrlyP),
label = lapply(eld_lab, HTML),
popup = lapply(eld_popup, HTML),
color = "gray40",
weight = 0.5,
fillOpacity = 0.8, smoothFactor = 0.5) %>%
addPolygons(
group = "Total Population",
stroke = TRUE,
fillColor = ~pal_tot(totPop),
label = lapply(tot_lab, HTML),
popup = lapply(tot_popup, HTML),
color = "gray40",
weight = 0.5,
fillOpacity = 0.8, smoothFactor = 0.5) %>%
addLegend(position = "bottomright", # location
group = "Total Population",
pal=pal_tot,
opacity = 0.8,
values=~totPop,
title = 'Total population') %>%
addLayersControl(overlayGroups = c("Positive Cases", "Elderly Population", "Total Population"), options = layersControlOptions(collapsed = FALSE)) %>%
hideGroup(c("Positive Cases", "Elderly Population"))
htmlwidgets::saveWidget(covid_elderly_map, "covid_elderly_total_map.html", selfcontained = TRUE)
covid_elderly_map