Map with ggplot2: Tracking the mildew prevalence in pearl millet production zones in Senegal

Data source

The data used for this tutorial are available here

Note

The senegal shape file was downloaded here

This tutorial is based on this published paper

Make sure you've installed all requiered packages.

Code

#Remove all objects from our current workspace

rm(list=ls())

#Libraries loading

library(rgdal)
## Loading required package: sp

## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/ANGE/Documents/R/win-library/3.5/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/ANGE/Documents/R/win-library/3.5/rgdal/proj
##  Linking to sp version: 1.3-1
library(mapdata)
## Loading required package: maps
library(mapproj)
library(maps)
library(ggplot2)
library(ggrepel)
library(legendMap)
## Loading required package: maptools

## Checking rgeos availability: TRUE

## Loading required package: grid
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#Set the working directory

setwd("C:/Users/ANGE/Documents/R MAP")

#Set the shapefile path

mySHP="C:/Users/ANGE/Documents/R MAP"

#Import the districts shapefile

myFile=readOGR (mySHP, layer="SEN_adm2", stringsAsFactors=FALSE)
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\ANGE\Documents\R MAP", layer: "SEN_adm2"
## with 45 features
## It has 11 fields
## Integer64 fields read as strings:  ID_0 ID_1 ID_2
#Focus on the districts of interest

myregion = subset(myFile,
myFile$ID_2 == "5" |
myFile$ID_2 == "6" |
myFile$ID_2 == "7" |
myFile$ID_2 == "8" |
myFile$ID_2 == "9" |
myFile$ID_2 == "10" |
myFile$ID_2 == "14" | 
myFile$ID_2 == "15" | 
myFile$ID_2 == "16" |
myFile$ID_2 == "17" |
myFile$ID_2 == "18" |
myFile$ID_2 == "19" |
myFile$ID_2 == "20" |
myFile$ID_2 == "21" |
myFile$ID_2 == "22" |
myFile$ID_2 == "23" |
myFile$ID_2 == "30" |
myFile$ID_2 == "31" |
myFile$ID_2 == "32" |
myFile$ID_2 == "36" |
myFile$ID_2 == "37" |
myFile$ID_2 == "38" |
myFile$ID_2 == "39" |
myFile$ID_2 == "40" |
myFile$ID_2 == "41" |
myFile$ID_2 == "42" )

#Transform into a dataframe

myreg = fortify(myregion)
## Regions defined for each Polygons
myDF = fortify(myFile)
## Regions defined for each Polygons
##Rename long et lat
myreg = rename(myreg, Longitude = long, Latitude= lat)

#Mildew incidence data importation
mydata = read.csv("dat_revu.csv", header=TRUE, sep=";")


#Make the plot

mybreaks = c(1.5, 10, 25, 50, 75, 100)
p <-  ggplot() +
      geom_polygon(data=myreg, linetype = "dotted", aes(x=Longitude, y=Latitude, group = group),colour="black", fill="white")+
      geom_point(data = mydata, aes(x=long, y=lat, size= Inc,color=Region)) +
      scale_size_continuous(name = "Downy\nMildew\nIncidence\n(%)", breaks = mybreaks) +
      coord_map(projection = "lagrange") +
      scale_bar(lon = -17, lat = 13, 
                               distance_lon = 40, distance_lat = 10, 
                               distance_legend = 25, dist_unit = "km", 
                               arrow_length = 10, arrow_distance = 50, arrow_north_size = 6)+
      theme_minimal() +
      theme(panel.grid.major = element_line(colour = "black", size = 0.5, linetype = "dotted")) +
      theme(plot.background = element_rect(colour = "white", size = 1)) 

#You can save using the following code: ggsave(pb, file = "mapinci.png", width = 10, height = 6.5, type = "cairo-png")

Map rendering

Just call the map variable p

Infected Areas Map following the productive zones in Senegal

mybreaks = c(100, 1000, 5000, 10000, 20000, 30000, 40000)
g <-  ggplot()+
      geom_polygon(data=myreg, linetype = "dotted", aes(x=Longitude, y=Latitude, group = group),colour="black", fill="white")+
      geom_point(data = mydata, aes(x=long, y=lat, size=Arae,color=Region))+
      scale_size_continuous(name = "Downy\nMildew\nInfected\nArea\n(meter square)", breaks = mybreaks)+
      coord_map(projection = "lagrange")+
      scale_bar(lon = -17, lat = 13, 
                               distance_lon = 40, distance_lat = 10, 
                               distance_legend = 25, dist_unit = "km", 
                               arrow_length = 10, arrow_distance = 50, arrow_north_size = 6)+
      theme_minimal()+
      theme(panel.grid.major = element_line(colour = "black", size = 0.5, linetype = "dotted"))+
      theme(plot.background = element_rect(colour = "white", size = 1))

Map rendering

Just call the map variable g

That is it!