How can I add a lat-long grid to a projected map?

Question

EDIT: question really pertains to how one can add lat/long gridlines to a projected map. I changed title to match.

I have some layers in geographic coords. I want to plot them in LCC projection, but have a geographic (lat/long) grid. From mapproj I can use map.grid() to add a grid, with the limits set by the lim argument. It takes a vector or a range object:

a vector of 4 numbers specifying limits: c(lon.low, lon.high, lat.low, lat.high). lim can also be a list with a component named range, such as the result of map, from which limits are taken.

I construct my map by clipping a large vector layer with a clipping polygon:

myPoly <- readOGR(dsn=".", layer="myPolygon")  # just a shapefile in geographic coords
library(raster)  # To convert an 'extent' object to a SpatialPolygons object
cp <- as(extent(146, 149, -39, -37.5), "SpatialPolygons")
proj4string(cp) <- CRS(proj4string(myPoly))  # copy from shapefile

# Transform and plot:
lcc <- CRS("+init=epsg:3111")
myPoly.proj <- spTransform(myPoly, lcc)
cp.proj <- spTransform(cp, lcc)  # transform the clip box
myPoly.proj.clip <- gIntersection(myPoly.proj, cp.proj, byid=TRUE)
plot(myPoly.proj.clip) 

# Then finally, add a lat/long grid:
map.grid(lim=as.vector(cp.proj@bbox), labels=TRUE)

That last line is not correct, as the @bbox returned is xmin, ymin, xmax, ymax, but needs to be in xmin, xmax, ymin, ymax. There must be a simple solution to all this, but as usual I am lost in the vortex. I could manually create a limits vector, but really?

don't use @ - try lim = as.vector(t(bbox(cp.proj))), but I don't know if map.grid plays well with Spatial objects after plotting, for example how does map.grid know that you want longlat grid, and what projection your data is in? — mdsumner
@mdsumner - thanks. I have no idea if map.grid() knows what to do! Just trying to work out anything seems hard for me at the moment. Any suggestions welcome. — a different ben
@mdsumner - ?map.grid() says: Draw a latitude/longitude grid on a projected map, but I suppose it's talking about a map produced by the map function. — a different ben
I think that's right. Pretty sure mapproj defaults to the unit sphere, so you're going to be pushing it uphill. See ?gridlines for working with sp. (This is a nightmare I know). — mdsumner
Ah-ha! rgdal has llgridlines, which does it I think! @mdsumner — a different ben

mdsumner mdsumner · Accepted Answer · 2014-07-08T05:54:39

EDIT: the OP points out rgdal::llgridlines which is a better solution.

You are using context from sp/rgdal which uses a different system to that of mapproj/maps.

Try this (untested):

library(rgdal)
gl <- gridlines(myPoly)
cp.gl <- spTransform(gl, lcc)
plot(cp.gl, add = TRUE)

See ?gridlines for more on using this with labels. I find it works well as long as you stay away from circumpolar maps.

How can I add a lat-long grid to a projected map?

1 Answers