libraries

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(leaflet)

load data

month1 <- read.csv('nj_filtered_monthly/filtered_1.csv')
month2 <- read.csv('nj_filtered_monthly/filtered_2.csv')
month3 <- read.csv('nj_filtered_monthly/filtered_3.csv')
month4 <- read.csv('nj_filtered_monthly/filtered_4.csv')
month5 <- read.csv('nj_filtered_monthly/filtered_5.csv')
month6 <- read.csv('nj_filtered_monthly/filtered_6.csv')
month7 <- read.csv('nj_filtered_monthly/filtered_7.csv')
month8 <- read.csv('nj_filtered_monthly/filtered_8.csv')
month9 <- read.csv('nj_filtered_monthly/filtered_9.csv')
month10 <- read.csv('nj_filtered_monthly/filtered_10.csv')
month11 <- read.csv('nj_filtered_monthly/filtered_11.csv')
month12 <- read.csv('nj_filtered_monthly/filtered_12.csv')

full2022data <- rbind(month1, month2, month3, month4, month5, month6, month7, month8, month9, month10, month11, month12)

sf_nj <- st_read('Downloads/tl_2022_34_bg/tl_2022_34_bg.shp')
## Reading layer `tl_2022_34_bg' from data source 
##   `/Users/erandoni/Downloads/tl_2022_34_bg/tl_2022_34_bg.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6599 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.56359 ymin: 38.78866 xmax: -73.88506 ymax: 41.35761
## Geodetic CRS:  NAD83

aggregate

block_summary <- full2022data %>%
  group_by(AREA) %>%
  summarise(months_represented = n(), total_stops = sum(RAW_STOP_COUNTS))

block_summary$avg_monthly_stops <- block_summary$total_stops / block_summary$months_represented

block_summary <- merge(block_summary, sf_nj, by.x = 'AREA', by.y = 'GEOID')

block_summary <- st_as_sf(block_summary)

block_summary$log_stops <- log1p(block_summary$avg_monthly_stops)

pal <- colorNumeric(palette = colorRampPalette(c('green', 'yellow', 'red'))(100),
                    domain = block_summary$log_stops,
                    na.color = "#cccccc"
                    )

leaflet(data = block_summary) %>%
  addTiles() %>%
  addPolygons(
    color = 'black',
    weight = 1,
    label = ~paste('GEOID:', AREA, 'Average Monthly Stops', avg_monthly_stops),
    fillColor = ~pal(log_stops),
    fillOpacity = 0.7
  ) %>%
  addLegend(
    pal = pal,
    values = ~log1p(avg_monthly_stops),
    labFormat = labelFormat(
      transform = function(x) round(expm1(x))
    ),
    title = "Average Monthly Stops by Census Block")
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'