R help - Nov 2012 - [lattice] how to overlay a geographical map on a levelplot?

r-help lattice adepts:

I have a question which is somewhat geospatial, so I posted to r-sig-geo
rather than here:

https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016757.html
> summary: How to overlay a geographical map on each panel in a lattice
> (or Trellis), e.g., of levelplot's? Note I am not inquiring about
> creating choropleth maps[,]
which Sarkar 2008 covers quite well (and which is supported by
latticeExtra::mapplot): I just want to appropriately overlay a
garden-variety boundary-line map on panel viewspaces defined by
longitude and latitude.

A relatively small, quite self-contained example follows the quote
above, in which I plot toy data in the sort of lattice layout I need ...
except that each panel lacks a map appropriate to the spatial domain. If
your competencies extend to that, your assistance would be appreciated.

TIA, Tom Roche <Tom_Roche at pobox.com>

Hello,

Your image is projected, not your panels.

Does it match what you are looking for?

wld <- map('world', xlim=c(west.lon.deg,east.lon.deg), 
ylim=c(south.lat.deg,north.lat.deg),plot=FALSE)

wld <- data.frame(lon=wld$x, lat=wld$y)

state <- map('state', xlim=c(west.lon.deg,east.lon.deg), 
ylim=c(south.lat.deg,north.lat.deg),plot=FALSE)

state <- data.frame(lon=state$x, lat=state$y)

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)
   )
) +
   xyplot(lat ~ lon, state, type='l', lty=1, lwd=0.5,
col='black') +
   xyplot(lat ~ lon, wld, type='l', lty=1, lwd=1, col='black')


Regards,
Pascal


Le 21/11/2012 08:40, Tom Roche a ?crit :
>
> r-help lattice adepts:
>
> I have a question which is somewhat geospatial, so I posted to r-sig-geo
> rather than here:
>
> https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016757.html
>> summary: How to overlay a geographical map on each panel in a lattice
>> (or Trellis), e.g., of levelplot's? Note I am not inquiring about
>> creating choropleth maps[,]
>
> which Sarkar 2008 covers quite well (and which is supported by
> latticeExtra::mapplot): I just want to appropriately overlay a
> garden-variety boundary-line map on panel viewspaces defined by
> longitude and latitude.
>
> A relatively small, quite self-contained example follows the quote
> above, in which I plot toy data in the sort of lattice layout I need ...
> except that each panel lacks a map appropriate to the spatial domain. If
> your competencies extend to that, your assistance would be appreciated.
>
> TIA, Tom Roche <Tom_Roche at pobox.com>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

https://stat.ethz.ch/pipermail/r-sig-geo/2012-November/016757.html
>> summary: How to overlay a geographical map on each panel in a lattice
>> (or Trellis), e.g., of levelplot's?
https://stat.ethz.ch/pipermail/r-help/2012-November/329714.html
> Does [this] match what you are looking for?
Alas, no: unless I'm doing something wrongly, it either throws an error
or redraws the viewport. (Note I omit map=state: the world map is good
enough.) Self-contained example follows:

# start example
library(lattice) # for plotting
library(maps)    # for ... maps

# 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)

# create geographical map: thanks, Pascal Oettli
wld.map <- map('world', plot=FALSE,
               xlim=c(west.lon.deg,east.lon.deg), 
               ylim=c(south.lat.deg,north.lat.deg))
wld.df <- data.frame(lon=wld.map$x, lat=wld.map$y)

sigdigs=3 # significant digits for lattice strips

# plot "appropriately" for atmospheric data: use
# * lattice::levelplot
# * one column, since atmospheric levels stack vertically

# This gets
#> Error in levelplot(concentration ~ longitude * latitude | level, data =
data.3d.df,  :
#>   non-numeric argument to binary operator
# 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)
#    )
# ) + xyplot(lat ~ lon, wld.df,   type='l', lty=1, lwd=1,
col='black')

# This throws no error,
# but merely redraws the entire viewport with the xyplot
# rather than overdrawing each panel of the levelplot.
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)
   )
)
xyplot(lat ~ lon, wld.df,   type='l', lty=1, lwd=1, col='black')