load and filter data, plot all pedestrian crashes

crashes <- read.csv('/Users/erandoni/mechkar/FARS_NJ-NY-CT-18-22.csv')
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(leaflet)

crashes <- crashes %>% filter(LONGITUD < 500)
nyc_crashes <- crashes %>% filter(STATENAME == "New York") %>% filter(COUNTYNAME %in% c("NEW YORK (61)", "BRONX (5)", "KINGS (47)", "QUEENS (81)"))

leaflet(data = nyc_crashes) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~LONGITUD,
    lat = ~LATITUDE,
    radius = 3,
    color= 'red',
    fillColor = 'red',
    fillOpacity = 0.5
  ) %>%
  setView(
    lng = mean(nyc_crashes$LONGITUD),
    lat = mean(nyc_crashes$LATITUDE),
    zoom = 10
  )

create heatmap using kernel density estimation

library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
library(sp)
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
nyc_shp <- read_sf('/Users/erandoni/mechkar/cb_2018_36_tract_500k/cb_2018_36_tract_500k.shp')
nyc_shp <- nyc_shp %>% filter(COUNTYFP %in% c('061', '005', '047', '081'))

raster_extent <- extent(nyc_shp)
grid <- raster(raster_extent, res = 0.0001)

coordinates(nyc_crashes) <- ~LONGITUD + LATITUDE
proj4string(nyc_crashes) <- CRS("+proj=longlat +datum=WGS84")
proj4string(grid) <- proj4string(nyc_crashes)

nyc_crashes$crash <- 1
crash_raster <- rasterize(nyc_crashes, grid, field = 'crash', fun = sum)
crash_raster[is.na(crash_raster)] <- 0

leafletplot <- function(input_raster) {
  turbo_palette <- colorRampPalette(c("lightgreen", "yellow", "orange", "red", "darkred"))
  turbo_colors <- turbo_palette(100)
  palette_min <- min(values(input_raster))
  palette_max <- max(values(input_raster))
  pal <- colorNumeric(
    palette = turbo_colors,
    na.color = "transparent",
    domain = c(palette_min, palette_max)
  )
  input_raster <- crop(input_raster, nyc_shp)
  input_raster <- mask(input_raster, nyc_shp)
  leaflet() %>%
    addProviderTiles(provider = providers$OpenStreetMap) %>%
    addRasterImage(input_raster, colors = pal, opacity = 0.7)
}

k <- focalWeight(crash_raster, d = 0.005, type = "Gauss")
h <- focal(crash_raster, w = k, fun = sum, na.rm = TRUE)
h[is.na(h)] <- 0

leafletplot(h)

create heatmap using Getis-Ord

library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
grid <- raster(raster_extent, res = 0.0025)
proj4string(grid) <- proj4string(nyc_crashes)

getisraster <- rasterize(nyc_crashes, grid, field = 'crash', fun = sum)
getisraster[is.na(getisraster)] <- 0
getisgrid <- rasterToPolygons(getisraster)

neighbors = poly2nb(getisgrid)
weighted_neighbors = nb2listw(neighbors, zero.policy=TRUE)

getisgrid$HOTSPOT = as.vector(localG(getisgrid$layer, weighted_neighbors))
breaks = c(-20, -1.96, -1, 1, 1.96, 20)
palette=c("#0000FF80", "#8080FF80", "#FFFFFF80", "#FF808080", "#FF000080")
getisgrid$color <- palette[cut(getisgrid$HOTSPOT, breaks = breaks, include.lowest = TRUE)]

leaflet(data = getisgrid) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addPolygons(
    fillColor = ~color,
    weight = 0,
    fillOpacity = 1,
    popup = ~paste("Z-Score:", round(HOTSPOT, 2), "P-value:", round(2 * (1-pnorm(abs(HOTSPOT))), 4))
  ) %>%
  addLegend(
    colors = palette,
    labels = c("Significant Cold Spot", "Cold Spot", "Neutral", "Hot Spot", "Significant Hot Spot"),
    title = "Getis-Ord Hotspot Significance"
  )

hotspots of low health areas (hot spot => low health, cold spot => good health)

sdoh <- read.csv('/Users/erandoni/mechkar/combined_health_crash.csv')
nyc_shp <- merge(nyc_shp, sdoh, by = 'GEOID', all.x = TRUE)
nyc_shp <- nyc_shp %>% filter(!is.na(health_intensity))

neighbors <- poly2nb(nyc_shp)
weights <- nb2listw(neighbors, style = "W", zero.policy = TRUE)
nyc_shp$hotspot <- as.vector(localG(nyc_shp$health_intensity, weights))

breaks = c(-20, -1.96, -1, 1, 1.96, 20)
palette=c("#0000FF80", "#8080FF80", "#FFFFFF80", "#FF808080", "#FF000080")
col = palette[cut(nyc_shp$hotspot, breaks)]

nyc_shp$color <- palette[cut(nyc_shp$hotspot, breaks = breaks, include.lowest = TRUE)]

leaflet(data = nyc_shp) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addPolygons(
    fillColor = ~color,
    color = 'black',
    weight = 1,
    fillOpacity = 1,
    popup = ~paste("Z-Score:", round(hotspot, 2), "P-value:", round(2 * (1-pnorm(abs(hotspot))), 4))
  ) %>%
  addLegend(
    colors = palette,
    labels = c("Significant Cold Spot", "Cold Spot", "Neutral", "Hot Spot", "Significant Hot Spot"),
    title = "Getis-Ord Hotspot Significance"
  )
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'

hotspots of asthma

nyc_shp <- nyc_shp %>% filter(!is.na(CASTHMA_CrudePrev))

neighbors <- poly2nb(nyc_shp)
weights <- nb2listw(neighbors, style = "W", zero.policy = TRUE)
nyc_shp$hotspot <- as.vector(localG(nyc_shp$CASTHMA_CrudePrev, weights))

breaks = c(-20, -1.96, -1, 1, 1.96, 20)
palette=c("#0000FF80", "#8080FF80", "#FFFFFF80", "#FF808080", "#FF000080")
col = palette[cut(nyc_shp$hotspot, breaks)]

nyc_shp$color <- palette[cut(nyc_shp$hotspot, breaks = breaks, include.lowest = TRUE)]

leaflet(data = nyc_shp) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addPolygons(
    fillColor = ~color,
    color = 'black',
    weight = 1,
    fillOpacity = 1,
    popup = ~paste("Z-Score:", round(hotspot, 2), "P-value:", round(2 * (1-pnorm(abs(hotspot))), 4))
  ) %>%
  addLegend(
    colors = palette,
    labels = c("Significant Cold Spot", "Cold Spot", "Neutral", "Hot Spot", "Significant Hot Spot"),
    title = "Getis-Ord Hotspot Significance"
  )
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'

hotspots of diabetes

nyc_shp <- nyc_shp %>% filter(!is.na(DIABETES_CrudePrev))

neighbors <- poly2nb(nyc_shp)
weights <- nb2listw(neighbors, style = "W", zero.policy = TRUE)
nyc_shp$hotspot <- as.vector(localG(nyc_shp$DIABETES_CrudePrev, weights))

breaks = c(-20, -1.96, -1, 1, 1.96, 20)
palette=c("#0000FF80", "#8080FF80", "#FFFFFF80", "#FF808080", "#FF000080")
col = palette[cut(nyc_shp$hotspot, breaks)]

nyc_shp$color <- palette[cut(nyc_shp$hotspot, breaks = breaks, include.lowest = TRUE)]

leaflet(data = nyc_shp) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addPolygons(
    fillColor = ~color,
    color = 'black',
    weight = 1,
    fillOpacity = 1,
    popup = ~paste("Z-Score:", round(hotspot, 2), "P-value:", round(2 * (1-pnorm(abs(hotspot))), 4))
  ) %>%
  addLegend(
    colors = palette,
    labels = c("Significant Cold Spot", "Cold Spot", "Neutral", "Hot Spot", "Significant Hot Spot"),
    title = "Getis-Ord Hotspot Significance"
  )
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'

hotspots of obesity

nyc_shp <- nyc_shp %>% filter(!is.na(OBESITY_CrudePrev))

neighbors <- poly2nb(nyc_shp)
weights <- nb2listw(neighbors, style = "W", zero.policy = TRUE)
nyc_shp$hotspot <- as.vector(localG(nyc_shp$OBESITY_CrudePrev, weights))

breaks = c(-20, -1.96, -1, 1, 1.96, 20)
palette=c("#0000FF80", "#8080FF80", "#FFFFFF80", "#FF808080", "#FF000080")
col = palette[cut(nyc_shp$hotspot, breaks)]

nyc_shp$color <- palette[cut(nyc_shp$hotspot, breaks = breaks, include.lowest = TRUE)]

leaflet(data = nyc_shp) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addPolygons(
    fillColor = ~color,
    color = 'black',
    weight = 1,
    fillOpacity = 1,
    popup = ~paste("Z-Score:", round(hotspot, 2), "P-value:", round(2 * (1-pnorm(abs(hotspot))), 4))
  ) %>%
  addLegend(
    colors = palette,
    labels = c("Significant Cold Spot", "Cold Spot", "Neutral", "Hot Spot", "Significant Hot Spot"),
    title = "Getis-Ord Hotspot Significance"
  )
## Warning: sf layer has inconsistent datum (+proj=longlat +datum=NAD83 +no_defs).
## Need '+proj=longlat +datum=WGS84'