Drawing beautiful maps programmatically with R, sf and ggplot2 — Part 1： Basics

转载自：http://feedproxy.google.com/~r/RBloggers/~3/n7DiModPhCo/

Mel Moreno and Mathieu Basille
发表于

view raw Rmd

Maps are used in a variety of fields to express data in an appealing andinterpretive way. Data can be expressed into simplified patterns, andthis data interpretation is generally lost if the data is only seenthrough a spread sheet. Maps can add vital context by incorporating manyvariables into an easy to read and applicable context. Maps are alsovery important in the information world because they can quickly allowthe public to gain better insight so that they can stay informed. It’scritical to have maps be effective, which means creating maps that canbe easily understood by a given audience. For instance, maps that needto be understood by children would be very different from maps intendedto be shown to geographers.

Knowing what elements are required to enhance your data is key intomaking effective maps. Basic elements of a map that should be consideredare polygon, points, lines, and text. Polygons, on a map, are closedshapes such as country borders. Lines are considered to be linear shapesthat are not filled with any aspect, such as highways, streams, orroads. Finally, points are used to specify specific positions, such ascity or landmark locations. With that in mind, one need to think aboutwhat elements are required in the map to really make an impact, andconvey the information for the intended audience. Layout and formattingare the second critical aspect to enhance data visually. The arrangementof these map elements and how they will be drawn can be adjusted to makea maximum impact.

A solution using R and its ecosystem of packages

Current solutions for creating maps usually involves GIS software, suchas ArcGIS, QGIS, eSpatial, etc., which allow to visually prepare a map,in the same approach as one would prepare a poster or a document layout.On the other hand, R, a free and open-source software developmentenvironment (IDE) that is used for computing statistical data andgraphic in a programmable language, has developed advanced spatialcapabilities over the years, and can be used to draw mapsprogrammatically.

R is a powerful and flexible tool. R can be used from calculating datasets to creating graphs and maps with the same data set. R is also free,which makes it easily accessible to anyone. Some other advantages ofusing R is that it has an interactive language, data structures,graphics availability, a developed community, and the advantage ofadding more functionalities through an entire ecosystem of packages. Ris a scriptable language that allows the user to write out a code inwhich it will execute the commands specified.

Using R to create maps brings these benefits to mapping. Elements of amap can be added or removed with ease — R code can be tweaked to makemajor enhancements with a stroke of a key. It is also easy to reproducethe same maps for different data sets. It is important to be able toscript the elements of a map, so that it can be re-used and interpretedby any user. In essence, comparing typical GIS software and R fordrawing maps is similar to comparing word processing software (e.g.Microsoft Office or LibreOffice) and a programmatic typesetting systemsuch as LaTeX, in that typical GIS software implement a WYSIWIG approach(“What You See Is What You Get”), while R implements a WYSIWYM approach(“What You See Is What You Mean”).

The package ggplot2 implements the grammar of graphics in R, as a wayto create code that make sense to the user: The grammar of graphics is aterm used to breaks up graphs into semantic components, such asgeometries and layers. Practically speaking, it allows (and forces!) theuser to focus on graph elements at a higher level of abstraction, andhow the data must be structured to achieve the expected outcome. Whileggplot2 is becoming the de facto standard for R graphs, it does nothandle spatial data specifically. The current state-of-the-art ofspatial objects in R relies on Spatial classes defined in the packagesp, but the new packagesf has recently implementedthe “simple feature” standard, and is steadily taking over sp.Recently, the package ggplot2 has allowed the use of simple featuresfrom the package sf as layers in a graph. The combination ofggplot2 and sf therefore enables to programmatically create maps,using the grammar of graphics, just as informative or visually appealingas traditional GIS software.

Drawing beautiful maps programmatically with R, sf and ggplot2

This tutorial is the first part in a series of three:

In this part, we will cover the fundamentals of mapping using ggplot2associated to sf, and presents the basics elements and parameters wecan play with to prepare a map.

Getting started

Many R packages are available from CRAN,the Comprehensive R Archive Network, which is the primary repository ofR packages. The full list of packages necessary for this series oftutorials can be installed with:

1
2
3
install.packages(c("cowplot", "googleway", "ggplot2", "ggrepel", 
 "ggspatial", "libwgeom", "sf", "rworldmap", "rworldxtra"))

We start by loading the basic packages necessary for all maps, i.e.ggplot2 and sf. We also suggest to use the classic dark-on-lighttheme for ggplot2 (theme_bw), which is appropriate for maps:

1
2
3
4
library("ggplot2")
theme_set(theme_bw())
library("sf")

The package rworldmap provides a map of countries of the entire world;a map with higher resolution is available in the package rworldxtra.We use the function getMap to extract the world map (the resolutioncan be set to "low", if preferred):

1
2
3
4
5
6
7
8
9
library("rworldmap")
library("rworldxtra")
world <- getMap(resolution = "high")
class(world)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

The world map is available as a SpatialPolygonsDataFrame from thepackage sp; we thus convert it to a simple feature using st_as_sffrom package sf:

1
2
3
4
5
world <- st_as_sf(world)
class(world)

## [1] "sf" "data.frame"

Data and basic plot (ggplot and geom_sf)

First, let us start with creating a base map of the world usingggplot2. This base map will then be extended with different mapelements, as well as zoomed in to an area of interest. We can check thatthe world map was properly retrieved and converted into an sf object,and plot it with ggplot2:

1
2
3
ggplot(data = world) +
 geom_sf()

This call nicely introduces the structure of a ggplot call: The firstpart ggplot(data = world) initiates the ggplot graph, and indicatesthat the main data is stored in the world object. The line ends upwith a + sign, which indicates that the call is not complete yet, andeach subsequent line correspond to another layer or scale. In this case,we use the geom_sf function, which simply adds a geometry stored in asf object. By default, all geometry functions use the main datadefined in ggplot(), but we will see later how to provide additionaldata.

Note that layers are added one at a time in a ggplot call, so theorder of each layer is very important. All data will have to be in ansf format to be used by ggplot2; data in other formats (e.g. classesfrom sp) will be manually converted to sf classes if necessary.

Title, subtitle, and axis labels (ggtitle, xlab, ylab)

A title and a subtitle can be added to the map using the functionggtitle, passing any valid character string (e.g. with quotationmarks) as arguments. Axis names are absent by default on a map, but canbe changed to something more suitable (e.g. “Longitude” and “Latitude”),depending on the map:

1
2
3
4
5
ggplot(data = world) +
 geom_sf() +
 xlab("Longitude") + ylab("Latitude") +
 ggtitle("World map", subtitle = paste0("(", length(unique(world$NAME)), " countries)"))

Map color (geom_sf)

In many ways, sf geometries are no different than regular geometries,and can be displayed with the same level of control on their attributes.Here is an example with the polygons of the countries filled with agreen color (argument fill), using black for the outline of thecountries (argument color):

1
2
3
ggplot(data = world) + 
 geom_sf(color = "black", fill = "lightgreen")

The package ggplot2 allows the use of more complex color schemes, suchas a gradient on one variable of the data. Here is another example thatshows the population of each country. In this example, we use the“viridis” colorblind-friendly palette for the color gradient (withoption = "plasma" for the plasma variant), using the square root ofthe population (which is stored in the variable POP_EST of the worldobject):

1
2
3
4
ggplot(data = world) +
 geom_sf(aes(fill = POP_EST)) +
 scale_fill_viridis_c(option = "plasma", trans = "sqrt")

Projection and extent (coord_sf)

The function coord_sf allows to deal with the coordinate system, whichincludes both projection and extent of the map. By default, the map willuse the coordinate system of the first layer that defines one (i.e.scanned in the order provided), or if none, fall back on WGS84(latitude/longitude, the reference system used in GPS). Using theargument crs, it is possible to override this setting, and project onthe fly to any projection. This can be achieved using any valid PROJ4string (here, the European-centric ETRS89 Lambert Azimuthal Equal-Areaprojection):

1
2
3
4
ggplot(data = world) +
 geom_sf() +
 coord_sf(crs = "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs ")

Spatial Reference System Identifier (SRID) or an European PetroleumSurvey Group (EPSG) code are available for the projection of interest,they can be used directly instead of the full PROJ4 string. The twofollowing calls are equivalent for the ETRS89 Lambert AzimuthalEqual-Area projection, which is EPSG code 3035:

1
2
3
4
5
6
7
8
ggplot(data = world) +
 geom_sf() +
 coord_sf(crs = "+init=epsg:3035")

ggplot(data = world) +
 geom_sf() +
 coord_sf(crs = st_crs(3035))

The extent of the map can also be set in coord_sf, in practiceallowing to “zoom” in the area of interest, provided by limits on thex-axis (xlim), and on the y-axis (ylim). Note that the limits areautomatically expanded by a fraction to ensure that data and axes don’toverlap; it can also be turned off to exactly match the limits providedwith expand = FALSE:

1
2
3
4
ggplot(data = world) +
 geom_sf() +
 coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE)

Scale bar and North arrow (package ggspatial)

Several packages are available to create a scale bar on a map (e.g.prettymapr, vcd, ggsn, or legendMap). We introduce here thepackage ggspatial, which provides easy-to-use functions…

scale_bar that allows to add simultaneously the north symbol and ascale bar into the ggplot map. Five arguments need to be set manually:lon, lat, distance_lon, distance_lat, and distance_legend. Thelocation of the scale bar has to be specified in longitude/latitude inthe lon and lat arguments. The shaded distance inside the scale baris controlled by the distance_lon argument. while its width isdetermined by distance_lat. Additionally, it is possible to change thefont size for the legend of the scale bar (argument legend_size, whichdefaults to 3). The North arrow behind the “N” north symbol can also beadjusted for its length (arrow_length), its distance to the scale(arrow_distance), or the size the N north symbol itself(arrow_north_size, which defaults to 6). Note that all distances(distance_lon, distance_lat, distance_legend, arrow_length,arrow_distance) are set to "km" by default in distance_unit; theycan also be set to nautical miles with “nm”, or miles with “mi”.

1
2
3
4
5
6
7
8
9
10
11
library("ggspatial")
ggplot(data = world) +
 geom_sf() +
 annotation_scale(location = "bl", width_hint = 0.5) +
 annotation_north_arrow(location = "bl", which_north = "true", 
 pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
 style = north_arrow_fancy_orienteering) +
 coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97))

## Scale on map varies by more than 10%, scale bar may be inaccurate

Note the warning of the inaccurate scale bar: since the map useunprojected data in longitude/latitude (WGS84) on an equidistantcylindrical projection (all meridians being parallel), length in(kilo)meters on the map directly depends mathematically on the degree oflatitude. Plots of small regions or projected data will often allow formore accurate scale bars.

Country names and other names (geom_text and annotate)

The world data set already contains country names and the coordinatesof the centroid of each country (among more information). We can usethis information to plot country names, using world as a regulardata.frame in ggplot2. We first check the country name information:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
head(world[, c("NAME", "LON", "LAT")])

## Simple feature collection with 6 features and 3 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -70.06164 ymin: -18.0314 xmax: 74.89231 ymax: 60.48075
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## NAME LON LAT geometry
## 1 Aruba -69.97345 12.51678 MULTIPOLYGON (((-69.87609 1...
## 2 Afghanistan 66.00845 33.83627 MULTIPOLYGON (((71.02458 38...
## 3 Angola 17.56405 -12.32934 MULTIPOLYGON (((13.98233 -5...
## 4 Anguilla -63.05667 18.22432 MULTIPOLYGON (((-63.0369 18...
## 5 Albania 20.05399 41.14258 MULTIPOLYGON (((20.06496 42...
## 6 Aland 19.94429 60.22851 MULTIPOLYGON (((19.91892 60...

The function geom_text can be used to add a layer of text to a mapusing geographic coordinates. The function requires the data needed toenter the country names, which is the same data as the world map. Again,we have a very flexible control to adjust the text at will on manyaspects:

  • The size (argument size);

  • The alignment, which is centered by default on the coordinatesprovided. The text can be adjusted horizontally or vertically usingthe arguments hjust and vjust, which can either be a numberbetween 0 (right/bottom) and 1 (top/left) or a character (“left”,“middle”, “right”, “bottom”, “center”, “top”). The text can also beoffset horizontally or vertically with the argument nudge_x andnudge_y;

  • The font of the text, for instance its color (argument color) orthe type of font (fontface);

  • The overlap of labels, using the argument check_overlap, whichremoves overlapping text. Alternatively, when there is a lot ofoverlapping labels, the package ggrepel provides ageom_text_repel function that moves label around so that they donot overlap.

Additionally, the annotate function can be used to add a singlecharacter string at a specific location, as demonstrated here to add theGulf of Mexico:

1
2
3
4
5
6
7
8
ggplot(data = world) +
 geom_sf() +
 geom_text(aes(LON, LAT, label = NAME), size = 4, hjust = "left", 
 color = "darkblue", fontface = "bold", check_overlap = TRUE) +
 annotate(geom = "text", x = -90, y = 26, label = "Gulf of Mexico", 
 fontface = "italic", color = "grey22", size = 6) +
 coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE)

Now to make the final touches, the theme of the map can be edited tomake it more appealing. We suggested the use of theme_bw for astandard theme, but there are many other themes that can be selectedfrom (see for instance ?ggtheme in ggplot2, or the packageggthemes which provideseveral useful themes). Moreover, specific theme elements can be tweakedto get to the final outcome:

  • Position of the legend: Although not used in this example, theargument legend.position allows to automatically place the legendat a specific location (e.g. "topright", "bottomleft", etc.);

  • Grid lines (graticules) on the map: by using panel.grid.major andpanel.grid.minor, grid lines can be adjusted. Here we set them toa gray color and dashed line type to clearly distinguish them fromcountry borders lines;

  • Map background: the argument panel.background can be used to colorthe background, which is the ocean essentially, with a light blue;

  • Many more elements of a theme can be adjusted, which would be toolong to cover here. We refer the reader to the documentation for thefunction theme.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
ggplot(data = world) +
 geom_sf(fill = "antiquewhite1") +
 geom_text(aes(LON, LAT, label = NAME), size = 4, hjust = "left",
 color = "darkblue", fontface = "bold", check_overlap = TRUE) +
 annotate(geom = "text", x = -90, y = 26, label = "Gulf of Mexico", 
 fontface = "italic", color = "grey22", size = 6) +
 annotation_scale(location = "bl", width_hint = 0.5) +
 annotation_north_arrow(location = "bl", which_north = "true", 
 pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
 style = north_arrow_fancy_orienteering) +
 coord_sf(xlim = c(-102.15, -74.12), ylim = c(7.65, 33.97), expand = FALSE) +
 xlab("Longitude") + ylab("Latitude") +
 ggtitle("Map of the Gulf of Mexico and the Caribbean Sea") +
 theme(panel.grid.major = element_line(color = gray(.5),
 linetype = "dashed", size = 0.5),
 panel.background = element_rect(fill = "aliceblue"))

Saving the map with ggsave

The final map now ready, it is very easy to save it using ggsave. Thisfunction allows a graphic (typically the last plot displayed) to besaved in a variety of formats, including the most common PNG (rasterbitmap) and PDF (vector graphics), with control over the size andresolution of the outcome. For instance here, we save a PDF version ofthe map, which keeps the best quality, and a PNG version of it for webpurposes:

1
2
3
ggsave("map.pdf")
ggsave("map_web.png", width = 6, height = 6, dpi = "screen")

Related