The original map can be found in my previous post, Spatial Data Visualization with R.

As a review, I use get_map function in the ggmap package to grab a map of Boston from Google Maps. The “crime report” data can be found at https://data.cityofboston.gov/. In the code below, I bring in crime data as a csv file. The data contains one observation per crime. It includes a description, crime category (drug, traffic violation, etc), and the longitude/latitude coordinates of the crime scene.

I added density areas using stat_density2d function. I feed this function a set of coordinates using the x= and y= parameters. The alpha parameter adjust transparency with 1 being perfect solid and 0 being fully transparent. I set the fill to vary with the number of points in the shaded area (..level..). I also specify bins=8, which gives us 7 shades of blue. The density areas can be interpreted as follows: all the shaded areas together contain 7/8 of drug crimes in the data. Each shade represents 1/8 of drug crimes in the data. Since all shades represent the same proportion of the data, the smaller the area of a particular shade, the higher the arrest density.

R code is given below.

//update: Thanks to Verbal for the meaningful distinction between crimes and arrest. It’s not really clear what this data actually tracks. I’m sure that this is reported crime, as called in by people. I don’t think every crime here leads to an arrest. I could be wrong.

#install.packages('ggmap')
library(ggmap)

## load data
lib='/...directory with CSV crime file.../'
setwd(lib)
boston=read.csv('Crime_Incident_Reports.csv')

## subset to drug crimes
drugs=boston[which(boston$INCIDENT_TYPE_DESCRIPTION=='DRUG CHARGES' & boston$Year=='2014'),]

## plot drug and shooting crimes
bos_plot=ggmap(get_map('Boston, Massachusetts',
                zoom=13,
                source='google',
                maptype='terrain'))

## Density areas
bos_plot +
  # density areas
  stat_density2d(data = drugs,
                 aes(x = drugs$Lat, y = drugs$Long, alpha=.75,fill=..level..),
                 bins = 8,
                 geom = 'polygon')+
  #density legend
  guides(fill = guide_colorbar(barwidth = 1, barheight = 10)) +
  # crime data points
  geom_point(data=drugs,aes(x=drugs$Lat,y=drugs$Long),
             col='gray', alpha=.5,size=1)+
  scale_alpha(guide = FALSE)+
  # Labels/Title
  xlab('')+ylab('')+
  ggtitle('Drug Crimes - 2014 - Boston, MA')