library(tidyverse)

The goal of this exercise is to visualize health code violations in New York City restaurants.

🚧 Following the instructions given in Section 18.6, Problem 1.

We need to limit the number of addresses we geocode because they take around 1 second each to process.

library(mdsr)
library(tidygeocoder)

violations_selected <- Violations %>%
  count(camis, dba, boro, building, street, zipcode) %>%
  filter(n > 90)  # n > 70 worked in spring 2022, but it's blowing out in spring 2023 (during the day).

violations_selected
## # A tibble: 8 × 7
##      camis dba                         boro      building street   zipcode     n
##      <int> <chr>                       <chr>        <int> <chr>      <int> <int>
## 1 41197758 CURRY LEAVES RESTAURANT     QUEENS       13531 40 ROAD    11354    92
## 2 41289382 JADE ASIAN RESTAURANT       QUEENS       13628 39 AVEN…   11354    91
## 3 41320205 DAI WAH YUMMY CITY          BROOKLYN      7218 18 AVEN…   11204    96
## 4 41510846 218 RESTAURANT              MANHATTAN   218220 GRAND S…   10013    94
## 5 41543033 NEW STAR SEAFOOD RESTAURANT BROOKLYN      1217 AVENUE U   11229    94
## 6 41586091 FENG CHENG YUAN RESTAURANT  MANHATTAN      100 BOWERY …   10013    94
## 7 41602559 RED CHOPSTICK               QUEENS       13617 41 AVEN…   11355    96
## 8 41690948 TIO POLLO RESTAURANT        QUEENS        1921 MOTT AV…   11691    93

Here we compute the geo-coordinates, filtering out an outlier that’s actually in Rochester, NY.

library(sf)
## Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
violations_geom <- violations_selected %>%
  mutate(address = paste(building, ", ", street, ", ", boro, ", NY,", zipcode)) %>%
  tidygeocoder::geocode(address, method = "osm") %>%
  na.omit() %>%
  filter(lat < 41.0) %>%
  st_as_sf(coords = c("long", "lat")) %>%
  st_set_crs(4326)  
## Passing 8 addresses to the Nominatim single address geocoder
## Query completed in: 8.9 seconds
violations_geom
## Simple feature collection with 3 features and 8 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -73.99915 ymin: 40.59877 xmax: -73.75442 ymax: 40.72011
## Geodetic CRS:  WGS 84
## # A tibble: 3 × 9
##      camis dba                       boro  building street zipcode     n address
## *    <int> <chr>                     <chr>    <int> <chr>    <int> <int> <chr>  
## 1 41510846 218 RESTAURANT            MANH…   218220 GRAND…   10013    94 218220…
## 2 41543033 NEW STAR SEAFOOD RESTAUR… BROO…     1217 AVENU…   11229    94 1217 ,…
## 3 41690948 TIO POLLO RESTAURANT      QUEE…     1921 MOTT …   11691    93 1921 ,…
## # ℹ 1 more variable: geometry <POINT [°]>
library(ggspatial)

violations_geom %>%
  ggplot() +
  aes(size = n) +
  annotation_map_tile(zoomin = 0, type = "osm") +
  geom_sf()
## Loading required namespace: raster
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
## Path to GDAL shared files: C:/Users/kvlinden/AppData/Local/R/win-library/4.3/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
## Path to PROJ shared files: C:/Users/kvlinden/AppData/Local/R/win-library/4.3/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.0-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Zoom: 12
## Fetching 12 missing tiles
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |======================================================================| 100%
## ...complete!

Now, we put them on an interactive map.

library(leaflet)

violations_geom %>%
  mutate(popup = paste(dba, address, paste("violations:", n), sep = "<br/>")) %>%
  leaflet() %>%
  addTiles() %>%
  addMarkers(popup = ~popup)