GTECH38520_2024Spring_week8_24075940

Explanation of the template

Update the title with your information. For example, it could be “GTECH_78520_2023S_Week12_JohnDoe25” for John Doe, 2022 Spring semester, week 12.

Also update the author name and date accordingly.

R Spaital Lab Assignment # 1

Reading in data:

acscovid <- st_read('/Users/mariyatrofymova/Desktop/Rspatial Study/R-Spatial_III_Lab/acsPopByZip.gpkg', layer="NYC_COVID_Pop_by_Zip")

## Reading layer `NYC_COVID_Pop_by_Zip' from data source 
##   `/Users/mariyatrofymova/Desktop/Rspatial Study/R-Spatial_III_Lab/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)

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.

display.brewer.all()

colorpal <- brewer.pal(7, "YlOrRd")

plot(acscovid["Positiv"], 
     main = "Covid Positives Per NYC Zipcode", 
     breaks = "quantile", nbreaks= 7,
     border = NA,
     pal = colorpal,
     graticule = st_crs(4326),
     axes = TRUE)

elderpal <- brewer.pal(7, "BrBG")

plot(acscovid["eldrlyP"], 
     main = "Elderly Populations Per NYC Zipcode", 
     breaks = "quantile", nbreaks = 7,
     border = NA,
     pal = elderpal,
     graticule = st_crs(4326),
     axes=TRUE,
     reset=FALSE)

### 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.

acscovid <- acscovid %>%
  mutate(Positiv = as.numeric(Positiv),
         totPop = as.numeric(totPop))

acscovid_long <- acscovid %>% 
  pivot_longer(cols = c(Positiv, totPop), names_to = "variables", values_to = "values") %>%
  st_as_sf()

breaks_lmh <- classIntervals(acscovid_long$values, n = 3, style = "quantile", na.rm = TRUE)$brks

## Warning in classIntervals(acscovid_long$values, n = 3, style = "quantile", :
## var has missing values, omitted in finding classes

breaks_lmh[1] <- breaks_lmh[1] - 0.0001

acscovid_long %<>% mutate(ValuesLevels = cut(values, 
                                             breaks_lmh, 
                                             labels = c('Low', 'Medium', 'High')))

ggplot(acscovid_long) + 
  geom_sf(data = acscovid_long %>% st_geometry(geom), col = rgb(0.6, 0.6, 0.98)) +
  geom_sf(aes(fill = ValuesLevels)) +
  scale_fill_brewer(palette = "OrRd", name = 'Levels') +
  labs(title = 'Covid Positive and Total Population Levels Per Zipcode',
       caption = 'American Community Survey') +
  facet_wrap(~variables) +
  coord_sf(datum = st_crs(4326)) +  
  theme(axis.text.x = element_blank(), 
        axis.text.y = element_blank(),
        legend.position = 'bottom')

###Create a web-based interactive map for COIVD-19 data using tmap, mapview, or leaflet package and save it as a HTML file.

library(leaflet) 

st_crs(acscovid)

## Coordinate Reference System:
##   User input: NAD83 / New York Long Island (ftUS) 
##   wkt:
## PROJCRS["NAD83 / New York Long Island (ftUS)",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 New York Long Island zone (US Survey feet)",
##         METHOD["Lambert Conic Conformal (2SP)",
##             ID["EPSG",9802]],
##         PARAMETER["Latitude of false origin",40.1666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-74,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",41.0333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",40.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["United States (USA) - New York - counties of Bronx; Kings; Nassau; New York; Queens; Richmond; Suffolk."],
##         BBOX[40.47,-74.26,41.3,-71.8]],
##     ID["EPSG",2263]]

acscovid <- st_transform(acscovid, 4326)

leaflet(acscovid) %>%
  addPolygons()

pal_cov <- colorQuantile("Oranges", NULL, n = 5)

pal_totpop <- colorQuantile("YloRed", 
                        domain = c(min(acscovid$totPop, na.rm = T), max(acscovid$totPop, na.rm = T)),
                        alpha = TRUE)

covid_popup <- paste0("<strong>Zip Code: </strong>", 
                  acscovid$ZIPCODE,
                   " <br/>",
                  "<strong> Place Name: </strong>",
                  acscovid$PO_NAME,
                  " <br/>",
                  "<strong> Confirmed COVID-19 Positives: </strong>",
                  acscovid$Positiv,
                  sep="")

highlights <- leaflet::highlightOptions(opacity = 1.0, fillColor = 'yellow')
labels <- leaflet::labelOptions(opacity = 1.0, textsize = '11px')

htmlMap <- leaflet(acscovid %>% sf::st_transform(4326)) %>%
        addPolygons(
                stroke = TRUE, 
                fillOpacity = 0.8,
                fillColor = ~pal_totpop(totPop),
                color = 'black', 
                weight = 1, 
                highlightOptions = highlights,
                label = paste(acscovid$totPop),
                group = "Total Population") %>%
        addPolygons(
                stroke = FALSE,
                fillColor = ~pal_cov(Positiv),
                fillOpacity = 0.8, smoothFactor = 0.5,
                popup = covid_popup,
                group = "NYC-COVID",
                label = paste0(acscovid$Positiv, " at ", acscovid$PO_NAME, ", ", acscovid$ZIPCODE),
                highlightOptions = highlights,
                labelOptions = labels) %>%
        addTiles(group = "OSM") %>%
        addProviderTiles("CartoDB.DarkMatter", group = "Carto") %>%
        addLayersControl(baseGroups = c("OSM", "Carto"), 
                         overlayGroups = c("NYC Covid Positives", "Total Population"))  

htmlMap

htmlwidgets::saveWidget(htmlMap, 'nyc_covid_pop_leaflet.html', 
                        selfcontained = FALSE, libdir = 'widget-lib')