Re: [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_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. : 283762Min. : 160844
# # 1st Qu.:26502441st Qu.:1054047
# # Median :3469204Median :1701052
# # Mean :3245997Mean :1643356
# # 3rd Qu.:43009693rd Qu.:2252531
# # Max. :4878260Max. :2993778
# # NA's :168NA's :168
# # end debugging
# Note above is not zero-centered, like our extents:
# extent : -2556000, 2952000, -1728000, 186 (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(
[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. : 283762Min. : 160844
# # 1st Qu.:26502441st Qu.:1054047
# # Median :3469204Median :1701052
# # Mean :3245997Mean :1643356
# # 3rd Qu.:43009693rd Qu.:2252531
# # Max. :4878260Max. :2993778
# # NA's :168NA's :168
# # end debugging
# Note above is not zero-centered, like our extents:
# extent : -2556000, 2952000, -1728000, 186 (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 +
