R Spaital Lab Assignment #3
# load in RData file from previous assignments with processed datasets as sf objects
load("C:/Users/jdhoc/Desktop/Hunter/SPRING 2024/GEOG 39300 R/R-spatial/Session_7/data/assignment1_data.RData")
load('C:/Users/jdhoc/Desktop/Hunter/SPRING 2024/GEOG 39300 R/R-spatial/Session_8/week8_data.RData')
# read in covid acs data
acs_covid <- sf::st_read('C:/Users/jdhoc/Desktop/Hunter/SPRING 2024/GEOG 39300 R/R-spatial/Session_9/acsPopByZip.gpkg', layer="NYC_COVID_Pop_by_Zip")
## Reading layer `NYC_COVID_Pop_by_Zip' from data source
## `C:\Users\jdhoc\Desktop\Hunter\SPRING 2024\GEOG 39300 R\R-spatial\Session_9\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)
task 1: Plot at least two high-quality static maps with one using
the COVID-19 data and one using a related factor. You can use either
plot method for sf or ggplot method.
# filter health facilities for nyc only
nyc_hf <- nys_hf %>% dplyr::filter(facility_county == "New York" |
facility_county == "Kings" | facility_county == "Queens" |
facility_county == "Bronx" | facility_county == "Richmond")
# get breaks for covid positive tests
breaks_qt <- classIntervals(c(min(acs_covid$Positiv) - .00001, acs_covid$Positiv), n = 6, style = "quantile")
acs_covid <- mutate(acs_covid, pos_cat = cut(Positiv, breaks_qt$brks, dig.lab = 4, digits=1))
# use ggrepel to disperse labels
labelCoords <- acs_covid %>% sf::st_centroid() %>%
dplyr::filter(Positiv > 1500) %>%
sf::st_coordinates();
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
labelData <- acs_covid %>% sf::st_centroid() %>%
dplyr::filter(Positiv > 1500) %>%
dplyr::mutate(x = labelCoords[,1], y=labelCoords[,2])
## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
# create map of NYC COVID Positive Cases by ZIP with health facilities
ggplot(acs_covid) +
geom_sf(aes(fill = pos_cat), color = "white") +
scale_fill_brewer(palette = "YlOrRd", name='Positive Tests') +
labs(x= 'Longitude', y= 'Latitude',
title= 'NYC COVID Positive Tests by ZIP 19 Apr 2020',
caption = 'Data Source: American Community Survey') +
geom_label_repel(data = labelData,
aes(x = x, y = y, label = Positiv),
label.size = .09,
size = 3,
segment.color = rgb(0.5,0.5,0.9),
segment.size = 0.8)
# get breaks for percent positive tests
breaks_qt <- classIntervals(c(min(acs_covid$Positiv / acs_covid$Total) - .00001, acs_covid$Positiv / acs_covid$Total), n = 4, style = "quantile")
acs_covid <- mutate(acs_covid,
perc_pos_cat = cut(Positiv / Total, breaks_qt$brks, dig.lab = 4, digits=1),
perc_pos = Positiv / Total)
# plot the % of positive test cases
ggplot(acs_covid) +
geom_sf(aes(fill = perc_pos_cat), color = "white") +
scale_fill_brewer(palette = "Purples", name='Percent Positive Tests') +
labs(x= 'Longitude', y= 'Latitude',
title= 'NYC COVID Positive Test Rate by ZIP 19 Apr 2020',
caption = 'Data Source: American Community Survey')
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, number of food stores, neighborhood racial composition, elderly
population, etc.). The maps should be put side by side on one single
page. Add graticule to at least one of those maps and label some of the
feature on the map where applicable and appropriate.
# calculate population density per sq mi
acs_covid$area <- st_area(acs_covid) %>% set_units(mi^2) %>% as.numeric()
acs_covid$pop_density <- acs_covid$POPULAT / acs_covid$area
# get breaks for density levels
breaks_lmh <- classIntervals(acs_covid$pop_density, n = 3, style = "quantile")
breaks_lmh$brks[1] %<>% magrittr::subtract(0.0001)
acs_covid %<>% mutate(density_lvl = cut(pop_density, breaks_lmh$brks, labels = c('Low','Medium','High')))
pop_dens_map <- ggplot(acs_covid) +
geom_sf(aes(fill = density_lvl), color = "white") +
scale_fill_brewer(palette = "Greens", name='Population Density') +
labs(title = "ZIP Population Density") +
theme(legend.position = 'bottom',
legend.direction = 'vertical')
# get breaks for covid positive levels
breaks_lmh <- classIntervals(acs_covid$Positiv, n = 3, style = "quantile")
breaks_lmh$brks[1] %<>% magrittr::subtract(0.0001)
acs_covid %<>% mutate(pos_lvl = cut(Positiv, breaks_lmh$brks, labels = c('Low','Medium','High')))
# edit covid pos map from previous task to fit with pop dens map next to it
covid_pos_map <- ggplot(acs_covid) +
geom_sf(aes(fill = pos_lvl), color = "white") +
scale_fill_brewer(palette = "YlOrRd", name='Positive Tests') +
labs(title = "Positive COVID Tests by ZIP") +
theme(legend.position = 'bottom',
legend.direction = 'vertical')
ggarrange(covid_pos_map, pop_dens_map)
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.
acs_covid <- st_transform(acs_covid, 4326)
pal_fun <- colorQuantile("Oranges", NULL, n = 4)
p_tip <- paste0("ZIP Code: ", acs_covid$ZIPCODE)
p_popup <- paste0("<strong>COVID Tests with Positive Result: </strong>",
(acs_covid$perc_pos * 100) %>% round(1), "%", "<br/>",
sep = "")
breaks_qt <- classIntervals(c(min(acs_covid$perc_pos) - .00001, acs_covid$perc_pos), n = 4, style = "quantile")
htmlMap <- leaflet() %>%
addTiles() %>%
addPolygons(
data = acs_covid,
stroke = FALSE,
fillColor = ~pal_fun(perc_pos), # set fill color with function from above and value
fillOpacity = 0.8, smoothFactor = 0.5, # make it nicer
popup = p_popup, # add popup
label = p_tip) %>%
addLegend("bottomright",
colors = brewer.pal(4, "Oranges"),
labels = paste0("up to ", format(breaks_qt$brks[-1], digits = 2)),
title = "NYC Percent of COVID Tests that are Positive")
htmlMap
# save as html file
htmlwidgets::saveWidget(htmlMap, 'nyc_percCovidPos_leaflet.html')