[R-sig-Geo] how to display a projected map on an rasterVis::layerplot?

[Tom Roche]

summary: How to convert from a simple (x,y) matrix (as produced by
M3::get.map.lines.M3.proj) to a SpatialLines, or to coordinates
consumable by latticeExtra::layer?

details:

Perhaps I should have asked the question above first, instead of
putting it at the end of

http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot

where it is obviously overlooked :-( The problem is also explained
there:

http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot (formatted for mail)
> ## Read in variable with daily ozone.
> 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?
> # start debugging
> o3.raster
> summary(coordinates(o3.raster)) # thanks, Felix Andrews!
> M3::get.grid.info.M3(lcc.file)
> #   end debugging

> # > o3.raster
> # 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 

> # > summary(coordinates(o3.raster))
> #        x                y         
> #  Min.   :  1.00   Min.   :  1.00  
> #  1st Qu.: 37.75   1st Qu.: 28.75  
> #  Median : 74.50   Median : 56.50  
> #  Mean   : 74.50   Mean   : 56.50  
> #  3rd Qu.:111.25   3rd Qu.: 84.25  
> #  Max.   :148.00   Max.   :112.00  

> # > M3::get.grid.info.M3(lcc.file)
> # $x.orig
> # [1] -2736000
> # $y.orig
> # [1] -2088000
> # $x.cell.width
> # [1] 36000
> # $y.cell.width
> # [1] 36000
> # $hz.units
> # [1] "m"
> # $ncols
> # [1] 148
> # $nrows
> # [1] 112
> # $nlays
> # [1] 1

> # The grid (`coordinates(o3.raster)`) is integers, because this
> # data is stored using the IOAPI convention. IOAPI makes the grid
> # Cartesian by using an (approximately) equiareal projection (LCC)
> # and abstracting grid positioning to metadata stored in netCDF
> # global attributes. Some of this spatial metadata can be accessed
> # by `M3::get.grid.info.M3` (above), and some can be accessed via
> # the coordinate reference system (`CRS`, see `lcc.proj4` above)

However I can

> ## Get state boundaries in projection units.
> map.lines <- M3::get.map.lines.M3.proj(
>   file=lcc.file, database="state", units="m")
> # start debugging
> class(map.lines)
> summary(map.lines)
> summary(map.lines$coords)
> #   end debugging

> # >     class(map.lines)
> # [1] "list"
> # >     summary(map.lines)
> #        Length Class  Mode     
> # coords 23374  -none- numeric  
> # units      1  -none- character
> # >     summary(map.lines$coords)
> #        x                  y           
> #  Min.   :-2272238   Min.   :-1567156  
> #  1st Qu.:   94244   1st Qu.: -673953  
> #  Median :  913204   Median :  -26948  
> #  Mean   :  689997   Mean   :  -84644  
> #  3rd Qu.: 1744969   3rd Qu.:  524531  
> #  Max.   : 2322260   Max.   : 1265778  
> #  NA's   :168        NA's   :168    

Note how these grid coordinates derive linearly from the IOAPI
metadata above, and to the raster bounds (farther above): e.g.,

map.lines.XY.min.raw <-
  c(-2272238,-1567156) # from summary(map.lines$coords)
x.orig       <- -2736000
x.cell.width <- 36000
y.orig       <- -2088000
y.cell.width <- 36000

map.lines.XY.min.IOAPI <- c(
  (map.lines.XY.min.raw[1] - x.orig)/x.cell.width,
  (map.lines.XY.min.raw[2] - y.orig)/y.cell.width
)
map.lines.XY.min.IOAPI
# [1] 12.88228 14.46789

Looking @ 

http://i.stack.imgur.com/eijZI.png (example 2)

eyeball the intersection of the latitude of south tip of Texas with
the westernmost longitude of California. Then look at the
corresponding location in

http://i.stack.imgur.com/z7OA0.png (example 1)

you will see that corresponding example1 coordinates are approximately
(12.88228,14.46789) as computed above. Hardly a proof :-) but the this
certainly shows that the transform from `map.lines`-coordinates to
IOAPI-coordinates is straightforward. So what I need to know, in order
to feed output from M3::get.map.lines.M3.proj(...) to
latticeExtra::layer(...), is how to set the coordinates of a
SpatialLines as produced by, e.g.,

> state.map.shp <-
>   maptools::map2SpatialLines(state.map, proj4string=lcc.crs)

or

> latticeExtra::layer(
>   sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))

or otherwise how to set the (x,y) coordinates being read by
latticeExtra::layer(...)

TIA, Tom Roche <Tom_Roche at pobox.com>