[R] [newbie] how to find and combine geographic maps with particular features?
Tom Roche
Tom_Roche at pobox.com
Fri Apr 26 03:38:35 CEST 2013
SUMMARY:
Specific problem: I'm regridding biomass-burning emissions from a
global/unprojected inventory to a regional projection (LCC over North
America). I need to have boundaries for Canada, Mexico, and US
(including US states), but also Caribbean and Atlantic nations
(notably the Bahamas). I would also like to add Canadian provinces and
Mexican states. How to put these together?
General problem: are there references regarding
* sources for different geographical and political features?
* combining maps for the different R graphics packages?
DETAILS:
(Apologies if this is a FAQ, but googling has not helped me with this.)
I'd appreciate help with a specific problem, as well as guidance
(e.g., pointers to docs) regarding the larger topic of combining
geographical maps (especially projected ones, i.e., not just lon-lat)
on plots of regional data (i.e., data that is multinational but not
global).
My specific problem is
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/downloads/GFED-3.1_2008_N2O_monthly_emissions_regrid_20130404_1344.pdf
which plots N2O concentrations from a global inventory of fire
emissions (GFED) regridded to a North American projection. (See
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na
for details.) The plot currently includes boundaries for Canada,
Mexico, and US (including US states, since this is being done for a US
agency), which are being gotten calling code from package=M3
http://cran.r-project.org/web/packages/M3/
like
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/95484c5d63502ab146402cedc3612dcdaf629bd7/vis_regrid_vis.r?at=master
> ## get projected North American map
> NorAm.shp <- project.NorAm.boundaries.for.CMAQ(
> units='m',
> extents.fp=template_input_fp,
> extents=template.extents,
> LCC.parallels=c(33,45),
> CRS=out.crs)
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/95484c5d63502ab146402cedc3612dcdaf629bd7/visualization.r?at=master
> # database: Geographical database to use. Choices include "state"
> # (default), "world", "worldHires", "canusamex", etc. Use
> # "canusamex" to get the national boundaries of the Canada, the
> # USA, and Mexico, along with the boundaries of the states.
> # The other choices ("state", "world", etc.) are the names of
> # databases included with the ‘maps’ and ‘mapdata’ packages.
> project.M3.boundaries.for.CMAQ <- function(
> database='state', # see `?M3::get.map.lines.M3.proj`
> units='m', # or 'km': see `?M3::get.map.lines.M3.proj`
> extents.fp, # path to extents file
> extents, # raster::extent object
> LCC.parallels=c(33,45), # LCC standard parallels: see https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain#wiki-EPA
> CRS # see `sp::CRS`
> ) {
> library(M3)
> ## Will replace raw LCC map's coordinates with:
> metadata.coords.IOAPI.list <- M3::get.grid.info.M3(extents.fp)
> metadata.coords.IOAPI.x.orig <- metadata.coords.IOAPI.list$x.orig
> metadata.coords.IOAPI.y.orig <- metadata.coords.IOAPI.list$y.orig
> metadata.coords.IOAPI.x.cell.width <- metadata.coords.IOAPI.list$x.cell.width
> metadata.coords.IOAPI.y.cell.width <- metadata.coords.IOAPI.list$y.cell.width
> library(maps)
> map.lines <- M3::get.map.lines.M3.proj(
> file=extents.fp, database=database, units="m")
> # dimensions are in meters, not cells. TODO: take argument
> map.lines.coords.IOAPI.x <-
> (map.lines$coords[,1] - metadata.coords.IOAPI.x.orig)
> map.lines.coords.IOAPI.y <-
> (map.lines$coords[,2] - metadata.coords.IOAPI.y.orig)
> map.lines.coords.IOAPI <-
> cbind(map.lines.coords.IOAPI.x, map.lines.coords.IOAPI.y)
> # # start debugging
> # class(map.lines.coords.IOAPI)
> # # [1] "matrix"
> # summary(map.lines.coords.IOAPI)
> # # map.lines.coords.IOAPI.x map.lines.coords.IOAPI.y
> # # Min. : 283762 Min. : 160844
> # # 1st Qu.:2650244 1st Qu.:1054047
> # # Median :3469204 Median :1701052
> # # Mean :3245997 Mean :1643356
> # # 3rd Qu.:4300969 3rd Qu.:2252531
> # # Max. :4878260 Max. :2993778
> # # NA's :168 NA's :168
> # # end debugging
> # Note above is not zero-centered, like our extents:
> # extent : -2556000, 2952000, -1728000, 1860000 (xmin, xmax, ymin, ymax)
> # So gotta add (xmin, ymin) below.
> ## Get LCC state map
> # see http://stackoverflow.com/questions/14865507/how-to-display-a-projected-map-on-an-rlatticelayerplot
> map.IOAPI <- maps::map(
> database="state", projection="lambert", par=LCC.parallels, plot=FALSE)
> # parameters to lambert: ^^^^^^^^^^^^^^^^^
> # see mapproj::mapproject
> map.IOAPI$x <- map.lines.coords.IOAPI.x + extents.xmin
> map.IOAPI$y <- map.lines.coords.IOAPI.y + extents.ymin
> map.IOAPI$range <- c(
> min(map.IOAPI$x),
> max(map.IOAPI$x),
> min(map.IOAPI$y),
> max(map.IOAPI$y))
> map.IOAPI.shp <-
> maptools::map2SpatialLines(map.IOAPI, proj4string=CRS)
> return(map.IOAPI.shp)
> } # end project.M3.boundaries.for.CMAQ
The maps are then overlaid on data and plotted using package=rasterVis
http://cran.r-project.org/web/packages/rasterVis/
with code like
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/95484c5d63502ab146402cedc3612dcdaf629bd7/vis_regrid_vis.r?at=master
> visualize.layers(
> nc.fp=monthly_out_fp,
> datavar.name=monthly_out_datavar_name,
> layer.dim.name=monthly_out_time_dim_name,
> sigdigs=sigdigs,
> brick=monthly.out.raster,
> map.list <- list(NorAm.shp),
> pdf.fp=monthly_out_pdf_fp,
> pdf.height=monthly_out_pdf_height,
> pdf.width=monthly_out_pdf_width
> )
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/src/95484c5d63502ab146402cedc3612dcdaf629bd7/visualization.r?at=master
> visualize.layers <- function(
> nc.fp, # path to netCDF datasource ...
> datavar.name, # name of the netCDF data variable
> layer.dim.name, # name of the datavar dimension indexing the layers (e.g., 'time')
> sigdigs=3, # precision to use for stats
> brick, # ... for data (a RasterBrick)
> map.list, # maps to overlay on data: list of SpatialLines or objects castable thereto
> pdf.fp, # path to PDF output
> pdf.height,
> pdf.width
> ) {
> nonplot.vis.layers(nc.fp, datavar.name, layer.dim.name, sigdigs)
> plot.vis.layers(brick, map.list, pdf.fp, pdf.height, pdf.width)
> } # end visualize.layers
...
> plot.vis.layers <- function(
> brick, # ... for data (a RasterBrick)
> map.list, # maps to overlay on data: list of SpatialLines or objects castable thereto
> pdf.fp, # path to PDF output
> pdf.height,
> pdf.width
> ) {
> # start plot driver
> cat(sprintf(
> '%s: plotting %s (may take awhile! TODO: add progress control)\n',
> 'plot.vis.layers', pdf.fp))
> pdf(file=pdf.fp, width=pdf.width, height=pdf.height)
> library(rasterVis)
> # I want to overlay the data with each map in the list, iteratively.
> # But the following does not work: only the last map in the list shows.
> # plots <- rasterVis::levelplot(brick,
> # # layers,
> # margin=FALSE,
> # # names.attr=names(global.df),
> # layout=c(1,length(names(brick))))
> # for (i.map in 1:length(map.list)) {
> # plots <- plots + latticeExtra::layer(
> # # why does this fail if map.shp is local? see 'KLUDGE' in callers :-(
> # sp::sp.lines(
> # as(map.list[[i.map]], "SpatialLines"),
> # lwd=0.8, col='darkgray'))
> # } # end for (i.map in 1:length(map.list))
> # plot(plots)
> # For now, kludge :-( handle lists of length 1 and 2 separately:
> if (length(map.list) == 1) {
> plot(
> rasterVis::levelplot(brick,
> margin=FALSE,
> layout=c(1,length(names(brick)))
> ) + latticeExtra::layer(
> sp::sp.lines(
> as(map.list[[1]], "SpatialLines"),
> lwd=0.8, col='darkgray')
> )
> )
> } else if (length(map.list) == 2) {
> plot(
> rasterVis::levelplot(brick,
> margin=FALSE,
> layout=c(1,length(names(brick)))
> ) + latticeExtra::layer(
> sp::sp.lines(
> as(map.list[[1]], "SpatialLines"),
> lwd=0.8, col='darkgray')
> ) + latticeExtra::layer(
> sp::sp.lines(
> as(map.list[[2]], "SpatialLines"),
> lwd=0.8, col='darkgray')
> )
> )
> } else {
> stop(sprintf('plot.vis.layers: ERROR: cannot handle (length(map.list)==%i', length(map.list)))
> }
> dev.off()
> } # end plot.vis.layers
This allows me to combine US state boundaries with national boundaries
of Canada and Mexico (though only via the kludge above). But the plot
https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/downloads/GFED-3.1_2008_N2O_monthly_emissions_regrid_20130404_1344.pdf
1. plainly also needs the Bahamas (see months Mar-Jun)
2. would benefit from Canadian provincial boundaries
3. ... and Mexican state boundaries
So I'd like to know specifically how to add those features. But I'd
also like to learn more about the general problems:
* How to combine maps when plotting using the different R graphics
packages (e.g., base, lattice, ggplot2)?
* How to find sources (i.e., R packages, map identifiers) for
arbitrary geographical and political features?
TIA, Tom Roche <Tom_Roche at pobox.com> <Roche.Tom at epa.gov>
More information about the R-help
mailing list