[R] [lattice] display a projected map on a layerplot

[Tom Roche]

summary: I can display a lon-lat map on a lattice::layerplot, and I
can display a Lambert conformal conic (LCC) map on a spam::image, but
I can't display an LCC map on a lattice::layerplot. Example follows.
What am I doing wrong?

details:

I've been using `lattice` (via `rasterVis`) successfully to display
global atmospheric data, which works well enough (though I am
definitely intrigued by ggplot/ggmap). Particularly, I am able to
overlay my plots with maps, which is essential for the sort of work
I'm doing. However I'm currently unable to use lattice::layerplot for
some regional data projected LCC: the data plots, but I cannot get a
map to overlay. I can do this by cruder means, but would prefer to
know how to do this in lattice/rasterVis or similar (or ggplot/ggmap).
Two nearly-self-contained examples follow.

The data, ozone_lcc.nc, comes with CRAN package=M3 ... except that M3
supplies it as

system.file("extdata/ozone_lcc.ncf", package="M3")

and that extension currently causes problems for CRAN package=raster.
You can either rename that file (and put it some current working
directory), or you can download (270 kB)

https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ozone_lcc.nc

You can then run the following code:

##### start example #####

library("M3")        # http://cran.r-project.org/web/packages/M3/
library("rasterVis") # http://cran.r-project.org/web/packages/rasterVis/

## Use an example file with projection=Lambert conformal conic.
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3")
lcc.file <- "./ozone_lcc.nc" # unfortunate problem with raster::raster
lcc.proj4 <- get.proj.info.M3(lcc.file)
lcc.proj4   # debugging
# [1] "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000"
# Note +lat_0=40 +lat_1=33 +lat_2=45 for maps::map at projection (below)
lcc.crs <- sp::CRS(lcc.proj4)
lcc.pdf <- "./ozone_lcc.pdf" # for output

## Read in variable with daily ozone.
# o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs)
# ozone_lcc.nc has 4 timesteps, choose 1 at random
o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs, level=1)
o3.raster at crs <- lcc.crs # why does the above assignment not take?
o3.raster # debugging
# class       : RasterLayer 
# band        : 1 
# dimensions  : 112, 148, 16576  (nrow, ncol, ncell)
# resolution  : 1, 1  (x, y)
# extent      : 0.5, 148.5, 0.5, 112.5  (xmin, xmax, ymin, ymax)
# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000 
# data source : .../ozone_lcc.nc 
# names       : O3 
# z-value     : 1 
# zvar        : O3 
# level       : 1 

## Get US state boundaries in projection units.
state.map <- maps::map(
  database="state", projection="lambert", par=c(33,45), plot=FALSE)
#                  parameters to lambert: ^^^^^^^^^^^^
#                  see mapproj::mapproject
state.map.shp <-
  maptools::map2SpatialLines(state.map, proj4string=lcc.crs)

pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
  sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))

dev.off()
# change this as needed to view PDFs on your system
system(sprintf("xpdf %s", lcc.pdf))
# data looks good, but there's no map.

## Try again, with lambert(40,33)
state.map <- maps::map(
  database="state", projection="lambert", par=c(40,33), plot=FALSE)
#                  parameters to lambert: ^^^^^^^^^^^^
#                  see mapproj::mapproject
state.map.shp <-
  maptools::map2SpatialLines(state.map, proj4string=lcc.crs)

pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
  sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
# still no map :-(

dev.off()
system(sprintf("xpdf %s", lcc.pdf))

#####   end example #####

The data looks right, in that it greatly resembles the original
example (from which the above is adapted), which follows my .sig.
However the original-example code *does* draw a map, while the code
above does not. How to fix?

TIA, Tom Roche <Tom_Roche at pobox.com> <Roche.Tom at epa.gov>
---------------original example follows to end of post---------------

# Following adapted from what is installed in my
# .../R/x86_64-pc-linux-gnu-library/2.14/m3AqfigExampleScript.r
# (probably by my sysadmin), which also greatly resembles
# https://wiki.epa.gov/amad/index.php/R_packages_developed_in_AMAD
# which is behind a firewall :-(

## EXAMPLE WITH LAMBERT CONIC CONFORMAL FILE.

library("M3")
library("aqfig") # http://cran.r-project.org/web/packages/aqfig/

## Use an example file with LCC projection: either local download ...
lcc.file <- "./ozone_lcc.nc"
## ... or as installed by package=M3:
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3")
## Choose the one that works for you.

##  READ AND PLOT OZONE FROM FILE WITH LCC PROJECTION.

## Read in variable with daily ozone. Note that we can give dates
## rather than date-times, and that will include all time steps
## anytime during those days.  Or, we can give lower and upper bounds
## and all time steps between these will be taken.
o3 <- M3::get.M3.var(
  file=lcc.file, var="O3",
  ldatetime=as.Date("2001-07-01"), udatetime=as.Date("2001-07-04"))

## Get colors and map boundaries for plot.
library("fields")
col.rng <- tim.colors(20)
detach("package:fields")

## Get state boundaries in projetion units.
map.lines <- M3::get.map.lines.M3.proj(file=lcc.file, database="state")

## Set color boundaries so that they incorporate the complete data range.
z.rng <- range(as.vector(o3$data))
## Make a color tile plot of the ozone for the first day (2001-07-01).
image(o3$x.cell.ctr, o3$y.cell.ctr, o3$data[,,1,1],
      col=col.rng, zlim=z.rng,
      xlab="Projection x-coord (km)", ylab="Projection y-coord (km)")

## Put date-time string and chemical name (O3) into a format I can use
## to label the actual figure.
date.str <- format(o3$datetime[1], "%Y-%m-%d")
title(main=bquote(paste(O[3], " on ", .(date.str), sep="")))

## Put the state boundaries on the plot.
lines(map.lines$coords)

## Add color bar to the right of plot to show what colors go with what
## ozone concentraton.  NOTE: AFTER USING vertical.image.legend(), YOU
## WON'T BE ABLE TO ADD NEW FEATURES TO THE PLOT.  Before making a new
## plot, you should open a new device or turn this device off.
vertical.image.legend(zlim=z.rng, col=col.rng)

dev.off() # close the plot if you haven't already