problem set 55

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE

library(osmdata)

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(ggspatial)
library(sfhotspot)

options(osmdata.use_cache = TRUE)
options(osmdata.base_url = "https://overpass-api.de/api/interpreter")

fixed_plot_aspect <- function() {
  coord_sf(datum = NA)
}

wards <- read_sf("https://mpjashby.github.io/crimemappingdata/nottingham_wards.gpkg") |>
  st_transform("EPSG:27700") |>
  dplyr::filter(ward_name %in% c("Castle",
                                 "Lenton & Wollaton East",
                                 "Meadows"))

burglaries <- readr::read_csv(
  "https://mpjashby.github.io/crimemappingdata/nottingham_burglary.csv.gz"
) |>
  st_as_sf(coords = c("longitude", "latitude"), crs = "EPSG:4326") |>
  st_transform("EPSG:27700") |>
  st_intersection(wards)

## Rows: 1795 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (2): location, lsoa_code
## dbl  (2): longitude, latitude
## date (1): month
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

nottingham_buildings <- wards |>
  st_transform("EPSG:4326") |>
  st_bbox() |>
  opq() |>
  add_osm_feature(key = "building") |>
  osmdata_sf()

library(tidyverse)
library(sf)
library(osmdata)
library(ggspatial)
library(sfhotspot)

options(osmdata.use_cache = TRUE)
options(osmdata.base_url = "https://overpass-api.de/api/interpreter")

fixed_plot_aspect <- function() {
  coord_sf(datum = NA)
}

wards <- read_sf("https://mpjashby.github.io/crimemappingdata/nottingham_wards.gpkg") |>
  st_transform("EPSG:27700")

robbery <- read_csv("https://mpjashby.github.io/crimemappingdata/nottingham_robbery.csv.gz") |>
  st_as_sf(coords = c("longitude", "latitude"), crs = "EPSG:4326") |>
  st_transform("EPSG:27700")

## Rows: 555 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (2): location, lsoa_code
## dbl  (2): longitude, latitude
## date (1): month
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

robbery_gistar <- robbery |>
  hotspot_gistar(cell_size = 100, bandwidth_adjust = 0.25, quiet = TRUE) |>
  dplyr::filter(gistar > 0, pvalue < 0.05) |>
  st_intersection(wards)

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

ggplot() +
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(aes(fill = kde), data = robbery_gistar, alpha = 0.8, colour = NA) +
  geom_sf(data = wards, colour = "grey70", fill = NA) +
  scale_fill_gradient(
    low = rgb(0.8, 0.8, 1),
    high = rgb(0, 0, 1),
    breaks = range(pull(robbery_gistar, kde)),
    labels = c("lower", "higher")
  ) +
  fixed_plot_aspect() +
  labs(
    title = "Nottingham Robbery Hotspots",
    subtitle = "Density of robbery at statistically significant hotspots",
    caption = "Contains public sector information licensed under the Open Government Licence v3.0. Map data from OpenStreetMap.",
    fill = "Robbery density"
  ) +
  theme_void() +
  theme(
    plot.caption = element_text(colour = "grey40", hjust = 0),
    plot.subtitle = element_text(margin = margin(t = 6, b = 6)),
    plot.title = element_text(colour = "grey50", face = "bold", size = 16)
  )