libraries

library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(ineq)
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)

data

poi_df <- read.csv('poi_nj.csv')
tracts_nj <- st_read('Downloads/tl_2022_34_tract/tl_2022_34_tract.shp')
## Reading layer `tl_2022_34_tract' from data source 
##   `/Users/erandoni/Downloads/tl_2022_34_tract/tl_2022_34_tract.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2181 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

attach tracts to POI data

poi_sf <- st_as_sf(poi_df, coords=c("LONGITUDE", "LATITUDE"), crs=4326)
tracts_nj <- st_transform(tracts_nj, crs = st_crs(poi_sf))

poi_tract_sf <- st_join(poi_sf, tracts_nj, join = st_nearest_feature)

calculate POI diversity

poi_tract_sf <- poi_tract_sf %>%
  mutate(NAICS_sector = substr(NAICS_CODE, 1, 2))

diversity_of_POIs <- poi_tract_sf %>%
  group_by(GEOID) %>%
  summarise(gini_coef = ineq(table(NAICS_sector), type = 'Gini'))

diversity_of_POIs <- st_drop_geometry(diversity_of_POIs)

diversity_of_POIs <- tracts_nj %>%
  left_join(diversity_of_POIs, by = 'GEOID')

calculate alcohol establishments, dining establishments, marijuana establishments

poi_tract_sf <- poi_tract_sf %>%
  mutate(NAICS_industry = substr(NAICS_CODE, 1, 4))

cannabis_categories = c('Cannabis Dispensary', 'Cannabis Clinic,Cannabis Dispensary', 'Cannabis Clinic,Cannabis Dispensary,Medical Cannabis Referrals')

establishment_summary <- poi_tract_sf %>%
  group_by(GEOID) %>%
  summarise(
    dining_establishments = sum(NAICS_industry == 7225),
    alcohol_establishments = sum(NAICS_industry %in% c(7224, 4248)),
    marijuana_establishments = sum(CATEGORY_TAGS %in% cannabis_categories))

establishment_summary <- st_drop_geometry(establishment_summary)

diversity_of_POIs <- merge(diversity_of_POIs, establishment_summary, by = 'GEOID')

leaflet colored by poi diversity, with count of specific establishments

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

diversity_of_POIs$gini <- round(diversity_of_POIs$gini_coef, 3)

leaflet(data = diversity_of_POIs) %>%
  addTiles() %>%
  addPolygons(
    color = 'black',
    weight = 1,
    label = ~paste("gini:", gini, "dining establishments:", dining_establishments, "alcohol establishments:", alcohol_establishments, "marijuana establishments:", marijuana_establishments),
    fillColor = ~pal(gini_coef),
    fillOpacity = 0.7
  ) %>%
  addLegend(pal = pal, values = ~gini_coef, title = "Diversity of POIs by Census Tract")