Maps, Coordinate Reference Systems and Visualising Geographic Data with mapmisc

The mapmisc package provides functions for visualising geospatial data, including fetching background map layers, producing colour scales and legends, and adding scale bars and orientation arrows to plots. Background maps are returned in the coordinate reference system of the dataset supplied, and inset maps and direction arrows reflect the map projection being plotted. This is a “light weight” package having an emphasis on simplicity and ease of use.

Patrick Brown (Cancer Care Ontario)
2016-07-28

1 Introduction

R has extensive facilities for managing, manipulating, and visualising spatial data, with the sp (Pebesma and Bivand 2005) and raster (Hijmans 2015b) packages providing a set of object classes and core functions which numerous other packages have built on. It is fairly straightforward to import spatial data of a variety of types and from a range of sources including: images for map backgrounds; high-resolution pixel grids of surface elevation; and polygons of administrative region boundaries. Large volumes of such data are available for download from sites such as worldgrids.org, gadm.org, and nhgis.org, and map images are freely available from OpenStreetMap.org and other online maps. The first issue often encountered after downloading and importing spatial data is reconciling different coordinate reference systems (CRS’s, or map projections). Most repositories of spatial data provide longitude-latitude coordinates, although single-country data sources often use a country-specific map projection (i.e. the UK’s Ordinance Survey National Grid) and online maps mostly use the Web Mercator projection. The suitability of a particular map projection will depend on the geographic region being considered and the specific problem at hand.

The mapmisc package (Brown 2016) provides tools for working with projected data which cover the following four areas:

This paper will cover each of these points in turn, working through examples and briefly describing the operations by the functions in the mapmisc package. An emphasis is given to tidy, intuitive, and reproducible code accessible for students and non-specialists.

The two most important packages required for using spatial data in R are the sp and raster packages, which provide tools and classes for vector data (spatial data on a continuous domain) and raster data (defined on a pixelated grid) respectively. Installing mapmisc with install.packages("mapmisc") or by using a menu item on a GUI will install sp and raster if they are not already present. A third important spatial package is rgdal (Bivand et al. 2016), which provides methods for re-projecting coordinates and importing spatial data in various file formats. The Geographic Data Abstraction Language (GDAL) underlies rgdal, aligning with R’s UNIX-like philosophy of combining separate and specialised pieces of software. On most UNIX-based systems, the GDAL and proj4 software must be installed separately prior to installing rgdal. All versions of Windows and most versions of MacOS have binary versions of rgdal which include the GDAL and proj4 binaries, and rgdal can be installed in the same manner as any other R package.

Additional packages used by mapmisc are: RColorBrewer (Neuwirth 2014), classInt (Bivand 2015), rgeos (Bivand and Rundel 2016), and geosphere (Hijmans 2015a). These four packages and rgdal are not always installed automatically with mapmisc, as they are marked as “suggested” packages with mapmisc being usable with a reduced level of functionality without them. Three further packages necessary for reproducing the examples in this paper are dismo (Hijmans et al. 2016), maptools (Bivand and Lewin-Koh 2016), and R.utils (Bengtsson 2016).

Loading the packages with

library("rgdal")
library("mapmisc")

also makes sp and raster available. The remaining packages do not need to be loaded explicitly and will be called by mapmisc as needed.

Getting started with spatial data in R

The getData function provided by raster is able to download a number of useful and interesting spatial datasets. The coastline and borders of Finland can be fetched with

finland <- raster::getData("GADM", country = "FIN", level = 0)

The object finland is a “SpatialPolygonsDataFrame”, and Bivand et al. (2013) contains a wealth of information on working with objects of this type. The command

plot(finland, axes = TRUE)

produces the plot in Figure 1a.

The choice of Finland as an example is due to its being far from the equator, and a useful contrast on the other side of the world is New Zealand. The coastline of New Zealand, obtained with

nz <- raster::getData("GADM", country = "NZL", level = 0)

includes a number of small outlying islands. The spatial extent of the nz object,

raster::extent(nz)
## class       : Extent 
## xmin        : -179 
## xmax        : 179 
## ymin        : -52.6 
## ymax        : -29.2

spans the entire globe in the \(x\)-direction since New Zealand has islands on both sides of the \(180^{\circ}\) meridian. Finding an appropriate axis limit through trial and error brings one to

plot(nz, xlim = c(167, 178), axes = TRUE)

and the map in Figure 1b.

The outlying islands can be removed from the nz object using the crop function in the raster package, which in turn calls rgeos. Using the locator function and a few iterations of trial and error leads to the discovery that a region spanning 160 to 180 degrees longitude and \(-47\) to \(-30\) degrees latitude boxes in the main islands of New Zealand fairly tightly. The parts of New Zealand contained within this box can be extracted by creating an extent object and passing it to crop.

nzClip <- raster::crop(nz, extent(160, 180, -47, -30))

The finland and nzClip objects will be used in the mapping examples which follow.

2 Working with map projections

This section covers mapping projected data and defining customised map projections. Adding background images, scale bars, and inset maps to plots with the map.new, openmap, scaleBar, and insetMap functions is demonstrated in the production of Figure 2. Map projections suitable for displaying the entire globe are constructed with the moll function, and along with the wrapPolys function Figure 5 is made. Map projections where Euclidean distances from \(x\)-\(y\) coordinates are useful approximations of shortest distances between points on the globe are obtained with the omerc function and used to produce Figure 7.

Spatial data with coordinate reference systems

The spatial coordinates in Figures 1a and 1b are angles of longitudes and latitudes; coordinates which would be equivalent to the two angles of the spherical coordinate system \((\rho, \theta, \phi)\) familiar to mathematicians were it not for the inconvenient fact that the Earth is not spherical. The Earth is rather an oblate spheroid, slightly “squashed” or pumpkin-shaped, and the angles of orientations of lines pointing directly “up” (with reference to the stars) and “down” (as defined by a plumb line pulled straight by gravity) differ. As a result, various types of long-lat coordinates are in use, with the World Geodesic System (WGS84) used by Global Positioning Systems being the most widespread. The European Petroleum Standards Group (EPSG) catalogue of Coordinate Reference Systems (or CRS’s) refers to this system by the code 4326, and this code can be used to create an R object of class “CRS” corresponding to the WGS84 system using the CRS function from the sp package.

sp::CRS("+init=epsg:4326")
## CRS arguments:
##  +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## +towgs84=0,0,0

The syntax of +argument=value for specifying a CRS comes from the PROJ.4 Cartographic Projections Library, with +proj=longlat indicating coordinates are angles and ellps=WGS84 specifying that the Earth is an ellipsoid with values of the major axis length, minor axis length, and flattening corresponding to the WGS84 specification.

The difficulty angular coordinates pose for interpreting maps or using spatial statistics is that Euclidean distance \(\sqrt{x^2 + y^2}\) is not always a useful measure of the distance separating two points. The shortest route between two points on a sphere follows a Great Circle which divides the globe into two equal halves. The distance between two points along this path, the Great Circle distance (see Wikipedia 2015a), can be computed with a trigonometric formula as implemented in the spDists function in sp. Euclidean distance will be roughly proportional to Great Circle distance for two points near the equator and reasonably close to one another. In Finland and New Zealand, however, Euclidean distance will over-emphasise the east-west direction since one degree of longitude is a much shorter distance (in kilometres) than one degree of latitude. It is for this reason that Greenland appears larger than India on many maps even though the opposite is true. Most R packages which perform spatial analyses rely on Euclidean distance, including this author’s geostatsp (Brown 2015), even though Great Circle distance would be straightforward to implement. Fitting a spatial model with geostatsp to data in long-lat coordinates from Finland might uncover directional effects with strong correlation in the east-west direction, which could well be an artefact arising from the over-estimation of east-west distances. The importance of transforming spherical coordinates to a coordinate system where the Euclidean distance is a reasonable approximation to the Great Circle distance should not be under-estimated.

Most countries have an “official” CRS which produces accurate Euclidean distances for specific areas of the globe, one of which is the Finland Uniform Coordinate System having EPSG code 2393. This projection is obtained in R by CRS with

CRS("+init=epsg:2393")
## CRS arguments:
##  +init=epsg:2393 +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +y_0=0
## +ellps=intl +towgs84=-96.062,-82.428,-121.753,4.801,0.345,-1.376,1.496
## +units=m +no_defs

This is a Transverse Mercator projection (+proj=tmerc) with \(x\) and \(y\) coordinates giving positions on a cylinder containing the earth. The entry +lon_0=27 indicates that the cylinder touches the Earth along the \(27^\circ\) meridian line. Provided two points are reasonably close to the \(27^\circ\) meridian, Euclidean distance between their Finland Uniform Coordinates will be very close to the true distance between them. The map of Finland can be converted to this coordinate system using spTransform from sp and rgdal with

finlandMerc <- spTransform(finland, CRS("+init=epsg:2393"))

Figure 1c is the result of plotting this object with plot(finlandMerc, axes = TRUE). Notice the projected map has a wider base and narrower top than the long-lat map in Figure 1a. The coordinates in Figure 1c refer to an origin where the \(27^\circ\) meridian intersects the equator, with the +x_0= argument above indicating that 3,500km are added to the \(x\) coordinates.

graphic without alt text graphic without alt text graphic without alt text graphic without alt text
  1. Finland, Longitude-Latitude
  1. New Zealand, Longitude-Latitude
  1. Finland Transverse Mercator
  1. New Zealand, Oblique Mercator
Figure 1: Basic maps of Finland and New Zealand in different Coordinate Reference Systems.

Cylindrical map projections can be constructed from any cylinder containing the Earth, and there is no mathematical requirement to use one of the standard transverse Mercator projections with an EPSG number. For example, a given user might consider the point \((170^\circ\text{E}, 45^\circ\text{S})\) to be an intuitive location for the origin of the transformed map of New Zealand, and thus decide to define a custom CRS with this centre. The cylinder can follow any Great Circle, it need not be a meridian line, and a Great Circle angled \(40^\circ\) clockwise would run the length of the two islands. Cylindrical projections with an angle are termed Oblique Mercator projections (see Snyder 1987 66), and can be constructed with the assistance of mapmisc’s omerc function. A custom projection with the \((170^\circ\text{E}, 45^\circ\text{S})\) origin and a \(40^\circ\) rotation is obtained by

nzCrs <- omerc(c(170, -45), angle = 40)
nzCrs
## CRS arguments:
##  +proj=omerc +lat_0=-45 +lonc=170 +alpha=40 +k=1 +x_0=0 +y_0=0 +gamma=0
## +ellps=WGS84 +units=m

The difference between the above and the projection for Finland is the omerc in place of tmerc, with the additional argument +alpha=40 specifying an angle. The New Zealand coastline can be projected to this CRS with spTransform.

nzRot <- spTransform(nzClip, nzCrs)

Figure 1d results from executing plot(nzRot, axes = TRUE), and New Zealand has been rotated \(40^\circ\) to a vertical position.

Finding a projection

When choosing a map projection for a dataset, a simple web search of a phrase such as “map projection finland epsg” will often give clear advice as to what the most commonly used national CRS is. A number of tools in rgdal can be used to obtain a projection in a more systematic manner. Below the make_EPSG function creates a table of all EPSG coded CRS’s which rgdal supports, and a grep command used to show all those projections with “Finland” in its description. The resulting list of projections below confirms that 2393 is a sensible choice.

allEpsg <- rgdal::make_EPSG()
allEpsg[grep("Finland", allEpsg$note), 1:2]
##      code                                                   note
## 859  2391                                 # KKJ / Finland zone 1
## 860  2392                                 # KKJ / Finland zone 2
## 861  2393              # KKJ / Finland Uniform Coordinate System
## 862  2394                                 # KKJ / Finland zone 4
## 1853 3386                                 # KKJ / Finland zone 0
## 1854 3387                                 # KKJ / Finland zone 5
## 4929 3901 # KKJ / Finland Uniform Coordinate System + N60 height

A second method for obtaining a CRS is to make a rough guess at a projection string and use showEPSG to attempt to find a corresponding EPSG code. Were one to use a Universal Transverse Mercator (or UTM) projection for a map of Finland, a web search for “UTM zone map” shows that Finland lies in UTM zone 35. A proj4 specification of a UTM zone 35 projection will contain +proj=utm and +zone=35, and showEPSG states that the EPSG code 32635 corresponds to an appropriate CRS.

rgdal::showEPSG("+proj=utm +zone=35")
## [1] "32635"
CRS("+init=epsg:32635")
## CRS arguments:
##  +init=epsg:32635 +proj=utm +zone=35 +datum=WGS84 +units=m +no_defs
## +ellps=WGS84 +towgs84=0,0,0

The definitive resource for information on national map projections is contained in the monthly bulletins of The Imaging and Geospatial Information Society. The GridsDatums data set in rgdal gives the year and month for each country’s entry at www.asprs.org/Grids-Datums.html. The entry for Finland appears in the October 2006 issue.

data("GridsDatums", package = "rgdal")
GridsDatums[grep("Finland", GridsDatums$country), ]
##                 country     month year ISO
## 100 Republic of Finland (October) 2006 FIN

Mapping projected data

The mapmisc functions openmap, map.new, scaleBar, and insetMap can be used together to improve on the basic maps in Figure 1, and they are used here to add background images, a scale bar, and an inset map to Figure 2.

Background map images are obtained from the openmap function, which downloads image files from OpenStreetMap.org or a number of other sources. The images used in Figure 2 were obtained with

nzBg <- openmap(nzRot, path = "mapquest-sat")
finlandBg <- openmap(finlandMerc, path = "landscape")
graphic without alt text graphic without alt text
  1. Finland. Background © OpenStreetMap
  1. New Zealand. Tiles courtesy of MapQuest and © OpenStreetMap
Figure 2: Maps produced using the mapmisc package, containing background images, inset maps, and scale bars.

The first argument of openmap is used to set both the spatial extent of the map to be retrieved and the CRS the map will be projected to. Any spatial object x for which extent(x) and projection(x) are defined can be provided to openmap. The path = argument specifies the source of the map, and sample maps from the various sources are shown in Figure 11 in the Appendix. The objects produced by openmap are raster objects, converted from image files downloaded from web map servers. The Finland map is an object of class “RasterLayer”.

finlandBg
## class       : RasterLayer 
## dimensions  : 768, 376, 288768  (nrow, ncol, ncell)
## resolution  : 2747, 2252  (x, y)
## extent      : 2889952, 3922969, 6183827, 7913059  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:2393 +proj=tmerc +lat_0=0 +lon_0=27 +k=1 +x_0=3500000 +...
## data source : in memory
## names       : landscape 
## values      : 1, 1022  (min, max)

Notice that the CRS is the same as for the finlandMerc object. The New Zealand map is a “RasterStack” with red, green and blue layers.

nzBg
## class       : RasterStack 
## dimensions  : 512, 289, 147968, 3  (nrow, ncol, ncell, nlayers)
## resolution  : 5237, 5190  (x, y)
## extent      : -1e+06, 510285, -977333, 1680003  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=omerc +lat_0=-45 +lonc=170 +alpha=40 +k=1 +x_0=0 +y_0=0 +gam...
## names       : mapquest.satRed, mapquest.satGreen, mapquest.satBlue 
## min values  :               0,                 0,                0 
## max values  :             248,               250,              247

The Finland map can be viewed with plot(finlandBg), whereas the three-layered New Zealand map needs the plotRGB function for plotting its red, green and blue values as colours.

Figure 2a is produced with the four function calls below.

map.new(finlandMerc, 0.8)
plot(finlandBg, add = TRUE)
plot(finlandMerc, add = TRUE)
scaleBar(finlandMerc, "bottomright")
insetMap(finlandMerc, "right", width = 0.3, cropInset = extent(0, 180, -50, 70))

The functions run are the following:

The New Zealand map in Figure 2b is produced with similar code.

map.new(nzRot, 0.8)
plotRGB(nzBg, add = TRUE)
rgdal::llgridlines(nzRot, col = "yellow")
plot(nzRot, add = TRUE, border = "red")
scaleBar(nzRot, "left")
insetMap(nzRot, "right", width = 0.3, cropInset = extent(0, 180, -50, 70))

The use of plotRGB in place of plot is used for the background map, and the scale bar has been placed at the centre-left. The llgridlines function from rgdal added latitude and meridian lines in yellow.

The images in Figure 2 have dimensions of 4 by 5 inches, saved as png files with 72 pixels per inch. Executing the code above in an interactive R session will likely produce maps with a slightly different appearance unless the graphics window has these same dimensions. This document is produced with knitr (Xie 2015), and the figure dimensions are set with the fig.height = and fig.width = options to code chunks.

The map images in finlandBg and nzBg were retrieved from OpenStreetMap.org and MapQuest respectively, and although they are free to use and reproduce they must be attributed. The openmapAttribution function produces an attribution for an object produced by openmap or a string valid as a path = argument for openmap. An attribution for nzBg (or "mapquest-sat"), as in the caption for Figure 2, is produced with

openmapAttribution(nzBg, short = TRUE)
##                                   mapquest.sat 
## "Tiles courtesy of MapQuest(www.mapquest.com)"

The Acknowledgements section of this paper has used this function without the short = TRUE argument.

openmapAttribution(finlandBg)
##   landscape 
## "copyright OpenStreetMap.org contributors. Data by OpenStreetMap.org
 #  available under the Open Database License (opendatacommons.org/licenses/odbl),
 #  cartography by Thunderforest.com"

The additional argument type = "latex" is used in the source code for this paper, and type = "markdown" is also available.

Projecting background maps

There are a number of potential pitfalls involved when using background map images with projected data, and this section will describe some additional options to openmap and map.new which can help in this regard.

Two undesirable features of Figure 2a are the white triangular section in the top left of the map, and the low resolution and lack of legibility of the names of towns and cities. How this arose can be understood by contrasting the map image retrieved from OpenStreetMap.org in Figure 3a with the Finland Uniform Coordinate System map in Figure 3b. The map in Figure 3a uses coordinates in the Spherical Mercator projection, where a vertically-oriented cylinder is wrapped around a spherical Earth at the equator,1 and the rectangular area covered by the original map becomes somewhat trapezoidal when projected to the Transverse Mercator coordinates in Figure 3b. The transformation has distorted the text on the image, and the image does not completely cover the black rectangle corresponding to the plotting region of Figure 2a.