[R-sig-Geo] how to overlay maps on lattice panels?

[Edzer Pebesma]

http://r-spatial.sourceforge.net/gallery/

shows a few options how this can be done with sp::spplot

On 11/20/2012 10:46 PM, Tom Roche wrote:
> 
> 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>
> 
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> 

-- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Weseler Straße 253, 48151 Münster, Germany. Phone: +49 251
8333081, Fax: +49 251 8339763  http://ifgi.uni-muenster.de
http://www.52north.org/geostatistics      e.pebesma at wwu.de