how to display a projected map on an rasterVis::layerplot?
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>