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

[Oscar Perpiñan]

Hello,

Yesterday I read your question in stackoverflow. I tried to display
the ozone file with the state boundaries overlaid without success. In
my opinion, the problem is that the extent of the RasterLayer as read
by raster is quite different from the state boundaries. If I am not
wrong, this is the same problem you asked and answered in another
stackoverflow question:

http://stackoverflow.com/questions/13005181/how-to-set-the-display-bounds-of-a-rasterplot

I download the file from
https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ozone_lcc.nc
and then:

library(raster)
library(maps)
library(maptools)
lcc.crs <- "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97
+a=6370000 +b=6370000"
o3.raster <- raster('ozone_lcc.nc', varname="O3", crs=lcc.crs, band=1)
> 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=longlat +datum=WGS84
data source : /home/oscar/temp/ozone_lcc.nc
names       : O3
z-value     : 1
zvar        : O3
level       : 1 o3.raster

state.map <- map(database="state", plot=FALSE, projection='lambert',
parameters=c(40, 33))
state.map.shp <- map2SpatialLines(state.map, proj4string=CRS(lcc.crs))

> state.map.shp
class       : SpatialLines
nfeatures   : 169
extent      : -0.2226344, 0.2112617, -0.9123794, -0.6458026  (xmin,
xmax, ymin, ymax)
coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97
+a=6370000 +b=6370000


Best,

Oscar.

Grupo de Sistemas Fotovoltaicos (IES-UPM)
Dpto. Ingeniería Eléctrica (EUITI-UPM)
URL: http://procomun.wordpress.com
Twitter: @oscarperpinan


2013/2/15 Tom Roche <Tom_Roche at pobox.com>:
>
> Please advise me regarding overlay of an LCC-projected map on `raster`
> data plotted by `rasterVis::layerplot`. What I mean, why I ask:
>
> As discussed in detail @
>
> http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot
>
> (more detail than is feasible for an email post) I'm trying to overlay
> LCC-projected data from a netCDF file with an LCC-projected map. The
> netCDF uses the IOAPI convention (described in link), and also manages
> to hit yet again the `raster` filename-extension problem. (Dunno why,
> but that's twice in one week for me--bad karma?) The data's horizontal
> extents are CONUS, so I'm looking for an LCC CONUS map.
>
> In one example (in the link above), I use CRAN package=M3 to get
> (basically) a matrix of coordinates
>
> map.lines <- M3::get.map.lines.M3.proj(
>   file=lcc.file, database="state", units="m")
>
> and feed that to old-style graphics. But I'd much prefer to make this
> work like the other example, which resembles other code I have that is
> using `rasterVis::layerplot`. Unfortunately that other code is using
> unprojected/lon-lat data, and I must also handle projected (probably all
> LCC) data.
>
> 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