[R] [newbie] how to find and combine geographic maps with particular features?

[Tom Roche]

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>