In JGCRI/gcammaptools: Plot GCAM Output on Maps

GCAM Mapping Tools Examples

Introduction

This vignette explains how to use the GCAM Mapping Tools package to display GCAM data in map form. Using this package has a couple of advantages over rolling your own maps. First, it's easy. By the time you've worked through these examples, you will be able to make maps of GCAM data with just a few simple commands. Second, we have defined some default projections and extents for common use cases. This will ensure that the maps you make look professional and conform to the GCAM house style.

if(!requireNamespace('rgcam', quietly=TRUE)) {
    stop('The rgcam package is required to build this vignette')
}

Setup

To get started with the GCAM map tools, attach the gcammaptools package.

library('gcammaptools')

Loading GCAM data

You will need to load your own GCAM data to work with. The best way to import GCAM results is to use the rgcam package, which will create a project file with your data in it. You can then load the project data with rgcam::loadProject() and retrieve tables for individual queries using rgcam::getQuery(). However, if you aren't using rgcam, you can start with any data frame that has a region column and one or more data columns.

library('rgcam')
### Load the example scenario data.
prj <- loadProject(system.file('sample-gcam-data','gcam-longform-sample.dat',                                              package='gcammaptools'))
listScenarios(prj)
listQueries(prj, 'Reference')

The sample data has one scenario with 19 queries. For demonstration purposes we will load the CO₂ emissions table. After loading the query you want to display, you must add the region identifiers used in the map data to the data frame using the add_region_ID() function.

co2 <- rgcam::getQuery(prj, 'CO2 emissions by region', 'Reference')
co2 <- dplyr::filter(co2, year==2050)
co2 <- add_region_ID(co2, lookupfile=rgn32, drops=rgn32)

### Show the first observations
head(co2)

The data is in long form, so using dplyr::filter() helps to get just the observations for the year you want to plot. At the end of this, co2 is a data frame that has CO₂ emissions by region for the year 2050. This structure can now be passed to plot_GCAM() to plot maps, as shown in the examples below.

Loading map data

When using plot_GCAM() for your GCAM data, you also need to provide the map data that it should be associated with. You should use one of the map datasets included in the package whenever possible, because they correspond to the geographical units reported by GCAM. The datasets available are map.rgn32, map.rgn14, map.basin235, map.usa, and map.chn. They correspond to the 32-region and 14-region GCAM region maps, the 235 global water basin map, and the 32-region map plus states/provinces for the USA and China. Additionally, the non-provincial datasets each have a simplified version such as map.rgn32.simple, which is helpful for faster plotting and less cluttered maps.

plot_GCAM(map.rgn32, title="Full 32-Region map")
plot_GCAM(map.rgn32.simple, title="Simplified 32-Region map")

It is possible, however, for you to load your own map data as well. The map data can be loaded by gcammaptools given the file path. Your map data can be in any of the following formats: sf objects, spatial data frames, ESRI Shapefiles, or GeoJSON files. You can either pass the file path to plot_GCAM() or load the data yourself first.

mapdata <- system.file("extdata/rgn32", "reg32_spart.shp", package = "gcammaptools")
mapdata
plot_GCAM(mapdata)

Sample maps

Example 1: Eckert III World Projection, Colored by Region

This example just plots the map data frame with the GCAM region name. Applying a scale with gcam32_colors, the default GCAM colors, you get each region colored according to a discrete color palette. For older data that uses 14-region GCAM, you can use the palette gcam14_colors.

plot_GCAM(map.rgn32.simple, col = 'region_name', proj = eck3) +
    ggplot2::scale_fill_manual(values = gcam32_colors, na.value=gray(0.75))

Example 2: Robinson World Projection, Colored by Regional CO₂ Emissions

In this example we plot the CO₂ data frame that we created above. We select the column value, which is the name of the column that contains the data. It is also necessary to specify how to join the this data frame to the map data, which is why we created the column id with the function add_region_id() above.

plot_GCAM(map.rgn32.simple, col='value', proj=robin, title="Robinson World", legend=T,
          gcam_df=co2) + 
    ggplot2::scale_fill_gradientn(colors = c("white", "red"),
                                  na.value = gray(0.75), name="CO2 Emissions (MTC)")

Example 3: U.S. Projection (Albers Equal Area)

This map is specialized to the continental USA. The na_aea and EXTENT_USA symbols are defined for convenience, but you can use any valid proj4 string (see proj4::project for how these strings are constructed) for the projection. The extent should be the bounding box of the plot area in the form c(lon.min, lon.max, lat.min, lat.max).

plot_GCAM(map.rgn32, col='region_id', proj=na_aea, extent=EXTENT_USA, 
          title="USA Albers Equal-Area")

Example 4: Africa Projection (Orthographic)

For superregions with a long north-south extent, the orthographic projection gives the best result. We have predefined one for Africa, but you can find others defined by proj4 strings from http://spatialreference.org as shown in the following example. Although the extent parameter is the best way to specify the view you want, you can fine-tune the final plot by adjusting the zoom.

plot_GCAM(map.rgn32, col='region_name', proj=af_ortho, extent=EXTENT_AFRICA,
          title="Africa Orthographic") +
    ggplot2::scale_fill_manual(values = gcam32_colors)

Example 5: Latin America Projection (Orthographic)

Orthographic projection of the Latin America superregion. Notice that projection strings can also be defined by specifying an EPSG, ESRI, or SR-ORG projection code.

plot_GCAM(map.rgn32, col='region_name', proj=7567, proj_type='SR-ORG',
          extent=EXTENT_LA, title="Latin America Orthographic") +
    ggplot2::scale_fill_manual(values = gcam32_colors)

Example 6: China Projection (Albers Equal Area)

A map of China. Although the projection is once again the Albers equal area projection, we have to have a different projection string because the string includes some information about the parallels the projection is based on.

plot_GCAM(map.rgn32, proj=ch_aea, extent=EXTENT_CHINA, title="China Albers Equal-Area")

Example 7: Global Water Basins

A map of the 235 global water basins. Because col is set to basin_name, each basin name is treated as a separate category, giving the colorful map below.

plot_GCAM(map.basin235, col='basin_name', proj=eck3)

Example 8: Gridded Data

The plot_GCAM_grid() function tiles gridded data over a base map.

co2grid <- rgcam::getQuery(prj, 'Cooling Degree Days', 'Reference')
plot_GCAM_grid(co2grid, col='value', proj=robin, extent=EXTENT_LA, legend=T) + 
    ggplot2::scale_fill_gradientn(colors=c("white", "red"), 
                                  guide=ggplot2::guide_colorbar(title="Cooling Degree Days", 
                                                                title.position="top"))

Example 9: Faceting

Using ggplot::facet_wrap() is often a good way to compare the results from several different scenarios. To do this with plot_GCAM(), it is important to make sure that your dataset contains the facet variable for each region you want to plot. By default, add_region_ID() puts in NA values for missing regions. The disaggregate parameter allows you to specify that the missing regions with no data should be included in each facet.

waterdata <- read.csv(system.file('extdata','facet_example.csv', package='gcammaptools'),
                      stringsAsFactors = F)
waterdata <- add_region_ID(waterdata, disaggregate = 'scenario')

plot_GCAM(map.rgn32.simple, col='value', legend=T, gcam_df=waterdata,
          title="Water Scarcity in 2100") + 
    ggplot2::scale_fill_gradientn(colors=c("deepskyblue", "firebrick1"),
                                  na.value=gray(0.9),
                                  guide=ggplot2::guide_colorbar(title="km^3",
                                                                barwidth=12)) +
    ggplot2::facet_wrap(~scenario)

Example 10: GCAM China

Data from GCAM China may contain extra regions that you want to drop as well as region name abbreviations that you want to translate into the full province name.

gcamchina <- read.csv(system.file('extdata', 'china_example.csv', package='gcammaptools'),
                      stringsAsFactors = F)
gcamchina <- add_region_ID(gcamchina, lookupfile = chn, provincefile = chn, drops = chn)

# Highlight the provinces by filtering out values from other regions
gcamchina[!is.na(gcamchina$id) & gcamchina$id <= 32, 'value'] <- NA
plot_GCAM(map.chn, col='value', gcam_df=gcamchina, proj=ch_aea,
          extent=EXTENT_CHINA)

Example 11: GCAM USA

Support is also provided for GCAM USA data. The map expects states to be identified by their two letter abbreviations, which can be added on using add_region_id() with lookupfile = usa, as long as the state full names are present.

gcamusa <- read.csv(system.file('extdata', 'usa_example.csv', package='gcammaptools'),
                    stringsAsFactors = F)
plot_GCAM(map.usa, col='value', gcam_df=gcamusa, gcam_key='region', proj=na_aea,
          extent=EXTENT_USA)

JGCRI/gcammaptools documentation built on May 7, 2019, 10:11 a.m.