2

I would like to plot some lat-lon data in a filled contour, and then overlay a map of the globe on top. I am trying something like this:

filled.contour(lons, lats, glb.data,
plot.axes={axis(1);axis(2);map(projection='rectangular',parameters=0,add=T)}
)

The contour is plotted fine, and the map shows up, but only as a very tiny little black rectangle in a small area of the contour plot. Of course, I want the map to stretch over the entire area of the filled contour plot.

I've also tried something like:

filled.contour(lons, lats, glb.data,
plot.axes={axis(1);axis(2);map(projection='rectangular',x=lons,y=lats,parameters=0,add=T)}
)

In this case, the map does not plot at all and I get an error message "nothing to draw: all regions out of bounds" People apparently do this successfully all the time, but I haven't been able to get it to work. Can anyone tell me what I am doing wrong? BTW, I don't want to use the lattice package. It will not work with the other stuff my group is doing. Thanks.

bob.sacamento
  • 6,283
  • 10
  • 56
  • 115
  • For the benefit of others, this was answered here: https://stat.ethz.ch/pipermail/r-help/2012-August/321207.html – mdsumner Aug 10 '12 at 01:44
  • Sorry, but that answer gives me a different but similar problem. The map still does not fit on the filled contour region. – bob.sacamento Aug 10 '12 at 15:32
  • that is because you choose a different projection for the map, different to the one your data are in – mdsumner Aug 11 '12 at 22:19

2 Answers2

3

OK. for anyone who cares, I have an answer. Maybe not the best answer, but an answer. The problem is that there are 2 common ways to represent latitude and longitude: degrees and radians. There is a third common way to represent longitude: local time, i.e. time zone time generalized as floating point numbers. And then, to make life more interesting, longitude in degrees is often given as -180 (east) to 180 (West) OR 0 to 360.

Now, when you try to display a map, R automatically picks the lat-lon units for the display based on the innards of the maps package. If you are working in one set of units, but R wants to display the map in another set, you get trouble. Another wrinkle is that different map projections want to be displayed with different units. As I have discovered, the "default" projection likes degrees; it wants to be displayed in a plot rectangle with limits -180 to 180 in longitude, and -90 to 90 in latitude. BUT, the "rectangular" projection, though it looks very much like the default projection, likes radians; it wants to be displayed in a rectangle stretching -pi to pi (longitude) and -pi/2 to pi/2 (latitude). These plot limits are not discussed at all in ?mapproject.

So here's what you need to do: 1) Determine which map projection you want and then determine what limits it is plotted in:

> map(projection=foo,parameters=foo)
> par('usr')   # will return vector of plot limits

2) Know the units your own data is expressed in. If you have just grabbed someone else's data set, you might be surprised here.

3) When you plot ...

filled.contour(lons, lats, glb.data)

if your "lons" and "lats" are in units consistent with the limits your projection expects, then you should be OK, and you should be able to simply add a map with

filled.contour(lons, lats, glb.data,
plot.axes={axis(1,...);axis(2,...);map(projection=foo,parameters=foo,add=T)})

But if your lons and lats are not what the projection needs, then you have a little more work to do. You have to redefine the limits of your plot with a call to par(). This is what worked for me:

filled.contour(lons, lats, glb.data,
plot.axes={axis(1,...);axis(2,...);par(usr=c(-180,180,-90,90));
map(add=T)})

And when I needed the prime meridian at the left end of the plot (i,e, longitude everywhere increasing on the x axis from 0 to 360) I made two calls to par() and map():

filled.contour(lons, lats, glb.data,
plot.axes={axis(1,...);axis(2,...);
par(usr=c(0,360,-90,90));map(add=T);  #eastern hemisphere
par(usr=c(-360,0,-90,90));map(add=T)}   #western
)

There are some other wrinkles too. (e.g. the 'rectangular' projection does not have limits of precisely -pi an +pi.) But if what I have said so far makes sense, you will be able to deal with them.

bob.sacamento
  • 6,283
  • 10
  • 56
  • 115
1

Why not just use image? 

image(lons, lats, glb.data)
library(maps)
map(add = TRUE)

Where did you get the idea to try to put map() in the plot.axes argument?

mdsumner
  • 29,099
  • 6
  • 83
  • 91
  • 1
    I'm not too familiar with image, but as I understand it, it doesn't give a color bar. And it doesn't interpolate, so the colors come out looking "blocky". I tried this and I get a different but related problem: the map appears, but runs off the edge of the window. This got me thinking, though. I can get an approximation to what I want by fiddling with par(usr). There must be a better way, but it works after I fine tune it enough. I got the idea because, AFAICT, "plot.axes" is the usual way to add level lines or coordinate lines to a filled contour. I figured a map would work too.Thanks. – bob.sacamento Aug 09 '12 at 22:43
  • 1
    it does interpolate with useRaster = TRUE, filled.contour sets up its own coordinate system that overrides the base one, which is why the map is offset, use Duncan Murdoch's answer, or use sp::spplot or raster package – mdsumner Aug 10 '12 at 01:45