Spatial Analytics

Geography, some economics and data analysis.

Crime data in Hamilton

Phil Donovan / 2019-05-15


I’ve recently discovered the crime data available on the police website. Furthermore, I’m presently looking into housing in Hamilton and thought that I could make good use of the information by creating a interactive chloropleth map of crime in Hamilton.

I’ve included all the code here for easy reuse.

Note that the meshblocks data set is loaded via my nzcensr package.

# Read data
crimes <- path(DATA_DIR, "nz_crime_data_20190513.csv") %>% read_csv

# Clean headers
header_names <- colnames(crimes) %>% 
  str_to_lower %>% 
  str_replace_all(" ", "_") %>% 
  str_replace_all("\\(|\\)", "")

colnames(crimes) <- header_names

After loading and sanitising the tables names, I then filter the crimes data for just Hamilton and the years 2018 and 2019.

crimes_hamilton <- crimes %>% 
  filter(str_detect(territorial_authority, "Hamilton"),
         str_detect(year_month, "18|19")) %>% 
  select(anzsoc_division, anzsoc_group, anzsoc_subdivision,
                          meshblock, year_month) 

Next, I count the number of crimes in each meshblock.

crimes_hamilton_summarise <-
  group_by(crimes_hamilton, anzsoc_division, meshblock) %>% 
  summarise(count = n())

I also create a table of the distinct crime types. This is needed to ensure that every meshblock has a row for each column, otherwise meshblocks would drop out. For example, if meshblock 1 had no burglaries it would disappear when I join the data.

crimes_anzsoc_division <- distinct(crimes, anzsoc_division) %>% mutate(key = 1)

Firstly, load the meshblocks with the meshblock geometry, select the meshblocks from Hamilton, simplify them for quicker plotting, join the distinct crimes to each meshblock, then join the actual crime counts and set all NAs to 0. Finally, split the dataframe into the 6 crimes types.

hamilton_mb <- meshblocks %>% 
  st_set_crs(2193) %>% 
  filter(str_detect(TA2013_NAM, "Hamilton")) %>% 
  st_simplify(dTolerance = 5, preserveTopology = T) %>% 
  select(MB2013, AU2013_NAM) %>% 
  mutate(mb_number = as.character(MB2013) %>% as.integer,
         key = 1) %>% 
  left_join(crimes_anzsoc_division, by = "key") %>% 
  select(-key) %>% 
  left_join(crimes_hamilton_summarise, by = c("mb_number" = "meshblock", "anzsoc_division")) %>% 
  mutate(count = ifelse(is.na(count), 0.0, count)) %>% 
  st_transform(4326) #%>% 
  # st_cast("POLYGON")

hamilton_mb_split <- split(hamilton_mb, hamilton_mb$anzsoc_division)

For each crime type, generate a different palette using the viridis ‘magma’ palette.

Now plot the data on a leaflet map. This is a bit laborious and I’d prefer to do it in a loop to make it more concise and easier to maintain but I was tired!

crime_names <- names(hamilton_mb_split)

leaflet() %>% 
  addProviderTiles("Stamen.Toner") %>% 
  
  # Add 1
  addPolygons(data = hamilton_mb_split[[1]], color = ~palettes[[1]](count), 
                  fillOpacity = 0.5, stroke = FALSE,
                  group=crime_names[1]) %>% 
      addLegend(pal = palettes[[1]], values = hamilton_mb_split[[1]]$count, opacity = 0.5, 
                position = "bottomright", title = "Number of incidents",
                group=crime_names[1]) %>% 
  
  # Add 2
  addPolygons(data = hamilton_mb_split[[2]], color = ~palettes[[2]](count), 
                  fillOpacity = 0.5, stroke = FALSE,
                  group=crime_names[2]) %>% 
      addLegend(pal = palettes[[2]], values = hamilton_mb_split[[2]]$count, opacity = 0.5, 
                position = "bottomright", title = "Number of incidents",
                group=crime_names[2]) %>% 
  
  # Add 3
  addPolygons(data = hamilton_mb_split[[3]], color = ~palettes[[3]](count), 
                  fillOpacity = 0.5, stroke = FALSE,
                  group=crime_names[3]) %>% 
      addLegend(pal = palettes[[3]], values = hamilton_mb_split[[3]]$count, opacity = 0.5, 
                position = "bottomright", title = "Number of incidents",
                group=crime_names[3]) %>% 
  
  #Add 4
  addPolygons(data = hamilton_mb_split[[4]], color = ~palettes[[4]](count), 
                  fillOpacity = 0.5, stroke = FALSE,
                  group=crime_names[4]) %>% 
      addLegend(pal = palettes[[4]], values = hamilton_mb_split[[4]]$count, opacity = 0.5, 
                position = "bottomright", title = "Number of incidents",
                group=crime_names[4]) %>% 
  
  # Add 5
  addPolygons(data = hamilton_mb_split[[5]], color = ~palettes[[5]](count), 
                  fillOpacity = 0.5, stroke = FALSE,
                  group=crime_names[5]) %>% 
      addLegend(pal = palettes[[5]], values = hamilton_mb_split[[5]]$count, opacity = 0.5, 
                position = "bottomright", title = "Number of incidents",
                group=crime_names[5]) %>% 
  
  # Add 6
  addPolygons(data = hamilton_mb_split[[6]], color = ~palettes[[6]](count), 
                  fillOpacity = 0.5, stroke = FALSE,
                  group=crime_names[6]) %>% 
      addLegend(pal = palettes[[6]], values = hamilton_mb_split[[6]]$count, opacity = 0.5, 
                position = "bottomright", title = "Number of incidents",
                group=crime_names[6]) %>% 
  
  # Add controls
  addLayersControl(overlayGroups = names(hamilton_mb_split), 
                   options = layersControlOptions(collapsed = F)) %>% 
  
  hideGroup(crime_names[1]) %>% 
  hideGroup(crime_names[2]) %>% 
  hideGroup(crime_names[3]) %>% 
  hideGroup(crime_names[4]) %>% 
  hideGroup(crime_names[5])