3

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?

a different ben
  • 3,900
  • 6
  • 35
  • 45
  • 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 Jul 08 '14 at 05:37
  • @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 Jul 08 '14 at 05:41
  • @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 Jul 08 '14 at 05:44
  • 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 Jul 08 '14 at 05:51
  • Ah-ha! `rgdal` has `llgridlines`, which does it I think! @mdsumner – a different ben Jul 08 '14 at 05:55
  • OK, llgridlines seems good, but labelling is fixed. Starting to shiver and sweat here. – a different ben Jul 08 '14 at 06:11

1 Answers1

1

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.

mdsumner
  • 29,099
  • 6
  • 83
  • 91