njmapping
Eran Doni
2024-08-20
import libraries
library(dplyr)
library(sf)
library(sp)
library(tidyr)
library(spdep)
library(stringr)
library(dbscan)
library(ggplot2)
library(h3jsr)
library(GGally)
library(leaflet)
library(cluster)
library(clusterSim)
library(clValid)
library(gstat)
library(raster)
library(viridis)
library(viridisLite)
library(tidycensus)load datasets
FARS_wide <- read.csv('FARS_NJ-NY-CT-18-22.csv')
SDOH_df <- read.csv('SDOH_NJ-NY-CT.csv')
combined_sf <- read.csv('combined_health_crash.csv')
sf_nj <- st_read('cb_2018_34_tract_500k/cb_2018_34_tract_500k.shp')## Reading layer `cb_2018_34_tract_500k' from data source
## `/Users/erandoni/mechkar/cb_2018_34_tract_500k/cb_2018_34_tract_500k.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 2006 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -75.55961 ymin: 38.92852 xmax: -73.89398 ymax: 41.35742
## Geodetic CRS: NAD83sf_ny <- st_read('cb_2018_36_tract_500k/cb_2018_36_tract_500k.shp')## Reading layer `cb_2018_36_tract_500k' from data source
## `/Users/erandoni/mechkar/cb_2018_36_tract_500k/cb_2018_36_tract_500k.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 4906 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -79.76215 ymin: 40.4961 xmax: -71.85621 ymax: 45.01585
## Geodetic CRS: NAD83sf_ct <- st_read('cb_2018_09_tract_500k/cb_2018_09_tract_500k.shp')## Reading layer `cb_2018_09_tract_500k' from data source
## `/Users/erandoni/mechkar/cb_2018_09_tract_500k/cb_2018_09_tract_500k.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 830 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.72777 ymin: 40.98014 xmax: -71.78699 ymax: 42.05059
## Geodetic CRS: NAD83tristate_sf <- rbind(sf_nj, sf_ny, sf_ct)
combined_sf$GEOID[combined_sf$STATEFP == 9] <- paste0(0, combined_sf$GEOID)
combined_sf <- merge(tristate_sf, combined_sf, by = 'GEOID', all.x = TRUE)health_quantiles <- quantile(combined_sf$health_intensity, probs = seq(0, 1, by = 0.2), na.rm=TRUE)
combined_sf$health_cat <- cut(combined_sf$health_intensity, breaks = health_quantiles, include.lowest = TRUE, labels=FALSE)#c('extremely high', 'high', 'moderate', 'low', 'extremely low'))
palette <- colorFactor(palette = colorRampPalette(c("green", "limegreen","orange", "red", "darkred"))(10), domain = combined_sf$health_cat)
nb <- combined_sf %>% filter(StateAbbr == "NJ" & CountyName == "Middlesex")
nb <- nb %>% filter((startsWith(GEOID, "34023005") | startsWith(GEOID, "34023009300") | startsWith(GEOID, "34023006002")) & (!startsWith(GEOID, "3402300500")))
leaflet(data = nb, options = leafletOptions(minZoom = 6)) %>%
addTiles() %>%
addPolygons(
fillColor = ~palette(health_cat),
fillOpacity = 0.7,
color = "black",
weight = 1,
opacity = 1,
highlightOptions = highlightOptions(
weight = 3,
color = "white",
bringToFront = TRUE
),
label = ~paste("Health Category:", health_cat)
) %>%
addLegend(
pal = palette,
values = ~health_cat,
title = "Health Category",
position = "bottomright"
) %>%
addControl(
html = "<h2>New Brunswick Health Category by Census Tract</h2>",
position = "topright",
className = "map-title"
)jc <- combined_sf %>% filter(COUNTYFP.x == "017")
jc <- jc %>% filter(startsWith(GEOID, "34017007") | startsWith(GEOID, "34017002") | startsWith(GEOID, "34017005") | startsWith(GEOID, "34017006") | startsWith(GEOID, "34017000") | startsWith(GEOID, "34017004") | startsWith(GEOID, "34017001") | startsWith(GEOID, "34017003") | startsWith(GEOID, "34017980100"))
leaflet(data = jc, options = leafletOptions(minZoom = 6)) %>%
addTiles() %>%
addPolygons(
fillColor = ~palette(health_cat),
fillOpacity = 0.7,
color = "black",
weight = 1,
opacity = 1,
highlightOptions = highlightOptions(
weight = 3,
color = "white",
bringToFront = TRUE
),
label = ~paste("Health Category:", health_cat)
) %>%
addLegend(
pal = palette,
values = ~health_cat,
title = "Health Category",
position = "bottomright"
) %>%
addControl(
html = "<h2>Jersey City Health Category by Census Tract</h2>",
position = "topright",
className = "map-title"
)jc_nb <- rbind(jc, nb)trainstations <- read.csv('/Users/erandoni/Downloads/NJTransit_Rail_Stations_-7781331129518987903.csv')
lrstations <- read.csv('/Users/erandoni/Downloads/NJTransit_Light_Rail_Stations_7467717369539091415.csv')
leaflet(data = jc_nb, options = leafletOptions(minZoom = 9)) %>%
addTiles() %>%
addPolygons(
fillColor = ~palette(health_cat),
fillOpacity = 0.7,
color = "black",
weight = 1,
opacity = 1,
label = ~paste("Health Category:", health_cat)
) %>%
addCircleMarkers(
data = trainstations,
~LONGITUDE, ~LATITUDE,
color = "darkblue",
fillColor = "darkblue",
fillOpacity = 1,
radius = 5,
popup = ~STATION_ID
) %>%
addCircleMarkers(
data = lrstations,
~LONGITUDE, ~LATITUDE,
color = "darkcyan",
fillColor = "darkcyan",
fillOpacity = 1,
radius = 5,
popup = ~STATION
) %>%
addLegend(
pal = palette,
values = ~health_cat,
title = "Health Category",
position = "bottomright"
) %>%
addControl(
html = "<h2>NB + JC Health Category by Tract, and Train/LR Stations</h2>",
position = "topright",
className = "map-title"
)blocks <- st_read('/Users/erandoni/Downloads/tl_2021_34_bg/tl_2021_34_bg.shp')## Reading layer `tl_2021_34_bg' from data source
## `/Users/erandoni/Downloads/tl_2021_34_bg/tl_2021_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: NAD83middlesex <- blocks %>% filter(COUNTYFP == "023")
ggplot(middlesex) +
geom_sf()
remote <- read.csv('/Users/erandoni/Downloads/remoteworkMiddlesex.csv')
middlesex <- merge(middlesex, remote, by = "GEOID", all.x = TRUE)
remote_quantiles <- quantile(middlesex$remotework_mayjune2024, probs = seq(0, 1, by = 0.25), na.rm=TRUE)
remote_quantiles## 0% 25% 50% 75% 100%
## 0.0000000 0.3246522 0.3880597 0.4746858 0.7517401middlesex$remote_cat <- cut(middlesex$remotework_mayjune2024, breaks = remote_quantiles, include.lowest = TRUE, labels=FALSE)
palette <- colorNumeric(
palette = colorRampPalette(c("red", "green"))(100),
domain = middlesex$remotework_mayjune2024
)
leaflet(data = middlesex, options = leafletOptions(minZoom = 6)) %>%
addTiles() %>%
addPolygons(
fillColor = ~palette(remotework_mayjune2024),
fillOpacity = 0.7,
color = "black",
weight = 1,
opacity = 1,
highlightOptions = highlightOptions(
weight = 3,
color = "white",
bringToFront = TRUE
),
label = ~paste("Remote Work Proportion:", remotework_mayjune2024)
) %>%
addLegend(
pal = palette,
values = ~remotework_mayjune2024,
title = "Remote Work Proportion",
position = "bottomright"
) %>%
addControl(
html = "<h2>Remote Work Proportion in Middlesex County</h2>",
position = "topright",
className = "map-title"
)palette <- colorFactor(palette = colorRampPalette(c("red", "orange","yellow", "green"))(10), domain = middlesex$remote_cat)
labels <- c("Low", "Medium", "High", "Very High", "Not Available")
leaflet(data = middlesex, options = leafletOptions(minZoom = 6)) %>%
addTiles() %>%
addPolygons(
fillColor = ~palette(remote_cat),
fillOpacity = 0.7,
color = "black",
weight = 1,
opacity = 1,
highlightOptions = highlightOptions(
weight = 3,
color = "white",
bringToFront = TRUE
),
label = ~paste("Remote Work Category:", remote_cat)
) %>%
addLegend(
pal = palette,
values = ~remote_cat,
title = "Remote Work Category",
position = "bottomright",
labFormat = function(type, cuts, p) {
paste(labels)
},
opacity = 1
) %>%
addCircleMarkers(
data = trainstations,
~LONGITUDE, ~LATITUDE,
color = "darkblue",
fillColor = "darkblue",
fillOpacity = 1,
radius = 5,
popup = ~STATION_ID
) %>%
addCircleMarkers(
data = lrstations,
~LONGITUDE, ~LATITUDE,
color = "darkcyan",
fillColor = "darkcyan",
fillOpacity = 1,
radius = 5,
popup = ~STATION
) %>%
addControl(
html = "<h2>Remote Work Category in Middlesex County</h2>",
position = "topright",
className = "map-title"
)