Plot maps with base mapping tools and ggmap in R

Plot maps with ‘base’ mapping tools in R

Understanding what kind of data you have (polygons or points?) and what you want to map is pivotal to start your mapping.

  1. First you need a shapefile of the area you want to plot, such as metropolitan France. There are various resources where to get them from: DIVA-GIS and EUROSTAT are those that I use the most. It’s always important to have a .prj file included, as your final map ‘should’ be projecte. I say “should” as sometimes it is just not possible, especially if you work with historical maps.
  2. Upload libraries

Load and prepare data

setwd(paste(mypath))
fr.prj <- readOGR(".", "FRA_adm2")
## OGR data source with driver: ESRI Shapefile
## Source: ".", layer: "FRA_adm2"
## with 96 features
## It has 18 fields
## NOTE: rgdal::checkCRSArgs: no proj_defs.dat in PROJ.4 shared files
map(fr.prj)
rplot
## Warning in SpatialPolygons2map(database, namefield = namefield): database
## does not (uniquely) contain the field 'name'.

head(fr.prj@data)
##   ID_0 ISO NAME_0 ID_1    NAME_1  ID_2         NAME_2   VARNAME_2
## 0   76 FRA France  989    Alsace 13755       Bas-Rhin  Unterelsaá
## 1   76 FRA France  989    Alsace 13756      Haut-Rhin   Oberelsaá
## 2   76 FRA France  990 Aquitaine 13757       Dordogne        <NA>
## 3   76 FRA France  990 Aquitaine 13758        Gironde Bec-D'Ambes
## 4   76 FRA France  990 Aquitaine 13759         Landes      Landas
## 5   76 FRA France  990 Aquitaine 13760 Lot-Et-Garonne        <NA>
##   NL_NAME_2 HASC_2 CC_2      TYPE_2  ENGTYPE_2 VALIDFR_2 VALIDTO_2
## 0      <NA>  FR.BR <NA> Département Department  17900226   Unknown
## 1      <NA>  FR.HR <NA> Département Department  17900226   Unknown
## 2      <NA>  FR.DD <NA> Département Department  17900226   Unknown
## 3      <NA>  FR.GI <NA> Département Department  17900226   Unknown
## 4      <NA>  FR.LD <NA> Département Department  17900226   Unknown
## 5      <NA>  FR.LG <NA> Département Department  17900226   Unknown
##   REMARKS_2 Shape_Leng Shape_Area
## 0      <NA>   4.538735  0.5840273
## 1      <NA>   3.214178  0.4198797
## 2      <NA>   5.012795  1.0389622
## 3      <NA>   9.200047  1.1489822
## 4      <NA>   5.531231  1.0372815
## 5      <NA>   4.489830  0.6062017
# load or create data
set.seed(100)
myvar <- rnorm(1:96)
# manipulate data for the plot
france.geodata  <- data.frame(id=rownames(fr.prj@data), mapvariable=myvar)
head(france.geodata)
##   id mapvariable
## 1  0  1.12200636
## 2  1  0.05912043
## 3  2 -1.05873510
## 4  3 -1.31513865
## 5  4  0.32392954
## 6  5  0.09152878

Use ggmap

# fortify prepares the shape data for ggplot
france.dataframe <- fortify(fr.prj) # convert to data frame for ggplot
## Regions defined for each Polygons
head(france.dataframe)
##       long      lat order  hole piece id group
## 1 7.847912 49.04728     1 FALSE     1  0   0.1
## 2 7.844539 49.04495     2 FALSE     1  0   0.1
## 3 7.852439 49.04510     3 FALSE     1  0   0.1
## 4 7.854333 49.04419     4 FALSE     1  0   0.1
## 5 7.855955 49.04431     5 FALSE     1  0   0.1
## 6 7.856299 49.03776     6 FALSE     1  0   0.1
#now combine the values by id values in both dataframes
france.dat <- join(france.geodata, france.dataframe, by="id")
head(france.dat)
##   id mapvariable     long      lat order  hole piece group
## 1  0    1.122006 7.847912 49.04728     1 FALSE     1   0.1
## 2  0    1.122006 7.844539 49.04495     2 FALSE     1   0.1
## 3  0    1.122006 7.852439 49.04510     3 FALSE     1   0.1
## 4  0    1.122006 7.854333 49.04419     4 FALSE     1   0.1
## 5  0    1.122006 7.855955 49.04431     5 FALSE     1   0.1
## 6  0    1.122006 7.856299 49.03776     6 FALSE     1   0.1
# Plot 3
p <- ggplot(data=france.dat, aes(x=long, y=lat, group=group))
p <- p + geom_polygon(aes(fill=mapvariable)) +
       geom_path(color="white",size=0.1) +
       coord_equal() +
       scale_fill_gradient(low = "#ffffcc", high = "#ff4444") +
       labs(title="Our map",fill="My variable")
# plot the map
p

image-22-02-2017-at-12-11

Use plot basic

nclassint <- 5 #number of colors to be used in the palette
cat <- classIntervals(myvar, nclassint,style = "jenks") #style refers to how the breaks are created
colpal <- brewer.pal(nclassint,"RdBu")
color <- findColours(cat,rev(colpal)) #sequential
bins <- cat$brks
lb <- length(bins)
plot(fr.prj, col=color,border=T)
legend("bottomleft",fill=rev(colpal),legend=paste(round(bins[-length(bins)],1),":",round(bins[-1],1)),cex=1, bg="white")

image-22-02-2017-at-12-23-copy

How to get good maps in R and avoid the expensive softwares

How to convey as much information as possible in a clear and simple way? Producing maps for social sciences is not difficult, there are a plethora of softwares that can help us. But there are a few issues to consider when choosing your to go program:
(1) Do I want to do all my analysis in one (or more) program(s) and then switch to another one to make those maps?
(2) Are those programs freely accessible to me?
1. Using more than one softwares usually implies spending time to learn different syntax:  why do your analysis in (insert name here ____) and then plot in R when you can do everything in R?
2. The availability of mapping softwares is no trivial issue. Not all researchers have powerful computers, not all institutes have bottomless funds to buy licences, and sometimes having the possibility to map on your laptop while bingeing on Netflix is way nicer than waiting for the one computer with the one licence.

Probably the best and most elegant mapping tool available to Geographers is ArcGIS (to my knowledge, but again, I use R and own a Mac), however it does not come for free. What to do? Well, R is a very good alternative, you can produce elegant maps, customizable to the very last detail. The only drawback I have encountered is the time you would spend to get the first map, but then you would have the syntax and any other map would be pretty quick to plot, and you can always for loop all graphics (although I do not recommend it). Moreover, R runs on your Mac (and Linux), it allows for way more control over features, and has great color palettes (see here and here).

Here are some useful libraries:
library(maps) #for creating geographical maps
library(RColorBrewer) #contains color palettes
library(classInt) #defines the class intervals for the color palettes
library(maptools) #tools for handling spatial objects
library(raster) #tools to deal with raster maps
library(ggplot2) #to create maps, quick and painless

Some stuff to keep in mind:
(1) add a scale with scale.map (or a nice  scalebar);
(2) it is sometimes required to add a north arrow, you can find many versions for that (see this document on page 4 for  examples, I use the same with no labels);
(3) locator() is a very useful tool to get the coordinates when adding labels, arrows, scales and so on.

Part 1: get a plain map.

Below is a very simple example produced using EUROSTAT shape files for world countries (world) and DIVA-GIS for Spain at NUTS3 level (spain). In this map I have removed the Canary Islands, but you can always cut it and paste it in the map using either par(fig=c(…)) or par(fin(…)), inset, or something more elaborated with layout, and framing it using box() or rectangle.

world is the shapefile for the whole world, where I select the neighboring countries I want to appear in the map, in this case Spain, France, Portugal, France, Morocco, and Algeria.
spain is the Spain NUTS3 shapefile where I remove the Canary Islands (45)


plot(spain[-c(45),], border=F) #this first line does not plot anything, it just centers my graph on Spain, the -c(45) removes the Canary Islands
plot(world[c(6, 67,74, 132, 177),], border="lightblue",add=T, col="beige") #plotting the countries appearing in the map
plot(spain[-c(45),], border="brown", lwd=0.2, add=T, col="lightblue") #plot spain, removing the Canary Islands
map.scale(3,35.81324, ratio=F, cex=0.7, relwidth=0.1) # scale map
map.axes(cex.axis=0.9)
northarrow(c(4.8,42.9),0.7, cex=0.8)

Spain

Part 2: Add labels

Using the function shadow text to avoid labels overlapping.

coords<- coordinates(spain) # get goordinates of the centroids, it's where you center your labels
# p.names is a data frame containing the coordinates and all the names of the provinces (remember to get rid of those you don't want to use if using only a selection). Usually you can find the names in the shapefile, but I didn't have them.
shadowtext(p.names[,1],p.names[,2], label=paste(p.names[,3]), cex=0.7,col="black", bg="white",r=0.1)

Spain3