R: simple world map (Robinson, ggplot)

Robinson projection seems a bit tricky to work with while using ggplot() in R. At least, at the time of this posting (Dec 2016) I could not find a smooth working solution with coord_map() as discussed here. However, I came across a post of Bob Rudis where he proposes coord_proj() which should do what coord_map does but using proj4::project instead of mapproject. It seems that coord_proj is not implemented in ggplot2. Nevertheless, Bob Rudis explains at the end of his post how to integrate it locally in ggplot2; or, use his package ggalt (though that failed for some updating reasons in ggplot2).

All in all, a general work around is to project everything outside of ggplot with rgdal::spTransform() function. Theoretically, this way should allow for any projection. A bit of care should be given to graticules and their labels as well. Helpful was this link and this SO link. Nevertheless, there are many more interesting projections out there – check this Wikipedia page for some examples and proj4.org for their PROJ.4 strings.

Below is some code example. Of course, things can be done in different ways and below is one perspective. Also some of the code below was helpful for creating this figure in Bennett, Joanne M., et al. “GlobTherm, a global database on thermal tolerances for aquatic and terrestrial organisms.” Scientific data 5 (2018): 180022.

	# ======================================================================================
	# Create a simple world map in Robinson projection with labeled graticules using ggplot
	# ======================================================================================
	# Set a working directory with setwd() or work with an RStudio project
	# __________ Set libraries
	library(rgdal) # for spTransform() & project()
	library(ggplot2) # for ggplot()
	# __________ Load ready to use data from GitHub
	load(url("https://github.com/valentinitnelav/RandomScripts/blob/master/NaturalEarth.RData?raw=true"))
	# This will load 6 objects:
	# xbl.X & lbl.Y are two data.frames that contain labels for graticule lines
	# They can be created with the code at this link:
	# https://gist.github.com/valentinitnelav/8992f09b4c7e206d39d00e813d2bddb1
	# NE_box is a SpatialPolygonsDataFrame object and represents a bounding box for Earth
	# NE_countries is a SpatialPolygonsDataFrame object representing countries
	# NE_graticules is a SpatialLinesDataFrame object that represents 10 dg latitude lines and 20 dg longitude lines
	# (for creating graticules check also the graticule package or gridlines fun. from sp package)
	# (or check this gist: https://gist.github.com/valentinitnelav/a7871128d58097e9d227f7a04e00134f)
	# NE_places – SpatialPointsDataFrame with city and town points
	# NOTE: data downloaded from http://www.naturalearthdata.com/
	# here is a sample script how to download, unzip and read such shapefiles:
	# https://gist.github.com/valentinitnelav/a415f3fbfd90f72ea06b5411fb16df16
	# __________ Project from long-lat (unprojected data) to Robinson projection
	# spTransform() is used for shapefiles and project() in the case of data frames
	# for more PROJ.4 strings check the followings
	# http://proj4.org/projections/index.html
	# https://epsg.io/
	PROJ <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
	# or use the short form "+proj=robin"
	NE_countries_rob <- spTransform(NE_countries, CRSobj = PROJ)
	NE_graticules_rob <- spTransform(NE_graticules, CRSobj = PROJ)
	NE_box_rob <- spTransform(NE_box, CRSobj = PROJ)
	# project long-lat coordinates for graticule label data frames
	# (two extra columns with projected XY are created)
	prj.coord <- project(cbind(lbl.Y$lon, lbl.Y$lat), proj=PROJ)
	lbl.Y.prj <- cbind(prj.coord, lbl.Y)
	names(lbl.Y.prj)[1:2] <- c("X.prj","Y.prj")
	prj.coord <- project(cbind(lbl.X$lon, lbl.X$lat), proj=PROJ)
	lbl.X.prj <- cbind(prj.coord, lbl.X)
	names(lbl.X.prj)[1:2] <- c("X.prj","Y.prj")
	# __________ Plot layers
	ggplot() +
	# add Natural Earth countries projected to Robinson, give black border and fill with gray
	geom_polygon(data=NE_countries_rob, aes(long,lat, group=group), colour="black", fill="gray80", size = 0.25) +
	# Note: "Regions defined for each Polygons" warning has to do with fortify transformation. Might get deprecated in future!
	# alternatively, use use map_data(NE_countries) to transform to data frame and then use project() to change to desired projection.
	# add Natural Earth box projected to Robinson
	geom_polygon(data=NE_box_rob, aes(x=long, y=lat), colour="black", fill="transparent", size = 0.25) +
	# add graticules projected to Robinson
	geom_path(data=NE_graticules_rob, aes(long, lat, group=group), linetype="dotted", color="grey50", size = 0.25) +
	# add graticule labels – latitude and longitude
	geom_text(data = lbl.Y.prj, aes(x = X.prj, y = Y.prj, label = lbl), color="grey50", size=2) +
	geom_text(data = lbl.X.prj, aes(x = X.prj, y = Y.prj, label = lbl), color="grey50", size=2) +
	# the default, ratio = 1 in coord_fixed ensures that one unit on the x-axis is the same length as one unit on the y-axis
	coord_fixed(ratio = 1) +
	# remove the background and default gridlines
	theme_void()
	# save to pdf and png file
	ggsave("map_draft_1.pdf", width=28, height=13.5, units="cm")
	ggsave("map_draft_1.png", width=28, height=13.5, units="cm", dpi=600)
	# REFERENCES:
	# This link was useful for graticule idea
	# http://stackoverflow.com/questions/38532070/how-to-add-lines-of-longitude-and-latitude-on-a-map-using-ggplot2
	# Working with shapefiles, projections and world maps in ggplot
	# http://rpsychologist.com/working-with-shapefiles-projections-and-world-maps-in-ggplot

view raw

World map Robinson projection – ggplot.R

hosted with ❤ by GitHub

Creating the graticule labels data.frames (lbl.Y & lbl.X) used in the code above can be done with:

	# __________ Create data.frame for graticule labels
	# use together with graticule of 10 dg step for latitude and 20 dg step for longitude
	lbl.Y <- data.frame(lon = rep(c(–175,175), each=17), # 175 dg used to fit latitude labels inside Earth; adjust this value to shift positions of labels
	lat = rep(seq(from=80, to=–80, by=–10), times=2))
	lbl.Y$dir <- ifelse(lbl.Y$lat == 0, "", ifelse(lbl.Y$lat > 0, "°N", "°S"))
	lbl.Y$lbl <- paste0(abs(lbl.Y$lat), lbl.Y$dir)
	lbl.X <- data.frame(lon = rep(seq(from=160, to=–160, by=–20), times=2),
	lat = rep(c(–85,85), each=17)) # 85 dg used to fit latitude labels inside Earth; adjust this value to shift positions of labels
	lbl.X$dir <- ifelse(lbl.X$lon == 0, "", ifelse(lbl.X$lon > 0, "°E", "°W"))
	lbl.X$lbl <- paste0(abs(lbl.X$lon), lbl.X$dir)

view raw

Graticule labels 20X10Y.R

hosted with ❤ by GitHub

Downloading data from www.naturalearthdata.com can be done as fallows: (update August 7, 2017: you might want to check the package rnaturalearth)

	# ==========================================================================
	# This is a simple script for downloading, unzipping and reading a shapefile
	# from http://www.naturalearthdata.com
	# ==========================================================================
	library(rgdal)
	# ~~~~~~~~~~~ Download shapefile from http://www.naturalearthdata.com ~~~~~~~~~~~ #
	# __ Example 1 – download countries data
	download.file(url = "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip",
	destfile = "ne_110m_admin_0_countries.zip")
	# unzip the shapefile in the directory mentioned with "exdir" argument
	unzip(zipfile="ne_110m_admin_0_countries.zip", exdir = "ne_110m_admin_0_countries")
	# delete the zip file
	file.remove("ne_110m_admin_0_countries.zip")
	# read the shapefile with readOGR
	NE_countries <- readOGR(dsn = "ne_110m_admin_0_countries", layer = "ne_110m_admin_0_countries")
	# __ Example 2 – download main cities data
	download.file(url = "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_populated_places.zip",
	destfile = "ne_110m_populated_places.zip")
	# unzip the shapefile in the directory mentioned with "exdir" argument
	unzip(zipfile="ne_110m_populated_places.zip", exdir = "ne_110m_populated_places")
	# delete the zip file
	file.remove("ne_110m_populated_places.zip")
	# read the shapefile with readOGR
	NE_places <- readOGR(dsn = "ne_110m_populated_places", layer = "ne_110m_populated_places", stringsAsFactors=FALSE)

view raw

Get Natural Earth data in R.R

hosted with ❤ by GitHub

 

Advertisement

• GIS
• R