Subscribe to R-spatial feed
R Spatial software blogs and ideas 2018-10-30T18:47:18+00:00 Jekyll v3.8.4 https://avatars1.githubusercontent.com/u/25086656
Updated: 13 hours 59 min ago

Drawing beautiful maps programmatically with R, sf and ggplot2 — Part 3: Layouts

Thu, 10/25/2018 - 02:00

view raw Rmd

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

After the presentation of the basic map concepts, and the flexible approach in layer implemented in ggplot2, this part illustrates how to achieve complex layouts, for instance with map insets, or several maps combined. Depending on the visual information that needs to be displayed, maps and their corresponding data might need to be arranged to create easy to read graphical representations. This tutorial will provide different approaches to arranges maps in the plot, in order to make the information portrayed more aesthetically appealing, and most importantly, convey the information better.

Getting started

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

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-light theme for ggplot2 (theme_bw), which is more appropriate for maps:

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 resolution can be set to "low", if preferred):

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 the package sp; we thus convert it to a simple feature using st_as_sf from package sf:

world <- st_as_sf(world) class(world) ## [1] "sf" "data.frame" General concepts

There are 2 solutions to combine sub-maps:

  • using Grobs (graphic objects, allow plots only in plot region, based on coordinates), which directly use ggplot2
  • using ggdraw (allows plots anywhere, including outer margins, based on relative position) from package cowplot

Example illustrating the difference between the two, and their use:

(g1 <- qplot(0:10, 0:10))

(g1_void <- g1 + theme_void() + theme(panel.border = element_rect(colour = "black", fill = NA)))

Graphs from ggplot2 can be saved, like any other R object. They can then be reused later in other functions.

Using grobs, and annotation_custom:

g1 + annotation_custom( grob = ggplotGrob(g1_void), xmin = 0, xmax = 3, ymin = 5, ymax = 10 ) + annotation_custom( grob = ggplotGrob(g1_void), xmin = 5, xmax = 10, ymin = 0, ymax = 3 )

Using ggdraw (note: used to build on top of initial plot; could be left empty to arrange subplots on a grid; plots are “filled” with their plots, unless the plot itself has a constrained ratio, like a map):

ggdraw(g1) + draw_plot(g1_void, width = 0.25, height = 0.5, x = 0.02, y = 0.48) + draw_plot(g1_void, width = 0.5, height = 0.25, x = 0.75, y = 0.09)

Several maps side by side or on a grid

Having a way show in a visualization, a specific area can be very useful. Many scientists usually create maps for each specific area individually. This is fine, but there are simpler ways to display what is needed for a report, or publication.

This exmaple is using two maps side by side, including the legend of the first one. It illustrates how to use a custom grid, which can be made a lot more complex with different elements.

First, simplify REGION for the legend:

levels(world$REGION)[7] <- "South America"

Prepare the subplots, #1 world map:

(gworld <- ggplot(data = world) + geom_sf(aes(fill = REGION)) + geom_rect(xmin = -102.15, xmax = -74.12, ymin = 7.65, ymax = 33.97, fill = NA, colour = "black", size = 1.5) + scale_fill_viridis_d(option = "plasma") + theme(panel.background = element_rect(fill = "azure"), panel.border = element_rect(fill = NA)))

And #2 Gulf map :

(ggulf <- ggplot(data = world) + geom_sf(aes(fill = REGION)) + 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) + scale_fill_viridis_d(option = "plasma") + theme(legend.position = "none", axis.title.x = element_blank(), axis.title.y = element_blank(), panel.background = element_rect(fill = "azure"), panel.border = element_rect(fill = NA)))

The command ggplotGrob signals to ggplot to take each created map, and how to arrange each map. The argument coord_equal can specify the length, ylim, and width, xlim, for the entire plotting area. Where as in annotation_custom, each maps’ xmin, xmax, ymin, and ymax can be specified to allow for complete customization.

#Creating a faux empty data frame df <- data.frame() plot1<-ggplot(df) + geom_point() + xlim(0, 10) + ylim(0, 10) plot2<-ggplot(df) + geom_point() + xlim(0, 10) + ylim(0, 10) ggplot() + coord_equal(xlim = c(0, 3.3), ylim = c(0, 1), expand = FALSE) + annotation_custom(ggplotGrob(plot1), xmin = 0, xmax = 1.5, ymin = 0, ymax = 1) + annotation_custom(ggplotGrob(plot2), xmin = 1.5, xmax = 3, ymin = 0, ymax = 1) + theme_void()

Below is the final map, using the same methodology as the exmaple plot above. Using ggplot to arrange maps, allows for easy and quick plotting in one function of R code.

ggplot() + coord_equal(xlim = c(0, 3.3), ylim = c(0, 1), expand = FALSE) + annotation_custom(ggplotGrob(gworld), xmin = 0, xmax = 2.3, ymin = 0, ymax = 1) + annotation_custom(ggplotGrob(ggulf), xmin = 2.3, xmax = 3.3, ymin = 0, ymax = 1) + theme_void()

In the second approach, using cowplot::plot_grid to arrange ggplot figures, is quite versatile. Any ggplot figure can be arranged just like the figure above. There are many commands that allow for the map to have different placements, such as nrow=1 means that the figure will only occupy one row and multiple columns, and ncol=1 means the figure will be plotted on one column and multiple rows. The command rel_widths establishes the width of each map, meaning that the first map gworld will have a relative width of 2.3, and the map ggulf has the relative width of 1.

library("cowplot") theme_set(theme_bw()) plot_grid(gworld, ggulf, nrow = 1, rel_widths = c(2.3, 1))

Some other commands can adjust the position of the figures such as adding align=v to align vertically, and align=h to align horiztonally.

Note also the existence of get_legend (cowplot), and that the legend can be used as any object.

This map can be save using,ggsave:

ggsave("grid.pdf", width = 15, height = 5) Map insets

For map insets directly on the background map, both solutions are viable (and one might prefer one or the other depending on relative or absolute coordinates).

Map example using map of the 50 states of the US, including Alaska and Hawaii (note: not to scale for the latter), using reference projections for US maps. First map (continental states) use a 10/6 figure:

usa <- subset(world, ADMIN == "United States of America") ## US National Atlas Equal Area (2163) ## http://spatialreference.org/ref/epsg/us-national-atlas-equal-area/ (mainland <- ggplot(data = usa) + geom_sf(fill = "cornsilk") + coord_sf(crs = st_crs(2163), xlim = c(-2500000, 2500000), ylim = c(-2300000, 730000)))

Alaska map (note: datum = NA removes graticules and coordinates):

## Alaska: NAD83(NSRS2007) / Alaska Albers (3467) ## http://www.spatialreference.org/ref/epsg/3467/ (alaska <- ggplot(data = usa) + geom_sf(fill = "cornsilk") + coord_sf(crs = st_crs(3467), xlim = c(-2400000, 1600000), ylim = c(200000, 2500000), expand = FALSE, datum = NA))

Hawaii map:

## Hawaii: Old Hawaiian (4135) ## http://www.spatialreference.org/ref/epsg/4135/ (hawaii <- ggplot(data = usa) + geom_sf(fill = "cornsilk") + coord_sf(crs = st_crs(4135), xlim = c(-161, -154), ylim = c(18, 23), expand = FALSE, datum = NA))

Using ggdraw from cowplot (tricky to define exact positions; note the use of the ratios of the inset, combined with the ratio of the plot):

(ratioAlaska <- (2500000 - 200000) / (1600000 - (-2400000))) ## [1] 0.575 (ratioHawaii <- (23 - 18) / (-154 - (-161))) ## [1] 0.7142857 ggdraw(mainland) + draw_plot(alaska, width = 0.26, height = 0.26 * 10/6 * ratioAlaska, x = 0.05, y = 0.05) + draw_plot(hawaii, width = 0.15, height = 0.15 * 10/6 * ratioHawaii, x = 0.3, y = 0.05)

This plot can be saved using ggsave:

ggsave("map-us-ggdraw.pdf", width = 10, height = 6)

The same kind of plot can be created using grobs, with ggplotGrob, (note the use of xdiff/ydiff and arbitrary ratios):

mainland + annotation_custom( grob = ggplotGrob(alaska), xmin = -2750000, xmax = -2750000 + (1600000 - (-2400000))/2.5, ymin = -2450000, ymax = -2450000 + (2500000 - 200000)/2.5 ) + annotation_custom( grob = ggplotGrob(hawaii), xmin = -1250000, xmax = -1250000 + (-154 - (-161))*120000, ymin = -2450000, ymax = -2450000 + (23 - 18)*120000 )

This plot can be saved using ggsave:

ggsave("map-inset-grobs.pdf", width = 10, height = 6)

The print command can also be used place multiple maps in one plotting area.

To specify where each plot is displayed with the print function, the argument viewport needs to include the maximum width and height of each map, and the minimum x and y coordinates of where the maps are located in the plotting area. The argument just will make a position on how the secondary maps will be displayed. All maps are defaulted the same size, until the sizes are adjusted with width and height.

vp <- viewport(width = 0.37, height = 0.10, x = 0.20, y =0.25, just = c("bottom")) vp1<- viewport(width = 0.37, height = 0.10, x = 0.35, y =0.25, just = c("bottom"))

Theprint function uses the previous specifications that were listed in each plots’ respective viewport, with vp=.

print(mainland) print(alaska, vp=vp) print(hawaii, vp=vp1)

Several maps connected with arrows

To bring about a more lively map arrangement, arrows can be used to bring the viewers eyes to specific areas in the plot. The next example will create a map with zoomed in areas, pointed to by arrows.

Firstly, we will create our main map, and then our zoomed in areas.

Site coordinates, same as Tutorial #1:

sites <- st_as_sf(data.frame(longitude = c(-80.15, -80.1), latitude = c(26.5, 26.8)), coords = c("longitude", "latitude"), crs = 4326, agr = "constant")

Mainlaind map of Florida, #1:

(florida <- ggplot(data = world) + geom_sf(fill = "antiquewhite1") + geom_sf(data = sites, size = 4, shape = 23, fill = "darkred") + annotate(geom = "text", x = -85.5, y = 27.5, label = "Gulf of Mexico", color = "grey22", size = 4.5) + coord_sf(xlim = c(-87.35, -79.5), ylim = c(24.1, 30.8)) + xlab("Longitude")+ ylab("Latitude")+ theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue"), panel.border = element_rect(fill = NA)))

A map for site A is created by layering the map and points we created earlier. ggplot layers geom_sf objects and plot them spatially.

(siteA <- ggplot(data = world) + geom_sf(fill = "antiquewhite1") + geom_sf(data = sites, size = 4, shape = 23, fill = "darkred") + coord_sf(xlim = c(-80.25, -79.95), ylim = c(26.65, 26.95), expand = FALSE) + annotate("text", x = -80.18, y = 26.92, label= "Site A", size = 6) + theme_void() + theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue"), panel.border = element_rect(fill = NA)))

A map for site B:

(siteB <- ggplot(data = world) + geom_sf(fill = "antiquewhite1") + geom_sf(data = sites, size = 4, shape = 23, fill = "darkred") + coord_sf(xlim = c(-80.3, -80), ylim = c(26.35, 26.65), expand = FALSE) + annotate("text", x = -80.23, y = 26.62, label= "Site B", size = 6) + theme_void() + theme(panel.grid.major = element_line(colour = gray(0.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue"), panel.border = element_rect(fill = NA)))

Coordinates of the two arrows will need to be specified before plotting. The argumemnts x1, and x2 will plot the arrow line from a specific starting x-axis location,x1, and ending in a specific x-axis,x2. The same applies for y1 and y2, with the y-axis respectively:

arrowA <- data.frame(x1 = 18.5, x2 = 23, y1 = 9.5, y2 = 14.5) arrowB <- data.frame(x1 = 18.5, x2 = 23, y1 = 8.5, y2 = 6.5)

Final map using (ggplot only). The argument geom_segment, will be the coordinates created in the previous script, to plot line segments ending with an arrow using arrow=arrow():

ggplot() + coord_equal(xlim = c(0, 28), ylim = c(0, 20), expand = FALSE) + annotation_custom(ggplotGrob(florida), xmin = 0, xmax = 20, ymin = 0, ymax = 20) + annotation_custom(ggplotGrob(siteA), xmin = 20, xmax = 28, ymin = 11.25, ymax = 19) + annotation_custom(ggplotGrob(siteB), xmin = 20, xmax = 28, ymin = 2.5, ymax = 10.25) + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowA, arrow = arrow(), lineend = "round") + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowB, arrow = arrow(), lineend = "round") + theme_void()

This plot can be saved using ggsave:

ggsave("florida-sites.pdf", width = 10, height = 7)

ggdraw could also be used for a similar result, with the argument draw_plot:

ggdraw(xlim = c(0, 28), ylim = c(0, 20)) + draw_plot(florida, x = 0, y = 0, width = 20, height = 20) + draw_plot(siteA, x = 20, y = 11.25, width = 8, height = 8) + draw_plot(siteB, x = 20, y = 2.5, width = 8, height = 8) + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowA, arrow = arrow(), lineend = "round") + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), data = arrowB, arrow = arrow(), lineend = "round")

Drawing beautiful maps programmatically with R, sf and ggplot2 — Part 2: Layers

Thu, 10/25/2018 - 02:00

view raw Rmd

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

In the previous part, we presented general concepts with a map with little information (country borders only). The modular approach of ggplot2 allows to successively add additional layers, for instance study sites or administrative delineations, as will be illustrated in this part.

Getting started

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

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-light theme for ggplot2 (theme_bw), which is more appropriate for maps:

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 resolution can be set to "low", if preferred):

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 the package sp; we thus convert it to a simple feature using st_as_sf from package sf:

world <- st_as_sf(world) class(world) ## [1] "sf" "data.frame" Adding additional layers: an example with points and polygons Field sites (point data)

We start by defining two study sites, according to their longitude and latitude, stored in a regular data.frame:

(sites <- data.frame(longitude = c(-80.144005, -80.109), latitude = c(26.479005, 26.83))) ## longitude latitude ## 1 -80.14401 26.47901 ## 2 -80.10900 26.83000

The quickest way to add point coordinates is with the general-purpose function geom_point, which works on any X/Y coordinates, of regular data points (i.e. not geographic). As such, we can adjust all characteristics of points (e.g. color of the outline and the filling, shape, size, etc.), for all points, or using grouping from the data (i.e defining their “aesthetics”). In this example, we add the two points as diamonds (shape = 23), filled in dark red (fill = "darkred") and of bigger size (size = 4):

ggplot(data = world) + geom_sf() + geom_point(data = sites, aes(x = longitude, y = latitude), size = 4, shape = 23, fill = "darkred") + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

A better, more flexible alternative is to use the power of sf: Converting the data frame to a sf object allows to rely on sf to handle on the fly the coordinate system (both projection and extent), which can be very useful if the two objects (here world map, and sites) are not in the same projection. To achieve the same result, the projection (here WGS84, which is the CRS code #4326) has to be a priori defined in the sf object:

(sites <- st_as_sf(sites, coords = c("longitude", "latitude"), crs = 4326, agr = "constant")) ## Simple feature collection with 2 features and 0 fields ## geometry type: POINT ## dimension: XY ## bbox: xmin: -80.14401 ymin: 26.479 xmax: -80.109 ymax: 26.83 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## geometry ## 1 POINT (-80.14401 26.47901) ## 2 POINT (-80.109 26.83) ggplot(data = world) + geom_sf() + geom_sf(data = sites, size = 4, shape = 23, fill = "darkred") + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

Note that coord_sf has to be called after all geom_sf calls, as to supersede any former input.

States (polygon data)

It would be informative to add finer administrative information on top of the previous map, starting with state borders and names. The package maps (which is automatically installed and loaded with ggplot2) provides maps of the USA, with state and county borders, that can be retrieved and converted as sf objects:

library("maps") states <- st_as_sf(map("state", plot = FALSE, fill = TRUE)) head(states) ## Simple feature collection with 6 features and 1 field ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -124.3834 ymin: 30.24071 xmax: -71.78015 ymax: 42.04937 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## geometry ID ## 1 MULTIPOLYGON (((-87.46201 3... alabama ## 2 MULTIPOLYGON (((-114.6374 3... arizona ## 3 MULTIPOLYGON (((-94.05103 3... arkansas ## 4 MULTIPOLYGON (((-120.006 42... california ## 5 MULTIPOLYGON (((-102.0552 4... colorado ## 6 MULTIPOLYGON (((-73.49902 4... connecticut

State names are part of this data, as the ID variable. A simple (but not necessarily optimal) way to add state name is to compute the centroid of each state polygon as the coordinates where to draw their names. Centroids are computed with the function st_centroid, their coordinates extracted with st_coordinates, both from the package sf, and attached to the state object:

states <- cbind(states, st_coordinates(st_centroid(states)))

Note the warning, which basically says that centroid coordinates using longitude/latitude data (i.e. WGS84) are not exact, which is perfectly fine for our drawing purposes. State names, which are not capitalized in the data from maps, can be changed to title case using the function toTitleCase from the package tools:

library("tools") states$ID <- toTitleCase(states$ID) head(states) ## Simple feature collection with 6 features and 3 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -124.3834 ymin: 30.24071 xmax: -71.78015 ymax: 42.04937 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## ID X Y geometry ## 1 Alabama -86.83042 32.80316 MULTIPOLYGON (((-87.46201 3... ## 2 Arizona -111.66786 34.30060 MULTIPOLYGON (((-114.6374 3... ## 3 Arkansas -92.44013 34.90418 MULTIPOLYGON (((-94.05103 3... ## 4 California -119.60154 37.26901 MULTIPOLYGON (((-120.006 42... ## 5 Colorado -105.55251 38.99797 MULTIPOLYGON (((-102.0552 4... ## 6 Connecticut -72.72598 41.62566 MULTIPOLYGON (((-73.49902 4...

To continue adding to the map, state data is directly plotted as an additional sf layer using geom_sf. In addition, state names will be added using geom_text, declaring coordinates on the X-axis and Y-axis, as well as the label (from ID), and a relatively big font size.

ggplot(data = world) + geom_sf() + geom_sf(data = states, fill = NA) + geom_text(data = states, aes(X, Y, label = ID), size = 5) + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

We can move the state names slightly to be able to read better “South Carolina” and “Florida”. For this, we create a new variable nudge_y, which is -1 for all states (moved slightly South), 0.5 for Florida (moved slightly North), and -1.5 for South Carolina (moved further South):

states$nudge_y <- -1 states$nudge_y[states$ID == "Florida"] <- 0.5 states$nudge_y[states$ID == "South Carolina"] <- -1.5

To improve readability, we also draw a rectangle behind the state name, using the function geom_label instead of geom_text, and plot the map again.

ggplot(data = world) + geom_sf() + geom_sf(data = states, fill = NA) + geom_label(data = states, aes(X, Y, label = ID), size = 5, fontface = "bold", nudge_y = states$nudge_y) + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

Counties (polygon data)

County data are also available from the package maps, and can be retrieved with the same approach as for state data. This time, only counties from Florida are retained, and we compute their area using st_area from the package sf:

counties <- st_as_sf(map("county", plot = FALSE, fill = TRUE)) counties <- subset(counties, grepl("florida", counties$ID)) counties$area <- as.numeric(st_area(counties)) head(counties) ## Simple feature collection with 6 features and 2 fields ## geometry type: MULTIPOLYGON ## dimension: XY ## bbox: xmin: -85.98951 ymin: 25.94926 xmax: -80.08804 ymax: 30.57303 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## geometry ID area ## 292 MULTIPOLYGON (((-82.66062 2... florida,alachua 2498863359 ## 293 MULTIPOLYGON (((-82.04182 3... florida,baker 1542466064 ## 294 MULTIPOLYGON (((-85.40509 3... florida,bay 1946587533 ## 295 MULTIPOLYGON (((-82.4257 29... florida,bradford 818898090 ## 296 MULTIPOLYGON (((-80.94747 2... florida,brevard 2189682999 ## 297 MULTIPOLYGON (((-80.89018 2... florida,broward 3167386973

County lines can now be added in a very simple way, using a gray outline:

ggplot(data = world) + geom_sf() + geom_sf(data = counties, fill = NA, color = gray(.5)) + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

We can also fill in the county using their area to visually identify the largest counties. For this, we use the “viridis” colorblind-friendly palette, with some transparency:

ggplot(data = world) + geom_sf() + geom_sf(data = counties, aes(fill = area)) + scale_fill_viridis_c(trans = "sqrt", alpha = .4) + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

Cities (point data)

To make a more complete map of Florida, main cities will be added to the map. We first prepare a data frame with the five largest cities in the state of Florida, and their geographic coordinates:

flcities <- data.frame(state = rep("Florida", 5), city = c("Miami", "Tampa", "Orlando", "Jacksonville", "Sarasota"), lat = c(25.7616798, 27.950575, 28.5383355, 30.3321838, 27.3364347), lng = c(-80.1917902, -82.4571776, -81.3792365, -81.655651, -82.5306527))

Instead of looking up coordinates manually, the package googleway provides a function google_geocode, which allows to retrieve geographic coordinates for any address, using the Google Maps API. Unfortunately, this requires a valid Google API key (follow instructions here to get a key, which needs to include “Places” for geocoding). Once you have your API key, you can run the following code to automatically retrieve geographic coordinates of the five cities:

library("googleway") key <- "put_your_google_api_key_here" # real key needed flcities <- data.frame(state = rep("Florida", 5), city = c("Miami", "Tampa", "Orlando", "Jacksonville", "Sarasota")) coords <- apply(flcities, 1, function(x) { google_geocode(address = paste(x["city"], x["state"], sep = ", "), key = key) }) flcities <- cbind(flcities, do.call(rbind, lapply(coords, geocode_coordinates)))

We can now convert the data frame with coordinates to sf format:

(flcities <- st_as_sf(flcities, coords = c("lng", "lat"), remove = FALSE, crs = 4326, agr = "constant")) ## Simple feature collection with 5 features and 4 fields ## Attribute-geometry relationship: 4 constant, 0 aggregate, 0 identity ## geometry type: POINT ## dimension: XY ## bbox: xmin: -82.53065 ymin: 25.76168 xmax: -80.19179 ymax: 30.33218 ## epsg (SRID): 4326 ## proj4string: +proj=longlat +datum=WGS84 +no_defs ## state city lat lng geometry ## 1 Florida Miami 25.76168 -80.19179 POINT (-80.19179 25.76168) ## 2 Florida Tampa 27.95058 -82.45718 POINT (-82.45718 27.95058) ## 3 Florida Orlando 28.53834 -81.37924 POINT (-81.37924 28.53834) ## 4 Florida Jacksonville 30.33218 -81.65565 POINT (-81.65565 30.33218) ## 5 Florida Sarasota 27.33643 -82.53065 POINT (-82.53065 27.33643)

We add both city locations and names on the map:

ggplot(data = world) + geom_sf() + geom_sf(data = counties, fill = NA, color = gray(.5)) + geom_sf(data = flcities) + geom_text(data = flcities, aes(x = lng, y = lat, label = city), size = 3.9, col = "black", fontface = "bold") + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

This is not really satisfactory, as the names overlap on the points, and they are not easy to read on the grey background. The package ggrepel offers a very flexible approach to deal with label placement (with geom_text_repel and geom_label_repel), including automated movement of labels in case of overlap. We use it here to “nudge” the labels away from land into the see, and connect them to the city locations:

library("ggrepel") ggplot(data = world) + geom_sf() + geom_sf(data = counties, fill = NA, color = gray(.5)) + geom_sf(data = flcities) + geom_text_repel(data = flcities, aes(x = lng, y = lat, label = city), fontface = "bold", nudge_x = c(1, -1.5, 2, 2, -1), nudge_y = c(0.25, -0.25, 0.5, 0.5, -0.5)) + coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)

Final map

For the final map, we put everything together, having a general background map based on the world map, with state and county delineations, state labels, main city names and locations, as well as a theme adjusted with titles, subtitles, axis labels, and a scale bar:

library("ggspatial") ggplot(data = world) + geom_sf(fill = "antiquewhite1") + geom_sf(data = counties, aes(fill = area)) + geom_sf(data = states, fill = NA) + geom_sf(data = sites, size = 4, shape = 23, fill = "darkred") + geom_sf(data = flcities) + geom_text_repel(data = flcities, aes(x = lng, y = lat, label = city), fontface = "bold", nudge_x = c(1, -1.5, 2, 2, -1), nudge_y = c(0.25, -0.25, 0.5, 0.5, -0.5)) + geom_label(data = states, aes(X, Y, label = ID), size = 5, fontface = "bold", nudge_y = states$nudge_y) + scale_fill_viridis_c(trans = "sqrt", alpha = .4) + annotation_scale(location = "bl", width_hint = 0.4) + 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(-88, -78), ylim = c(24.5, 33), expand = FALSE) + xlab("Longitude") + ylab("Latitude") + ggtitle("Observation Sites", subtitle = "(2 sites in Palm Beach County, Florida)") + theme(panel.grid.major = element_line(color = gray(0.5), linetype = "dashed", size = 0.5), panel.background = element_rect(fill = "aliceblue"))

This example fully demonstrates that adding layers on ggplot2 is relatively straightforward, as long as the data is properly stored in an sf object. Adding additional layers would simply follow the same logic, with additional calls to geom_sf at the right place in the ggplot2 sequence.

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

Thu, 10/25/2018 - 02:00

view raw Rmd

Beautiful maps in a beautiful world

Maps are used in a variety of fields to express data in an appealing and interpretive way. Data can be expressed into simplified patterns, and this data interpretation is generally lost if the data is only seen through a spread sheet. Maps can add vital context by incorporating many variables into an easy to read and applicable context. Maps are also very important in the information world because they can quickly allow the public to gain better insight so that they can stay informed. It’s critical to have maps be effective, which means creating maps that can be easily understood by a given audience. For instance, maps that need to be understood by children would be very different from maps intended to be shown to geographers.

Knowing what elements are required to enhance your data is key into making effective maps. Basic elements of a map that should be considered are polygon, points, lines, and text. Polygons, on a map, are closed shapes such as country borders. Lines are considered to be linear shapes that are not filled with any aspect, such as highways, streams, or roads. Finally, points are used to specify specific positions, such as city or landmark locations. With that in mind, one need to think about what elements are required in the map to really make an impact, and convey the information for the intended audience. Layout and formatting are the second critical aspect to enhance data visually. The arrangement of these map elements and how they will be drawn can be adjusted to make a maximum impact.

A solution using R and its ecosystem of packages

Current solutions for creating maps usually involves GIS software, such as 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 development environment (IDE) that is used for computing statistical data and graphic in a programmable language, has developed advanced spatial capabilities over the years, and can be used to draw maps programmatically.

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

Using R to create maps brings these benefits to mapping. Elements of a map can be added or removed with ease — R code can be tweaked to make major enhancements with a stroke of a key. It is also easy to reproduce the same maps for different data sets. It is important to be able to script the elements of a map, so that it can be re-used and interpreted by any user. In essence, comparing typical GIS software and R for drawing maps is similar to comparing word processing software (e.g. Microsoft Office or LibreOffice) and a programmatic typesetting system such 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 way to create code that make sense to the user: The grammar of graphics is a term used to breaks up graphs into semantic components, such as geometries and layers. Practically speaking, it allows (and forces!) the user to focus on graph elements at a higher level of abstraction, and how the data must be structured to achieve the expected outcome. While ggplot2 is becoming the de facto standard for R graphs, it does not handle spatial data specifically. The current state-of-the-art of spatial objects in R relies on Spatial classes defined in the package sp, but the new package sf has recently implemented the “simple feature” standard, and is steadily taking over sp. Recently, the package ggplot2 has allowed the use of simple features from the package sf as layers in a graph1. The combination of ggplot2 and sf therefore enables to programmatically create maps, using the grammar of graphics, just as informative or visually appealing as 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 ggplot2 associated to sf, and presents the basics elements and parameters we can 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 of R packages. The full list of packages necessary for this series of tutorials can be installed with:

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-light theme for ggplot2 (theme_bw), which is appropriate for maps:

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 resolution can be set to "low", if preferred):

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 the package sp; we thus convert it to a simple feature using st_as_sf from package sf:

world <- st_as_sf(world) class(world) ## [1] "sf" "data.frame" General concepts illustrated with the world map Data and basic plot (ggplot and geom_sf)

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

ggplot(data = world) + geom_sf()

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

Note that layers are added one at a time in a ggplot call, so the order of each layer is very important. All data will have to be in an sf format to be used by ggplot2; data in other formats (e.g. classes from 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 function ggtitle, passing any valid character string (e.g. with quotation marks) as arguments. Axis names are absent by default on a map, but can be changed to something more suitable (e.g. “Longitude” and “Latitude”), depending on the map:

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 a green color (argument fill), using black for the outline of the countries (argument color):

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

The package ggplot2 allows the use of more complex color schemes, such as a gradient on one variable of the data. Here is another example that shows the population of each country. In this example, we use the “viridis” colorblind-friendly palette for the color gradient (with option = "plasma" for the plasma variant), using the square root of the population (which is stored in the variable POP_EST of the world object):

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, which includes both projection and extent of the map. By default, the map will use 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 the argument crs, it is possible to override this setting, and project on the fly to any projection. This can be achieved using any valid PROJ4 string (here, the European-centric ETRS89 Lambert Azimuthal Equal-Area projection):

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 Petroleum Survey Group (EPSG) code are available for the projection of interest, they can be used directly instead of the full PROJ4 string. The two following calls are equivalent for the ETRS89 Lambert Azimuthal Equal-Area projection, which is EPSG code 3035:

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 practice allowing to “zoom” in the area of interest, provided by limits on the x-axis (xlim), and on the y-axis (ylim). Note that the limits are automatically expanded by a fraction to ensure that data and axes don’t overlap; it can also be turned off to exactly match the limits provided with expand = FALSE:

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 the package ggspatial, which provides easy-to-use functions…

scale_bar that allows to add simultaneously the north symbol and a scale bar into the ggplot map. Five arguments need to be set manually: lon, lat, distance_lon, distance_lat, and distance_legend. The location of the scale bar has to be specified in longitude/latitude in the lon and lat arguments. The shaded distance inside the scale bar is controlled by the distance_lon argument. while its width is determined by distance_lat. Additionally, it is possible to change the font size for the legend of the scale bar (argument legend_size, which defaults to 3). The North arrow behind the “N” north symbol can also be adjusted 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; they can also be set to nautical miles with “nm”, or miles with “mi”.

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 use unprojected data in longitude/latitude (WGS84) on an equidistant cylindrical projection (all meridians being parallel), length in (kilo)meters on the map directly depends mathematically on the degree of latitude. Plots of small regions or projected data will often allow for more accurate scale bars.

Country names and other names (geom_text and annotate)

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

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 map using geographic coordinates. The function requires the data needed to enter 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 many aspects:

  • The size (argument size);
  • The alignment, which is centered by default on the coordinates provided. The text can be adjusted horizontally or vertically using the arguments hjust and vjust, which can either be a number between 0 (right/bottom) and 1 (top/left) or a character (“left”, “middle”, “right”, “bottom”, “center”, “top”). The text can also be offset horizontally or vertically with the argument nudge_x and nudge_y;
  • The font of the text, for instance its color (argument color) or the type of font (fontface);
  • The overlap of labels, using the argument check_overlap, which removes overlapping text. Alternatively, when there is a lot of overlapping labels, the package ggrepel provides a geom_text_repel function that moves label around so that they do not overlap.

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

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)

Final map

Now to make the final touches, the theme of the map can be edited to make it more appealing. We suggested the use of theme_bw for a standard theme, but there are many other themes that can be selected from (see for instance ?ggtheme in ggplot2, or the package ggthemes which provide several useful themes). Moreover, specific theme elements can be tweaked to get to the final outcome:

  • Position of the legend: Although not used in this example, the argument legend.position allows to automatically place the legend at a specific location (e.g. "topright", "bottomleft", etc.);
  • Grid lines (graticules) on the map: by using panel.grid.major and panel.grid.minor, grid lines can be adjusted. Here we set them to a gray color and dashed line type to clearly distinguish them from country borders lines;
  • Map background: the argument panel.background can be used to color the background, which is the ocean essentially, with a light blue;
  • Many more elements of a theme can be adjusted, which would be too long to cover here. We refer the reader to the documentation for the function theme.
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. This function allows a graphic (typically the last plot displayed) to be saved in a variety of formats, including the most common PNG (raster bitmap) and PDF (vector graphics), with control over the size and resolution of the outcome. For instance here, we save a PDF version of the map, which keeps the best quality, and a PNG version of it for web purposes:

ggsave("map.pdf") ggsave("map_web.png", width = 6, height = 6, dpi = "screen")
  1. Note: Support of sf objects is available since version 3.0.0 of ggplot2, recently released on CRAN. 

RSAGA 1.0.0

Sat, 10/20/2018 - 02:00

[view raw Rmd]

RSAGA 1.0.0 has been released on CRAN. The RSAGA package provides an interface between R and the open-source geographic information system SAGA, which offers a variety of geoscientific methods to analyse spatial data. SAGA GIS is supported on Windows, Linux, Mac OS and FreeBSD and is available from sourceforge.

After a long break in the development, RSAGA now supports SAGA GIS 2.3 LTS to 6.2.0. The main issue with the maintenance of RSAGA lies in the changing names of the parameters used by the SAGA command line program. These parameters are essential to provide the user with wrapper functions like rsaga.slope.asp.curv and rsaga.wetness.index. We were able to update all parameter names, while keeping the support for older SAGA versions.

RSAGA is able to find SAGA installations on Windows, Linux and Mac OS automatically again with the build-in function rsaga.env. For new users we would like to refer to the updated vignette.

Further development will go into the integration of Travis CI to test new SAGA versions automatically. We would like to reduce the effort of maintaining the list of parameter names. A solution could be the collaboration with the recently started Rsagacmd project, which uses the html help files of SAGA GIS to get the parameter names.

The development of RSAGA is now integrated into the r-spatial.org community. If you find any bugs please report them at our new development repository r-spatial/RSAGA. We are happy about every project that involves RSAGA. Please help us with bug reports, new ideas and feedback. Stay tuned for further updates.

RSAGA 1.0.0

Mon, 08/20/2018 - 02:00

mapedit and leaflet.js &gt; 1.0

Sun, 07/15/2018 - 02:00

mapedit and leaflet.js &gt; 1.0

Sun, 07/15/2018 - 02:00

RSAGA 1.0.0

Wed, 06/20/2018 - 02:00

Geocomputation with R: workshop at eRum

Thu, 05/31/2018 - 02:00

Geocomputation with R: workshop at eRum

Thu, 05/31/2018 - 02:00

RSAGA 1.0.0

Sun, 05/20/2018 - 02:00

Plotting and subsetting stars objects

Thu, 03/22/2018 - 01:00

Plotting and subsetting stars objects

Thu, 03/22/2018 - 01:00

Pages