Adaptive bandwidth kernel density estimate about vacant-lots-proj HOT 2 CLOSED

# 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