If someone else hasn't suggested it already, you will probably get more/better help on the R-sig-geo mailing list.
(if you decide to repost there, just mention up front that it's a repost and why) -Don -- Don MacQueen Lawrence Livermore National Laboratory 7000 East Ave., L-627 Livermore, CA 94550 925-423-1062 On 4/25/13 6:38 PM, "Tom Roche" <tom_ro...@pobox.com> wrote: > >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/95484c5d635 >02ab146402cedc3612dcdaf629bd7/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/95484c5d635 >02ab146402cedc3612dcdaf629bd7/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-ma >>p-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/95484c5d635 >02ab146402cedc3612dcdaf629bd7/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/95484c5d635 >02ab146402cedc3612dcdaf629bd7/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_ro...@pobox.com> <roche....@epa.gov> > >______________________________________________ >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. ______________________________________________ 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.