how to overlay maps on lattice panels?

summary: How to overlay a geographical map on each panel in a lattice?
(Note I am not seeking to create a choropleth map.)

details:

I have atmospheric data specifying gas concentrations over a 3D space
with dimensions longitude, latitude, and (vertical) level. Since the
data is numeric, and I need to plot it by level, lattice::levelplot
seems to be the "tool for the job." However I also need to spatially
manipulate the data, and therefore to display data before-and-after, and
the best way to do that seems to be to overlay the data with a map.
(Am I missing something?)

From the excellent open-access book on package=lattice

ISBN 978-0-387-75968-5
e-ISBN 978-0-387-75969-2
http://dx.doi.org/10.1007/978-0-387-75969-2

I know how to plot a choropleth map with latticeExtra::mapplot, but
that's not what I want, which is just to overlay a lattice (or Trellis)
of levelplot's with the lines of a geographical map. How to do that?

A relatively small, quite self-contained example follows, in which I
plot toy data in the sort of lattice layout I need. The toy data covers
a horizontal space over North America; each panel has bounds=

longitudes: -130? to -59.5?
latitudes:  23.5? to 58.5?

corresponding to this geography

https://github.com/TomRoche/cornbeltN2O/wiki/images/Cmaq_aqmeii_domain.png

So how to (programmatically) overlay something like that image over each
panel in the lattice generated below? Your assistance is appreciated!

# start example
library(lattice)

# Not too many points, for this example
# -130? to -59.5: (130-59.5)/1.5=47
# 23.5? to 58.5?: 58.5-23.5     =35

# define longitudes
west.lon.deg=-130
#east.lon.deg=-115    # debugging
east.lon.deg=-59.5
step.lon.deg=1.5
lons  <- c(seq(west.lon.deg, east.lon.deg, step.lon.deg))
n.lon <- length(lons)

# define latitudes
south.lat.deg=23.5
#north.lat.deg=30.5    # debugging
north.lat.deg=58.5
step.lat.deg=1.0
lats  <- c(seq(south.lat.deg, north.lat.deg, step.lat.deg))
n.lat <- length(lats)

# define vertical levels
bot.lev.mb=1000
top.lev.mb=50
n.lev=4
levs <- c(seq(from=bot.lev.mb, to=top.lev.mb, length.out=n.lev))

# define concentrations
grid.3d.len <- n.lon * n.lat * n.lev
concs <- c(1:grid.3d.len)

# define grid
grid.3d.df <- expand.grid(longitude=lons, latitude=lats, level=levs)
# add concentrations as data
data.3d.df <- cbind(grid.3d.df, concentration=concs)
# debugging
head(data.3d.df)
tail(data.3d.df)

sigdigs=3 # significant digits
# plot "appropriately" for atmospheric data: use
# * lattice::levelplot
# * one column, since atmospheric levels stack vertically
levelplot(
  concentration ~ longitude * latitude | level, 
  data=data.3d.df, layout=c(1,n.lev),
  levs=as.character(round(data.3d.df[['level']], 1)),
  strip=FALSE,
  strip.left=strip.custom(
    factor.levels= # thanks, David Winsemius
      as.character(signif(unique(data.3d.df[['level']]), sigdigs)),
    strip.levels=TRUE,
    horizontal=TRUE,
    strip.names=FALSE,
    # gotta shrink strip text size to fit strip width:
    # more on that separately
    par.strip.text=list(cex=0.5)
  )
)
# end example

TIA, Tom Roche <Tom_Roche at pobox.com>