r - Using stat_density2d with coord_map in ggplot2

In my example I am plotting a shapefile for all England & Wales, however I'm interested in the greater London area. If I allow ggplot to simply default and display the entire shapefile then the heatmapping generated by stat_density2d isn't very useful for me, I need to change the scale to show the Greater London area.

To address this I need to set the xlim and ylim to effectively zoom in. I can do this in one of two ways:

1. Use the xlim & ylim standalone functions in ggplot


2. Use xlim & ylim within coord_map

The former causes the heatmapping to scale correctly (ie size of the tiles reduce to give a meaningful indication of density in a smaller area) but the shapefile data gets truncated to the plot leading to some very odd looking polygons.

The latter causes the heatmapping to remain at the scale of the whole country (ie the tiles are massive and so not useful for showing activity density when the map is zoomed in) but the polygons display correctly.

Is there a way to get ggplot to change the xlim/ylim values of a plot without truncating the shapefile/polygon data so that it also rescales the tiles generated by geom_density2d?

I know I could fudge this by only plotting a subset of the shapefile polygons but this seems...messy. It would be way cooler if I could choose anywhere in the country to zoom/centre on and have everything scale correctly without having to chop and change the shapefile data.

This is a follow up to my previous question: How do I create a heatmap in ggplot?

Edit 1: I'm currently investigating a few undocumented arguments that stat_density2d passes to the kde2d 2d density function from the MASS package. There are a few of these but I get the impression a combination of n and lims may help. Basically if I pass a higher value to n (this dictates the number of rows and columns there are in the matrix produced) then the 'resolution' of the 2d density estimation increases at the expense of processing time (this is applied to the entire bounding area of the shapefile rather than just the window I'm limiting the plot to using coord_map). I then use lims to try and limit the matrix to just the plotting area but that doesn't seem to have much effect.