summary: I can display a lon-lat map on a lattice::layerplot, and I
can display a Lambert conformal conic (LCC) map on a spam::image, but
I can't display an LCC map on a lattice::layerplot. Example follows.
What am I doing wrong?
details:
I've been using `lattice` (via `rasterVis`) successfully to display
global atmospheric data, which works well enough (though I am
definitely intrigued by ggplot/ggmap). Particularly, I am able to
overlay my plots with maps, which is essential for the sort of work
I'm doing. However I'm currently unable to use lattice::layerplot for
some regional data projected LCC: the data plots, but I cannot get a
map to overlay. I can do this by cruder means, but would prefer to
know how to do this in lattice/rasterVis or similar (or ggplot/ggmap).
Two nearly-self-contained examples follow.
The data, ozone_lcc.nc, comes with CRAN package=M3 ... except that M3
supplies it as
system.file("extdata/ozone_lcc.ncf", package="M3")
and that extension currently causes problems for CRAN package=raster.
You can either rename that file (and put it some current working
directory), or you can download (270 kB)
https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ozone_lcc.nc
You can then run the following code:
##### start example #####
library("M3") # http://cran.r-project.org/web/packages/M3/
library("rasterVis") # http://cran.r-project.org/web/packages/rasterVis/
## Use an example file with projection=Lambert conformal conic.
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3")
lcc.file <- "./ozone_lcc.nc" # unfortunate problem with raster::raster
lcc.proj4 <- get.proj.info.M3(lcc.file)
lcc.proj4 # debugging
# [1] "+proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000"
# Note +lat_0=40 +lat_1=33 +lat_2=45 for maps::map@projection (below)
lcc.crs <- sp::CRS(lcc.proj4)
lcc.pdf <- "./ozone_lcc.pdf" # for output
## Read in variable with daily ozone.
# o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs)
# ozone_lcc.nc has 4 timesteps, choose 1 at random
o3.raster <- raster::raster(x=lcc.file, varname="O3", crs=lcc.crs, level=1)
o3.raster@crs <- lcc.crs # why does the above assignment not take?
o3.raster # debugging
# 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
## Get US state boundaries in projection units.
state.map <- maps::map(
database="state", projection="lambert", par=c(33,45), plot=FALSE)
# parameters to lambert: ^^^^^^^^^^^^
# see mapproj::mapproject
state.map.shp <-
maptools::map2SpatialLines(state.map, proj4string=lcc.crs)
pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
dev.off()
# change this as needed to view PDFs on your system
system(sprintf("xpdf %s", lcc.pdf))
# data looks good, but there's no map.
## Try again, with lambert(40,33)
state.map <- maps::map(
database="state", projection="lambert", par=c(40,33), plot=FALSE)
# parameters to lambert: ^^^^^^^^^^^^
# see mapproj::mapproject
state.map.shp <-
maptools::map2SpatialLines(state.map, proj4string=lcc.crs)
pdf(file=lcc.pdf)
rasterVis::levelplot(o3.raster, margin=FALSE
) + latticeExtra::layer(
sp::sp.lines(state.map.shp, lwd=0.8, col='darkgray'))
# still no map :-(
dev.off()
system(sprintf("xpdf %s", lcc.pdf))
##### end example #####
The data looks right, in that it greatly resembles the original
example (from which the above is adapted), which follows my .sig.
However the original-example code *does* draw a map, while the code
above does not. How to fix?
TIA, Tom Roche <tom_ro...@pobox.com> <roche....@epa.gov>
---------------original example follows to end of post---------------
# Following adapted from what is installed in my
# .../R/x86_64-pc-linux-gnu-library/2.14/m3AqfigExampleScript.r
# (probably by my sysadmin), which also greatly resembles
# https://wiki.epa.gov/amad/index.php/R_packages_developed_in_AMAD
# which is behind a firewall :-(
## EXAMPLE WITH LAMBERT CONIC CONFORMAL FILE.
library("M3")
library("aqfig") # http://cran.r-project.org/web/packages/aqfig/
## Use an example file with LCC projection: either local download ...
lcc.file <- "./ozone_lcc.nc"
## ... or as installed by package=M3:
# lcc.file <- system.file("extdata/ozone_lcc.ncf", package="M3")
## Choose the one that works for you.
## READ AND PLOT OZONE FROM FILE WITH LCC PROJECTION.
## Read in variable with daily ozone. Note that we can give dates
## rather than date-times, and that will include all time steps
## anytime during those days. Or, we can give lower and upper bounds
## and all time steps between these will be taken.
o3 <- M3::get.M3.var(
file=lcc.file, var="O3",
ldatetime=as.Date("2001-07-01"), udatetime=as.Date("2001-07-04"))
## Get colors and map boundaries for plot.
library("fields")
col.rng <- tim.colors(20)
detach("package:fields")
## Get state boundaries in projetion units.
map.lines <- M3::get.map.lines.M3.proj(file=lcc.file, database="state")
## Set color boundaries so that they incorporate the complete data range.
z.rng <- range(as.vector(o3$data))
## Make a color tile plot of the ozone for the first day (2001-07-01).
image(o3$x.cell.ctr, o3$y.cell.ctr, o3$data[,,1,1],
col=col.rng, zlim=z.rng,
xlab="Projection x-coord (km)", ylab="Projection y-coord (km)")
## Put date-time string and chemical name (O3) into a format I can use
## to label the actual figure.
date.str <- format(o3$datetime[1], "%Y-%m-%d")
title(main=bquote(paste(O[3], " on ", .(date.str), sep="")))
## Put the state boundaries on the plot.
lines(map.lines$coords)
## Add color bar to the right of plot to show what colors go with what
## ozone concentraton. NOTE: AFTER USING vertical.image.legend(), YOU
## WON'T BE ABLE TO ADD NEW FEATURES TO THE PLOT. Before making a new
## plot, you should open a new device or turn this device off.
vertical.image.legend(zlim=z.rng, col=col.rng)
dev.off() # close the plot if you haven't already
______________________________________________
R-help@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.
