geom_contour_filled() and ggmap

The question was asked about whether ggplot2’s new geom_contour_filled() can be used with ggmap. In this brief article I illustrate a basic example that shows the answer is “yes”.

library("tidyverse"); theme_set(theme_minimal())
## ── Attaching packages ────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.4
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0
## ── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library("ggmap")
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.

Preliminaries

The first thing I do is create a clean dataset to illustrate the point on:

`%notin%` <- function(lhs, rhs) !(lhs %in% rhs)

violent_crime <- crime %>% 
  filter(
    offense %notin% c("auto theft", "theft", "burglary"),
    -95.39681 <= lon & lon <= -95.34188,
     29.73631 <= lat & lat <=  29.78400
  ) %>% 
  as_tibble()

violent_crime contains this data, among others:

Manually computing the density

geom_contour_filled() requires x, y, and z aesthetics, which works for ggplot2::faithfuld since it already has a pre-computed density column, but most datasets don’t have this, including violent_crime. Of course, we can compute this manually:

kdeout <- violent_crime %>% 
  with( 
    MASS::kde2d(lon, lat, n = 101,
      lims = c(
        scales::expand_range(range(lon), .20),
        scales::expand_range(range(lat), .20)
      )
    )
  )

kde_df <- kdeout %>% 
  .[c("x", "y")] %>% 
  cross_df() %>% 
  rename("lon" = "x", "lat" = "y") %>% 
  mutate(density = as.vector(kdeout$z))

We can use this 2d KDE to make a density surface plot:

ggplot() +
  geom_raster(aes(lon, lat, fill = density), kde_df)

And the same data frame can be used with geom_contour_filled():

ggplot() +
  geom_contour_filled(aes(lon, lat, z = density), kde_df)

This can be used as an overlay of a map with ggmap:

(bbox <- make_bbox(lon, lat, data = violent_crime))
##      left    bottom     right       top 
## -95.39886  29.73408 -95.33919  29.78628
map <- get_stamenmap(bbox, zoom = 14, maptype = "toner-lite")
## Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

ggmap(map, extent = "device") +
  geom_contour_filled(aes(lon, lat, z = density), kde_df, alpha = .65)
## Warning: Removed 3960 rows containing non-finite values (stat_contour_filled).

With stat_density_2d()

It is likely that something like the above is also possible using stat_density_2d(). Outside of ggmap, something like this would work:

ggmap(map, extent = "device") +
  stat_density_2d(
    mapping = aes(lon, lat, fill = stat(level)), 
    data = violent_crime, geom = "raster", 
    alpha = .5, contour = FALSE
  )

However, this doesn’t work here because ggplot2’s geom_raster() demands a Cartesian coordinate system, and ggmap() uses ggmap:::coord_map2(), like ggplot2::coord_map(). Of course, the request is not for a raster overlay, but one made by geom_contour_filled() using isoband. It is not immediately clear to me how to do that. geom_polygon() works, but the semi-transparent polygons stack to become near opaque:

ggmap(map, extent = "device") +
  stat_density_2d(
    mapping = aes(lon, lat, fill = stat(level)), 
    data = violent_crime, geom = "polygon", alpha = .5
  )