[R] [lattice] display a projected map on a layerplot
Felix Andrews
felix at nfrac.org
Thu Feb 14 13:38:32 CET 2013
Hi
I'm a bit lost in the labyrinthine world of map projections, but
clearly something is mismatched between the two coordinate systems
being plotted:
1. the raster - the plot scale here seems to be cell numbers starting
from cell (1,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
2. the state boundaries - presumably projected correctly as lambert -
coordinate values less than 1
> summary(do.call("rbind", unlist(coordinates(state.map.shp), recursive = FALSE)))
V1 V2
Min. :-0.234093 Min. :-0.9169
1st Qu.:-0.000333 1st Qu.:-0.8289
Median : 0.080378 Median :-0.7660
Mean : 0.058492 Mean :-0.7711
3rd Qu.: 0.162993 3rd Qu.:-0.7116
Max. : 0.221294 Max. :-0.6343
Regards
Felix
On 14 February 2013 10:32, Tom Roche <Tom_Roche at pobox.com> wrote:
>
> 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 at 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 at 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_Roche at pobox.com> <Roche.Tom at 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 at 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.
--
Felix Andrews / 安福立
http://www.neurofractal.org/felix/
More information about the R-help
mailing list