require(tidyverse); require(magrittr)
require(sp); require(sf); require(rgdal)
require(classInt); require(RColorBrewer)
require(ggplot2); require(ggmap); require(leaflet)
require(ggrepel); require(ggsn)
library(ggpubr); library(devtools)
require(htmlwidgets); require(mapview)
Coviddata <- sf::st_read('C:/Users/iiStr/Documents/R course/Week_09/acsPopByZip.gpkg',
layer = "NYC_COVID_Pop_by_Zip")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source
## `C:\Users\iiStr\Documents\R course\Week_09\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)
Nycdata <- Coviddata %>% sf::st_transform(2263)
nycwgs <- Coviddata %>% sf::st_transform(4326)
Nycdata <- Nycdata %>%
mutate(Positiv_cat = cut(Positiv,
breaks = classInt::classIntervals(Positiv, n = 5,
style = "quantile")$brks,
include.lowest = TRUE))
map1 <- ggplot(data = Nycdata) +
geom_sf(aes(fill = Positiv_cat), color = "black", linewidth = 0.2) +
scale_fill_brewer(palette ="OrRd",name ="COVID-19 Cases"
) +
coord_sf(
crs = sf::st_crs(Nycdata),
datum = sf::st_crs(4326),
expand = FALSE
) +
labs(
title = "COVID-19 Confirmed Cases by Zip Code",
x = "Longitude (°)",
y = "Latitude (°)"
)
print(map1)
ggsave("map1.png", plot = map1, width = 8, height = 9, dpi = 300)
Nycdata <- Nycdata %>%
mutate(eldrlyP_cat = cut(eldrlyP,
breaks = classInt::classIntervals(eldrlyP, n = 5,
style = "quantile", na.rm = T)$brks,
include.lowest = TRUE))
map2 <- ggplot(data = Nycdata) +
geom_sf(aes(fill = eldrlyP_cat), color = "black", linewidth = 0.2) +
scale_fill_brewer(palette ="YlOrRd",name ="Population\nAged 65+ (%)")+
coord_sf(
crs = sf::st_crs(Nycdata),
datum = sf::st_crs(4326),
expand = FALSE
) +
labs(
title = "Elderly Population by Zip Code",
x = "Longitude (°)",
y = "Latitude (°)"
)
print(map2)
ggsave("map2.png", plot = map2, width = 8, height = 9, dpi = 300)
Nycdata <- Nycdata %>%
mutate(Positiv_cat = cut(Positiv,
breaks = classInt::classIntervals(Positiv, n = 5,
style = "quantile")$brks,
include.lowest = TRUE))
Nycdata <- Nycdata %>%
mutate(eldrlyP_cat = cut(eldrlyP,
breaks = classInt::classIntervals(eldrlyP, n = 5,
style = "quantile", na.rm = T)$brks,
include.lowest = TRUE))
bbox <- sf::st_bbox(Nycdata)
pleft <- ggplot(data = Nycdata) +
geom_sf(aes(fill = Positiv_cat), color = "black", linewidth = 0.15) +
scale_fill_brewer(palette ="OrRd",name ="COVID-19 Cases"
) +
coord_sf(
xlim = c(bbox["xmin"], bbox["xmax"]),
ylim = c(bbox["ymin"], bbox["ymax"]),
crs = sf::st_crs(Nycdata),
datum = sf::st_crs(4326),
expand = FALSE
) +
labs(title = "A. COVID-19 Confirmed Cases",
x = "Longitude (°)", y = "Latitude (°)")
pright <- ggplot(data = Nycdata) +
geom_sf(aes(fill = eldrlyP_cat), color = "black", linewidth = 0.15) +
scale_fill_brewer(palette ="YlOrRd",name ="Population\nAged 65+ (%)")+
coord_sf(
xlim = c(bbox["xmin"], bbox["xmax"]),
ylim = c(bbox["ymin"], bbox["ymax"]),
crs = sf::st_crs(Nycdata),
datum = sf::st_crs(4326),
expand = FALSE
) +
labs(title = "B. Elderly Population (Age 65+)",
x = "Longitude (°)", y = "Latitude (°)")
combine <- ggarrange(pleft,pright, nrow = 1,ncol = 2)
combine
ggsave("combinedmap.png", plot = combine, width = 20, height = 12, dpi = 300)
polyHighlight <- leaflet::highlightOptions(
opacity = 1.0,
fillColor = "black",
fillOpacity = 0.3,
bringToFront = TRUE
)
labelOpts <- leaflet::labelOptions(opacity = 1.0, textsize = "11px")
pal_covid <- colorQuantile("YlOrRd", domain = nycwgs$Positiv, n = 5)
pal_elderly <- colorQuantile("YlGnBu", domain = nycwgs$eldrlyP, n = 5)
pal_white <- colorQuantile("Purples", domain = nycwgs$whitePp, n = 5)
popupcovid <- paste0(
"<strong>Zip Code: </strong>",nycwgs$ZIPCODE, "<br/>",
"<strong>Place: </strong>", nycwgs$PO_NAME, "<br/>",
"<strong>Confirmed Cases: </strong>", nycwgs$Positiv, "<br/>",
"<strong>Elderly Pop (%): </strong>", round(nycwgs$eldrlyP, 1), "%"
)
popupelderly <- paste0(
"<strong>Zip Code: </strong>",nycwgs$ZIPCODE, "<br/>",
"<strong>Place: </strong>",nycwgs$PO_NAME, "<br/>",
"<strong>Elderly Pop (%): </strong>",round(nycwgs$eldrlyP, 1), "%", "<br/>",
"<strong>Total Population: </strong>",format(nycwgs$totPop, big.mark = ",")
)
popupwhite <- paste0(
"<strong>Zip Code: </strong>",nycwgs$ZIPCODE, "<br/>",
"<strong>Place: </strong>",nycwgs$PO_NAME, "<br/>",
"<strong>White Pop (%): </strong>", round(nycwgs$whitePp, 1), "%", "<br/>",
"<strong>Confirmed Cases: </strong>",nycwgs$Positiv
)
popup <- paste0(popupcovid,
"<br/><hr/>",
popupelderly,
"<br/><hr/>",
popupwhite)
htmlMap <- leaflet(nycwgs) %>%
addTiles(group = "OSM (default)") %>%
addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
addProviderTiles("Esri.WorldImagery", group = "ESRI Imagery") %>%
addPolygons(
fillColor = ~pal_covid(Positiv),
fillOpacity = 0.75,
color = "grey",
weight = 0.8,
stroke = TRUE,
smoothFactor = 0.3,
popup = popup,
label = ~paste0(Positiv, " cases — ", PO_NAME, " (", ZIPCODE, ")"),
labelOptions = labelOpts,
highlightOptions = polyHighlight,
group = "COVID-19 Cases"
) %>%
addLegend(
position = "bottomright",
pal = pal_covid,
values = ~Positiv,
title = "COVID-19 Cases<br>(quantile)",
opacity = 0.8,
group = "COVID-19 Cases"
) %>%
addPolygons(
fillColor = ~pal_elderly(eldrlyP),
fillOpacity = 0.75,
color = "grey",
weight = 0.8,
stroke = TRUE,
smoothFactor = 0.3,
popup = popup,
label = ~paste0(round(eldrlyP, 1), "% elderly — ", PO_NAME),
labelOptions = labelOpts,
highlightOptions = polyHighlight,
group = "Elderly Population"
) %>%
addLegend(
position = "bottomright",
pal = pal_elderly,
values = ~eldrlyP,
title = "Elderly Pop. (%)<br>(quantile)",
opacity = 0.8,
group = "Elderly Population"
) %>%
addPolygons(
fillColor = ~pal_white(whitePp),
fillOpacity = 0.75,
color = "grey",
weight = 0.8,
stroke = TRUE,
smoothFactor = 0.3,
popup = popup,
label = ~paste0(round(whitePp, 1), "% white — ", PO_NAME),
labelOptions = labelOpts,
highlightOptions = polyHighlight,
group = "White Population"
) %>%
addLegend(
position = "bottomright",
pal = pal_white,
values = ~whitePp,
title = "White Pop. (%)<br>(quantile)",
opacity = 0.8,
group = "White Population"
) %>%
addLayersControl(
baseGroups = c("OSM (default)", "Carto", 'ESRI Imagery'),
overlayGroups = c("COVID-19 Cases", "Elderly Population", "White Population"),
options = layersControlOptions(collapsed = T)
) %>%
hideGroup("Elderly Population") %>%
hideGroup("White Population")
htmlMap