Truncated Density Polygons with ggmap

Question

I am having trouble plotting density maps using R/ggmap. My data looks like this:

> head(W)
        date    lat     lon    dist
1 2010-01-01 31.942 -86.659 292.415
2 2010-01-10 32.970 -84.174  89.121
3 2010-01-17 31.000 -85.363 319.552
4 2010-01-17 31.457 -83.951 258.501
5 2010-01-17 31.073 -81.987 373.915
6 2010-01-17 33.210 -83.149 129.927

And I am plotting it using this:

ggmap(atlanta.map, extent = "panel") +
  geom_point(data = W, aes(x = lon, y = lat)) +
  geom_density2d(data = W, aes(x = lon, y = lat)) +
  stat_density2d(data = W, aes(x = lon, y = lat, fill = ..level..,
    alpha = ..level..), size = 0.01, bins = 8, geom = 'polygon') +
  theme(axis.title = element_blank()) +
  scale_fill_gradient(low = "yellow", high = "red") +
  scale_alpha(range = c(.5, .75), guide = FALSE)

But, although the contours look fine and SOME of the polygons are okay, the polygons which intersect the boundaries get broken into two components: the correct "smooth" portion and a straight line portion which closes the polygon. These two portions meet on the boundary of the map.

My data extends beyond the boundaries of the map, so there is ample information for KDE to derive meaningful estimates of the density all the way up to the boundaries.

Does anybody have an idea what the problem might be? And, more importantly, how I can go about fixing it?

Thanks, Andrew.

Answer 1

With some further googling and hunting around on SO, I have put together a solution.

ggmap(map, extent = "normal", maprange=FALSE) %+% W + aes(x = lon, y = lat) +
    geom_density2d() +
    stat_density2d(aes(fill = ..level.., alpha = ..level..),
                   size = 0.01, bins = 16, geom = 'polygon') +
    scale_fill_gradient(low = "green", high = "red") +
    scale_alpha(range = c(.00, .25), guide = FALSE) +
    coord_map(projection="mercator", 
              xlim=c(attr(map, "bb")$ll.lon, attr(map, "bb")$ur.lon),
              ylim=c(attr(map, "bb")$ll.lat, attr(map, "bb")$ur.lat)) +
    theme(legend.position = "none", axis.title = element_blank())