Tom Roche
2013-Apr-26 01:38 UTC
[R] [newbie] how to find and combine geographic maps with particular features?
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.CMAQThe 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.layersThis 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>
MacQueen, Don
2013-Apr-26 18:53 UTC
[R] [newbie] how to find and combine geographic maps with particular features?
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_Roche at 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_Roche at pobox.com> <Roche.Tom at epa.gov> > >______________________________________________ >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.