ggplot map plot fails when limits set with coord_map

Question

I'm trying to limit a map plot to a specific area. coord_map is preferred to scale_x_continuous and y equivalent, since the latter mess up the polygons. But here I'm finding it fails for some reason. Here's my code (downloads a 57kb shapefile):

require(maptools)
require(ggplot2)
download.file('https://dl.dropboxusercontent.com/u/46043231/UK.zip', "uk.zip", method="internal", mode="wb")
unzip('uk.zip')
uk = readShapePoly('uk_outline_1000m.shp')
print(bbox(uk))
        min       max
x  259.9625  655566.4
y 7211.7025 1218558.9
uk2 = fortify(uk)
(p = ggplot(uk2, aes(x=long, y=lat, group=group)) + geom_polygon() + coord_equal())

But when coord_map is used the plot disappears:

p + coord_map(xlim=c(0, 630000), ylim=c(0, 1000000))

Any idea what's going on??

Answer 1

I would try something like this to test a few options.

library(maptools)
library(ggplot2)
library(rgdal)
library(raster)
library(latticeExtra)

Download and read the data

download.file('https://dl.dropboxusercontent.com/u/46043231/UK.zip',
              "uk.zip", method="internal", mode="wb")
unzip('uk.zip')
uk <- readOGR(dsn = getwd(), layer = 'uk_outline_1000m')

Data is assumed to use OSGB 1936 / British National Grid.

More about at SpatialReference

proj4string(uk) <- CRS('+init=epsg:27700') # EPSG 27700
extent(uk)
bb.uk <- as(extent(uk), 'SpatialPolygons') # a spatial object
proj4string(bb.uk) <- CRS('+init=epsg:27700')

Write projected shapefile of uk bbox. I'll write it out to map layers with QGIS. It will be my reference system.

writeOGR(as(bb.uk, 'SpatialPolygonsDataFrame'),
         dsn = getwd(),
         layer = 'bbuk2_bng', driver = 'ESRI Shapefile')

The desired bounding box. Using projected coordinates

bb.uk2 <- as(extent(c(0, 630000), c(0, 1000000)), 'SpatialPolygons')
proj4string(bb.uk2) <- CRS('+init=epsg:27700')

Write projected shapefile of user bbox

writeOGR(as(bb.uk2, 'SpatialPolygonsDataFrame'),
         dsn = getwd(),
         layer = 'bbuk2user_bng', driver = 'ESRI Shapefile')

QGis map using British National Grid EPSG:27700

Plot projected layers

Base plot

plot(uk, col = 'grey50', axes = T, xlim=c(-50000, 705566.4),
     ylim=c(-50000, 1325000))
plot(bb.uk, add = T)
plot(bb.uk2, border = 'red', add = T)

spplot

I took a arbitrary window to expand plot area.

sp::spplot(uk, zcol = 'NAME_ISO', scales = list(draw = TRUE), 
           xlim=c(-50000, 705566.4), ylim=c(-50000, 1325000),
           col.regions="grey90") +
  latticeExtra::layer(sp.polygons(bb.uk, fill = NA, col = 'blue')) +
  latticeExtra::layer(sp.polygons(bb.uk2, fill = NA, col = 'red'))

ggmap with projected layers

uk.df = fortify(uk) # admin 
bbuk.df <- fortify(bb.uk) # country bbox extent
bbuk2.df <- fortify(bb.uk2) # user bbox extent

plot it

p <- ggplot() + geom_polygon(data = uk.df, aes(x=long, y=lat, group=group)) + 
  geom_polygon(data = bbuk.df, aes(x=long, y=lat, group=group),
               colour = 'blue', fill = NA) +
  geom_polygon(data = bbuk2.df, aes(x=long, y=lat, group=group),
               colour = 'red', fill = NA) +
  coord_equal() # cartesian
p

plot it with user bounding box

p + coord_equal(xlim=c(0, 630000), ylim=c(0, 1000000))

Now ggplot with geographic (unprojected) coordinates

WGS84 Unprojected Coordinate Reference System

p.wgs84 <- CRS("+init=epsg:4326") # WGS84 Long Lat

Convert projected layer to WGS84

uk.wgs89 <- spTransform(uk, p.wgs84)

Geographic bbox uk

bbuk.wgs84 <- as(as(extent(uk.wgs89), 'SpatialPolygons'),
                 'SpatialPolygonsDataFrame')

Geographic bbox of user extent

bbuk2.wgs84 <- spTransform(bb.uk2, p.wgs84)
bbuk2.wgs84 <- as(bbuk2.wgs84, 'SpatialPolygonsDataFrame')

Plot it with ggplot and cartesian map

uk.df = fortify(uk.wgs89) # admin 
bbuk.df <- fortify(bbuk.wgs84) # country bbox extent
bbuk2.df <- fortify(bbuk2.wgs84) # user bbox extent

The result is not what I'd expect. I don't figure out why the red box is distorted.