This is a tutorial for rendering geographic map using R
Welcome to the training module designed for postdocs and alumni fellows in the Climap Africa programme. We will work today with the open data science tool R.
Note:
The powerpoint file, script, and datafiles are freely available. Just click on the corresponding hyperlink.
Course objective: Enhance R users for informative geographic maps rendering using R
Notice: Please visit the hyperlink (in blue) for installation and settings.
Pre-requisites:
-
Install R and RStudio on Windows 7, 8 or 10. A tutorial for a beginner is here.
-
Install the following packages before the course:
A tutorial for package installation in RStudio is here.
It is more easier to install CRAN packages at once by typing
install.packages(c("rgdal", "mapdata", "mapproj" ,"maps" ,"ggplot2", "ggrepel", "dplyr", "scales", "ggmap", "gpclib", "rgeos"))For installation of the legendMap package, you need to install the remotes package first.
install.packages("remotes")Then install legendMap
remotes::install_github("3wen/legendMap", force = TRUE)Please skip this part if you do not encounter package installation troubles. Continue here
Based on our recent experience on the online session, some participants have some troubles regarding this warning message:
Error in maptools::unionSpatialPolygons(cp, attr[, region]) :
isTRUE(gpclibPermitStatus()) is not TRUEDon't worry at all!
It is easy to solve this issue.
Just run the following codes:
- You need to install gpclib package by doing
install.packages("gpclib", dependencies=TRUE)- In addition you need to install an other package call rgeos
install.packages("rgeos", dependencies=TRUE)- After you've run the previous codes, please check if the libraries are well installed:
library(gpclib)
library(rgeos)- Then run your R code for a nice map.
- Download the data for exercise here
| # | Section | Duration |
|---|---|---|
| 1. | Working directory setting, data preparation and shapefile importation | 5 min |
| 2. | Rendering a basic map in R using ggplot2 | 10 min |
| 3. | Rendering a choropleth map | 20 min |
| 4. | Add a scale bar and North Arrow | 5 min |
| 5. | Tips | 5 mn |
| 6. | Q & A | 15 min |
| Total | 60 min |
The working directory is the folder named R MAP. Please put all the shapefiles and data in your working disrectory. To set your working directory, type:
setwd("C:/Users/ANGE/Documents/R MAP")rm(list=ls())mySHP="C:/Users/ANGE/Documents/R MAP"myFile=readOGR (mySHP, layer="SEN_adm1", stringsAsFactors=FALSE)This result should appear
OGR data source with driver: ESRI Shapefile
Source: "C:\Users\ANGE\Documents\R MAP", layer: "SEN_adm1"
with 14 features
It has 9 fields
Integer64 fields read as strings: ID_0 ID_1 class(myFile)The class should display as follows
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"names(myFile)The name will display as follows:
[1] "ID_0" "ISO" "NAME_0" "ID_1" "NAME_1" "TYPE_1"
[7] "ENGTYPE_1" "NL_NAME_1" "VARNAME_1"
print(myFile$NAME_1)
The regions names should appear as follows:
[1] "Dakar" "Diourbel" "Fatick" "Kédougou" "Kaffrine"
[6] "Kaolack" "Kolda" "Louga" "Matam" "Sédhiou"
[11] "Saint-Louis" "Tambacounda" "Thiès" "Ziguinchor"
You can import the packages like this:
library(rgdal)
library(mapdata)
library(mapproj)
library(maps)
library(ggplot2)
library(ggrepel)
library(legendMap)
library(dplyr)
library(scales)
library(ggmap)
library(gpclib)
library(rgeos)
Or, to load multiple packages at once, type:
Packages = "rgdal", "mapdata", "mapproj" ,"maps" ,"ggplot2", "ggrepel", "legendMap", "dplyr", "scales", "ggmap", "gpclib","rgeos")
lapply(Packages, library, character.only = TRUE)
myDF = fortify(myFile, region = "NAME_1")
Type:
head(myDF, 4)
| long | lat | order | hole | piece | id | group |
|---|---|---|---|---|---|---|
| -17.16056 | 14.89375 | 1 | FALSE | 1 | Dakar | Dakar.1 |
| -17.16004 | 14.89333 | 2 | FALSE | 1 | Dakar | Dakar.1 |
| -17.16000 | 14.89335 | 3 | FALSE | 1 | Dakar | Dakar.1 |
| -17.15683 | 14.89042 | 4 | FALSE | 1 | Dakar | Dakar.1 |
myDF = rename(myDF, Longitude = long, Latitude= lat)
Type:
head(myDF, 4)
| Longitude | Latitude | order | hole | piece | id | group |
|---|---|---|---|---|---|---|
| -17.16056 | 14.89375 | 1 | FALSE | 1 | Dakar | Dakar.1 |
| -17.16004 | 14.89333 | 2 | FALSE | 1 | Dakar | Dakar.1 |
| -17.16000 | 14.89335 | 3 | FALSE | 1 | Dakar | Dakar.1 |
| -17.15683 | 14.89042 | 4 | FALSE | 1 | Dakar | Dakar.1 |
p <- ggplot() +
geom_polygon(data = myDF,
aes(x = Longitude, y = Latitude, group = group),
color = "black", size = 0.25) +
coord_map() +
theme_minimal() +
ggtitle("Basic map with ggplot2")The map is in p.
p
You should get this output:
3.1. Import the data we want to plot on the map.Here that is the production of pearl millet per region
mydata = read.csv("production_data.csv", header=TRUE, sep=";")
mydata1 = read.csv("region_names.csv", header=TRUE, sep=";")
Type:
head(mydata, 4)
| long | lat | id | Production |
|---|---|---|---|
| -17.33 | 14.75 | Dakar | 0 |
| -16.25 | 14.75 | Diourbel | 46231 |
| -16.53 | 14.36 | Fatick | 80000 |
| -12.18 | 12.80 | Kédougou | 152 |
Type:
head(mydata1, 4)
| Region | long | lat |
|---|---|---|
| DAKAR | -17.33 | 14.75 |
| DIOURBEL | -16.25 | 14.75 |
| FATICK | -16.53 | 14.36 |
| KAOLACK | -16.00 | 14.00 |
plotData <- left_join(myDF, mydata)Key point: Note that myDF and mydata has id as a common variable
Type:
head(plotData)You will get this table
| Longitude | Latitude | order | hole | piece | id | group | long | lat | Production |
|---|---|---|---|---|---|---|---|---|---|
| -17.16056 | 14.89375 | 1 | FALSE | 1 | Dakar | Dakar.1 | -17.33 | 14.75 | 0 |
| -17.16004 | 14.89333 | 2 | FALSE | 1 | Dakar | Dakar.1 | -17.33 | 14.75 | 0 |
| -17.16000 | 14.89335 | 3 | FALSE | 1 | Dakar | Dakar.1 | -17.33 | 14.75 | 0 |
| -17.15683 | 14.89042 | 4 | FALSE | 1 | Dakar | Dakar.1 | -17.33 | 14.75 | 0 |
p <- ggplot() +
# English: Plot the shapefile geographic information
# French: Projeter les contours geographiques du shapefile
geom_polygon(data = plotData,
aes(x = Longitude,
y = Latitude,
group = group,
fill = Production),
color = "black", size = 0.25) +
#-------------------------------------------------------------------
#English: Map projection type. Default is "mercator".
#See help(coord_map() for other choice)
coord_map()+
#------------------------------------------------------------------
# English: Color choice, direction = 1 means colored from lowest to highest value
# French: Choix de couleur. direction=1 veut dire colorer par ordre croissant de valeur
scale_fill_distiller(palette = "Greens",direction=1) +
#------------------------------------------------------------------
# English: Plot the localization points corresponding to each region
# French: Projeter les points geographiques de chaque region
geom_point(data=mydata1,
aes(x=long, y=lat),
shape = 21,
fill = "white",
size = 3,
color = "black") +
#---------------------------------------------------------------------
# English: Avoid overlapping text for the points annotation
# French: Eviter la superposition des noms de regions
geom_label_repel(data=mydata1,
aes(x=long,
y=lat,
label=Region),
fontface = 'bold',
color = 'black',
box.padding = 0.35,
point.padding = 0.5,
segment.color = 'grey10')+
# English: Theme option
# French: Mise forme de l'arriere plan
theme_minimal() + # Simple background
theme(panel.grid.major =
element_line(
colour = "black",
size = 0.5,
linetype = "dotted")) + # Customize the grid line type, size and color
theme(plot.background =
element_rect(
colour = "white",
size = 1)) + # Customize the background line type, color and size
# Add a title
ggtitle("Map of Pearl Millet Production in Senegal (Rainy season 2017)")
Just call the map variable p
p
That is it!
p <- ggplot() +
geom_polygon(data = plotData, aes(x = Longitude, y = Latitude, group = group,
fill = Production), color = "black", size = 0.25) +
coord_map()+
scale_fill_distiller(palette = "Greens",direction=1) +
geom_point(
data=mydata1,
aes(x=long, y=lat),
shape = 21,
fill = "white",
size = 3,
color = "black") +
geom_label_repel(data=mydata1, aes(x=long, y=lat, label=Region),
fontface = 'bold', color = 'black',
box.padding = 0.35, point.padding = 0.5,
segment.color = 'grey10')+
# English: Scale bar and north arrow
# French: Ajouter le north geographique et l'echelle
scale_bar(lon = -12, lat = 16, #--longitude and latitude of the scale bar position
distance_lon = 40, distance_lat = 10, #--legnth and width of each rectangle
distance_legend = 25, #--distance between legend rectangles and legend texts
dist_unit = "km", #-- Unit
arrow_length = 10, #-- Arrow length
arrow_distance = 50, #-- Arrow distance to the scale bar
arrow_north_size = 6) + #-- Arrow size
theme_minimal() +
theme(panel.grid.major =
element_line(
colour = "black",
size = 0.5,
linetype = "dotted")) +
theme(plot.background =
element_rect(
colour = "white",
size = 1)) +
ggtitle("Map of Pearl Millet Production in Senegal (Rainy season 2017)")
Just call the map variable p
p
Export a high quality map by typing:
PDF format
ggsave(p, file = "carte.pdf", limitsize = FALSE, width = 12, height = 10.5, dpi=500 )
PNG format
ggsave(p, file = "carte.png", limitsize = FALSE, width = 10, height = 6.5, type = "cairo-png", dpi=500)
To find out a desirable position for scale bar or any adjustment, it is possible to plot the map with the basic R by typing:
plot(myFile, axes=T, col="aliceblue")You will get this output
and then type:
locator(n=2)
2 is just an example. You can define many number as much as possible. Using your mouse, click on the position you want. You will get the coordinates.
- Pebesma EJ, Bivand RS (2005). Classes and methods for spatial data in R. R News, 5(2), 9–13. https://CRAN.R-project.org/doc/Rnews/.
- Bivand RS, Pebesma E, Gomez-Rubio V (2013). Applied spatial data analysis with R, Second edition. Springer, NY. https://asdar-book.org/.
- Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.
Scripts in this repository are licensed under General Public License v3.
For information inquiry please email Yedomon Ange Bovys Zoclanclounon PhD Candidate (angez9914@gmail.com)