Adaptive bandwidth kernel density estimate

[nlebovits]

Currently, the gun crime layer in this project is based on a fixed bandwidth kernel density estimate. It creates a raster, the values of which are then extracted to the centroids of vacant properties in the primary dataset. Before launching this dashboard, we need to switch to an adaptive bandwidth kernel density estimate. Full justification can be found in this memo from the Philadelphia District Attorney's Office.

However, there is no implementation in Python that I have been able to find. PySal does allow for the creation of adaptive bandwidth weights, but I have not figured out how to use this in our specific case. The best example I can find is implemented in the spatstat package in R and includes references to the paper that explains the methodology. My suggestion would be to use that as the basis for writing a custom function in Python.

[nlebovits]

Reprex here based on the code that I was given by the DAO (this is in R):

# Load in packages

library(tidyverse)
library(tidylog)
library(tmap)
library(rgdal);

# Load in Data

## Criminal Incidents 

crime_incidents <- read_csv("https://phl.carto.com/api/v2/sql?filename=incidents_part1_part2&format=csv&skipfields=cartodb_id,the_geom,the_geom_webmercator&q=SELECT%20*%20,%20ST_Y(the_geom)%20AS%20lat,%20ST_X(the_geom)%20AS%20lng%20FROM%20incidents_part1_part2%20WHERE%20dispatch_date_time%20%3E=%20%272015-01-01%27%20AND%20dispatch_date_time%20%3C%20%272021-12-31%27")

## Firearm related incidents

firearm_incid <- crime_incidents %>% 
  # Create year variable for filtering to study period
  mutate(crime_year = year(dispatch_date)) %>% 
  # Remove NA lat and lng coords since that will violate Sf conversion
  filter(!is.na(lat),
         !is.na(lng),
         # grab all firearm related incidents 
         str_detect(text_general_code,"Firearm"),
         # Filter to study period
         crime_year >= 2015 & crime_year <= 2021) %>% 
  # Transform to spatial object
  st_as_sf(coords = c("lng","lat"), crs = 4326) %>%
  # Project data to proper crs (NAD83 / UTM zone 18N)
  st_transform(crs = 26918) 

## City Limits

city_lim <- read_sf("https://opendata.arcgis.com/datasets/405ec3da942d4e20869d4e1449a2be48_0.geojson") %>% 
  # Project data to proper crs (NAD83 / UTM zone 18N)
  st_transform(crs = 26918)

# Convert all data from sf to Spatialpolygons format. 
# This needs to be done due to spatstat package (Adaptive-Bandwidth KDE) only 
# taking this spatial format.

## City Limits
city_limit_sp <- city_lim %>% 
  as(Class = "Spatial") %>% 
  spTransform(CRS("+init=epsg:26918"))

## Firearm related incidents

gun_crime <- firearm_incid %>% 
  as(Class = "Spatial") %>% 
  spTransform(CRS("+init=epsg:26918"))

# Kernel Density with fixed bandwidth

## Creates a bounding box based on the extents of the city limits polygon, and will
## be used for making the raster output later when mapping

bounding_box <- city_limit_sp@bbox

## Creates kernel density raster
kde_output <- kernelUD(gun_crime, h="href", grid = 1000)

# converts to raster
kde <- raster(kde_output)

# sets projection to (NAD83 / UTM zone 18N)

projection(kde) <- CRS("+init=EPSG:26918")

# mask the raster by the output area polygon
masked_kde <- mask(kde, city_limit_sp)

# maps the masked raster, also maps white output area boundaries
kde_map <- tm_shape(masked_kde, bbox = bounding_box) +
  tm_raster("ud", style = "fisher", n = 6,
            legend.show = TRUE, palette = "viridis") +
  tm_shape(city_limit_sp) +
  tm_borders(alpha=.3, col = "white") +
  tm_layout(frame = FALSE) +
  tm_legend(legend.outside=TRUE)

# Kernel Density with Adaptive Bandwidth 

## Creates window or bounding box for the point to be contained
w <- as.owin(city_limit_sp)

## makes the gun crime object into a spatstat object and confines it to the bbox
pts_ppp <- as.ppp.data.frame(coordinates(gun_crime), w) %>% 
  ## "Applies independent random displacements to each point in a point pattern."
  rjitter(retry=TRUE, nsim=1, drop=TRUE)

## Shows the objects together (city limits and gun crimes)
plot(pts_ppp)

## Identify global bandwidth
## "Computes adaptive kernel estimates of spatial density/intensity using a 3D 
##  FFT for multiple global bandwidth scales."

gun_multi <- multiscale.density(pts_ppp,h0=1)

plot(gun_multi)

## calculate kernel density estimates using adaptive bandwidth
Z <- densityAdaptiveKernel(pts_ppp, h0=0.1, at=c("pixel"))

plot(Z, main="Adaptive kernel estimate")

## convert kde IM (Image) format to raster
kde_convert <- raster(Z)

## set crs (NAD83 / UTM zone 18N)
projection(kde_convert) <- CRS("+init=EPSG:26918")

## mask the raster by the output area polygon
masked_kde2 <- mask(kde_convert, city_limit_sp)

## maps the masked raster, also maps white output area boundaries
akde_map <- tm_shape(masked_kde2, bbox = bounding_box) +
  tm_raster("layer", style = "fisher", n = 6,
            legend.show = TRUE, palette = "viridis") +
  tm_shape(city_limit_sp) +
  tm_borders(alpha=.3, col = "white") +
  tm_layout(frame = FALSE) +
  tm_legend(legend.outside=TRUE)

# Plot Fixed and Adaptive Bandwidth Kernell Density Estimations

## Kernel Density (Fixed bandwidth)

kde_map

## Kernel Density (adapative bandwidth)

akde_map

^{Created on 2023-08-07 with reprex v2.0.2}

[brandonfcohen1]

closing this for now, will re-open when @nlebovits wants to deep dive here