[R] [newbie] how to find and combine geographic maps with particular features?
MacQueen, Don
macqueen1 at llnl.gov
Fri Apr 26 20:53:15 CEST 2013
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.
More information about the R-help
mailing list