Introduction

Welcome to Module V of our Spatial Statistics and Disease Mapping course! This module focuses on the practical aspects of creating maps in R using various mapping packages. We’ll explore different libraries, learn how to create thematic maps using survey data, customize map styles, and work with different datasets. This module will equip you with the skills to visualize spatial data effectively.

1. Introduction to Different R Mapping Libraries

R offers a variety of packages for creating maps, each with its own strengths and use cases. We will focus on the following:

2. Basic Maps with maptools and geodata

We’ll start by using the maptools and geodata packages, which can get us some basic maps.

2.1 Using maptools

# Load the necessary library
library(maptools)
## Loading required package: sp
## Please note that 'maptools' will be retired during October 2023,
## plan transition at your earliest convenience (see
## https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
## for guidance);some functionality will be moved to 'sp'.
##  Checking rgeos availability: FALSE
## 
## Attaching package: 'maptools'
## The following object is masked from 'package:sp':
## 
##     sp2Mondrian
# Load the wrld_simpl dataset, which has country outlines
data(wrld_simpl)

# Look at the country names
print(wrld_simpl$NAME)
##   [1] Antigua and Barbuda                      
##   [2] Algeria                                  
##   [3] Azerbaijan                               
##   [4] Albania                                  
##   [5] Armenia                                  
##   [6] Angola                                   
##   [7] American Samoa                           
##   [8] Argentina                                
##   [9] Australia                                
##  [10] Bahrain                                  
##  [11] Barbados                                 
##  [12] Bermuda                                  
##  [13] Bahamas                                  
##  [14] Bangladesh                               
##  [15] Belize                                   
##  [16] Bosnia and Herzegovina                   
##  [17] Bolivia                                  
##  [18] Burma                                    
##  [19] Benin                                    
##  [20] Solomon Islands                          
##  [21] Brazil                                   
##  [22] Bulgaria                                 
##  [23] Brunei Darussalam                        
##  [24] Canada                                   
##  [25] Cambodia                                 
##  [26] Sri Lanka                                
##  [27] Congo                                    
##  [28] Democratic Republic of the Congo         
##  [29] Burundi                                  
##  [30] China                                    
##  [31] Afghanistan                              
##  [32] Bhutan                                   
##  [33] Chile                                    
##  [34] Cayman Islands                           
##  [35] Cameroon                                 
##  [36] Chad                                     
##  [37] Comoros                                  
##  [38] Colombia                                 
##  [39] Costa Rica                               
##  [40] Central African Republic                 
##  [41] Cuba                                     
##  [42] Cape Verde                               
##  [43] Cook Islands                             
##  [44] Cyprus                                   
##  [45] Denmark                                  
##  [46] Djibouti                                 
##  [47] Dominica                                 
##  [48] Dominican Republic                       
##  [49] Ecuador                                  
##  [50] Egypt                                    
##  [51] Ireland                                  
##  [52] Equatorial Guinea                        
##  [53] Estonia                                  
##  [54] Eritrea                                  
##  [55] El Salvador                              
##  [56] Ethiopia                                 
##  [57] Austria                                  
##  [58] Czech Republic                           
##  [59] French Guiana                            
##  [60] Finland                                  
##  [61] Fiji                                     
##  [62] Falkland Islands (Malvinas)              
##  [63] Micronesia, Federated States of          
##  [64] French Polynesia                         
##  [65] France                                   
##  [66] Gambia                                   
##  [67] Gabon                                    
##  [68] Georgia                                  
##  [69] Ghana                                    
##  [70] Grenada                                  
##  [71] Greenland                                
##  [72] Germany                                  
##  [73] Guam                                     
##  [74] Greece                                   
##  [75] Guatemala                                
##  [76] Guinea                                   
##  [77] Guyana                                   
##  [78] Haiti                                    
##  [79] Honduras                                 
##  [80] Croatia                                  
##  [81] Hungary                                  
##  [82] Iceland                                  
##  [83] India                                    
##  [84] Iran (Islamic Republic of)               
##  [85] Israel                                   
##  [86] Italy                                    
##  [87] Cote d'Ivoire                            
##  [88] Iraq                                     
##  [89] Japan                                    
##  [90] Jamaica                                  
##  [91] Jordan                                   
##  [92] Kenya                                    
##  [93] Kyrgyzstan                               
##  [94] Korea, Democratic People's Republic of   
##  [95] Kiribati                                 
##  [96] Korea, Republic of                       
##  [97] Kuwait                                   
##  [98] Kazakhstan                               
##  [99] Lao People's Democratic Republic         
## [100] Lebanon                                  
## [101] Latvia                                   
## [102] Belarus                                  
## [103] Lithuania                                
## [104] Liberia                                  
## [105] Slovakia                                 
## [106] Liechtenstein                            
## [107] Libyan Arab Jamahiriya                   
## [108] Madagascar                               
## [109] Martinique                               
## [110] Mongolia                                 
## [111] Montserrat                               
## [112] The former Yugoslav Republic of Macedonia
## [113] Mali                                     
## [114] Morocco                                  
## [115] Mauritius                                
## [116] Mauritania                               
## [117] Malta                                    
## [118] Oman                                     
## [119] Maldives                                 
## [120] Mexico                                   
## [121] Malaysia                                 
## [122] Mozambique                               
## [123] Malawi                                   
## [124] New Caledonia                            
## [125] Niue                                     
## [126] Niger                                    
## [127] Aruba                                    
## [128] Anguilla                                 
## [129] Belgium                                  
## [130] Hong Kong                                
## [131] Northern Mariana Islands                 
## [132] Faroe Islands                            
## [133] Andorra                                  
## [134] Gibraltar                                
## [135] Isle of Man                              
## [136] Luxembourg                               
## [137] Macau                                    
## [138] Monaco                                   
## [139] Palestine                                
## [140] Montenegro                               
## [141] Mayotte                                  
## [142] Aaland Islands                           
## [143] Norfolk Island                           
## [144] Cocos (Keeling) Islands                  
## [145] Antarctica                               
## [146] Bouvet Island                            
## [147] French Southern and Antarctic Lands      
## [148] Heard Island and McDonald Islands        
## [149] British Indian Ocean Territory           
## [150] Christmas Island                         
## [151] United States Minor Outlying Islands     
## [152] Vanuatu                                  
## [153] Nigeria                                  
## [154] Netherlands                              
## [155] Norway                                   
## [156] Nepal                                    
## [157] Nauru                                    
## [158] Suriname                                 
## [159] Nicaragua                                
## [160] New Zealand                              
## [161] Paraguay                                 
## [162] Peru                                     
## [163] Pakistan                                 
## [164] Poland                                   
## [165] Panama                                   
## [166] Portugal                                 
## [167] Papua New Guinea                         
## [168] Guinea-Bissau                            
## [169] Qatar                                    
## [170] Reunion                                  
## [171] Romania                                  
## [172] Republic of Moldova                      
## [173] Philippines                              
## [174] Puerto Rico                              
## [175] Russia                                   
## [176] Rwanda                                   
## [177] Saudi Arabia                             
## [178] Saint Kitts and Nevis                    
## [179] Seychelles                               
## [180] South Africa                             
## [181] Lesotho                                  
## [182] Botswana                                 
## [183] Senegal                                  
## [184] Slovenia                                 
## [185] Sierra Leone                             
## [186] Singapore                                
## [187] Somalia                                  
## [188] Spain                                    
## [189] Saint Lucia                              
## [190] Sudan                                    
## [191] Sweden                                   
## [192] Syrian Arab Republic                     
## [193] Switzerland                              
## [194] Trinidad and Tobago                      
## [195] Thailand                                 
## [196] Tajikistan                               
## [197] Tokelau                                  
## [198] Tonga                                    
## [199] Togo                                     
## [200] Sao Tome and Principe                    
## [201] Tunisia                                  
## [202] Turkey                                   
## [203] Tuvalu                                   
## [204] Turkmenistan                             
## [205] United Republic of Tanzania              
## [206] Uganda                                   
## [207] United Kingdom                           
## [208] Ukraine                                  
## [209] United States                            
## [210] Burkina Faso                             
## [211] Uruguay                                  
## [212] Uzbekistan                               
## [213] Saint Vincent and the Grenadines         
## [214] Venezuela                                
## [215] British Virgin Islands                   
## [216] Viet Nam                                 
## [217] United States Virgin Islands             
## [218] Namibia                                  
## [219] Wallis and Futuna Islands                
## [220] Samoa                                    
## [221] Swaziland                                
## [222] Yemen                                    
## [223] Zambia                                   
## [224] Zimbabwe                                 
## [225] Indonesia                                
## [226] Guadeloupe                               
## [227] Netherlands Antilles                     
## [228] United Arab Emirates                     
## [229] Timor-Leste                              
## [230] Pitcairn Islands                         
## [231] Palau                                    
## [232] Marshall Islands                         
## [233] Saint Pierre and Miquelon                
## [234] Saint Helena                             
## [235] San Marino                               
## [236] Turks and Caicos Islands                 
## [237] Western Sahara                           
## [238] Serbia                                   
## [239] Holy See (Vatican City)                  
## [240] Svalbard                                 
## [241] Saint Martin                             
## [242] Saint Barthelemy                         
## [243] Guernsey                                 
## [244] Jersey                                   
## [245] South Georgia South Sandwich Islands     
## [246] Taiwan                                   
## 246 Levels: Aaland Islands Afghanistan Albania Algeria ... Zimbabwe
# Extract the data for Somalia
somalia <- wrld_simpl[wrld_simpl$NAME == "Somalia", ]

# Print the regions and subregions information of Somalia
print(somalia$REGION)
## [1] 2
print(somalia$SUBREGION)
## [1] 14
# Extract the data for Djibouti to compare
Djibouti<-wrld_simpl[wrld_simpl$NAME=="Djibouti",]

# Print the regions and subregions information of Djibouti
print(Djibouti$SUBREGION)
## [1] 14
print(Djibouti$REGION)
## [1] 2
# Plot Somalia 

plot(somalia, main = "Somalia Outline")

# Plot Djibouti 

plot(Djibouti, main = "Djibouti Outline")

# We can change the color and appearance:

plot(somalia, border = "black", col = "blue", main = "Somalia: Administrative Boundaries")

plot(Djibouti,border="red",col="green", main="Djibouti Regions")

2.2 Using geodata

Now, let’s get more detailed maps using the geodata package. This will provide regional and district level maps.

# Load the necessary libraries
library(raster)
library(geodata)
## Loading required package: terra
## terra 1.7.71
# Get Somalia's country-level map
Somalia1 <- getData('GADM', country = 'SO', level = 0)
## Warning in getData("GADM", country = "SO", level = 0): getData will be removed in a future version of raster
## . Please use the geodata package instead
plot(Somalia1, main = "Somalia Country Level")

# Get Somalia's regional map
Somalia2 <- getData('GADM', country = 'SO', level = 1)
## Warning in getData("GADM", country = "SO", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package instead
plot(Somalia2, main = "Somalia Regional Level")

# Get Somalia's district map
Somalia3 <- getData('GADM', country = 'SO', level = 2)
## Warning in getData("GADM", country = "SO", level = 2): getData will be removed in a future version of raster
## . Please use the geodata package instead
plot(Somalia3, main = "Somalia District Level")

2.3 Plotting with Names

Now we are going to plot the maps we have and add labels to the regions and districts.

# Print region names
print(Somalia2$NAME_1)
##  [1] "Awdal"             "Bakool"            "Banaadir"         
##  [4] "Bari"              "Bay"               "Galguduud"        
##  [7] "Gedo"              "Hiiraan"           "Jubbada Dhexe"    
## [10] "Jubbada Hoose"     "Mudug"             "Nugaal"           
## [13] "Sanaag"            "Shabeellaha Dhexe" "Shabeellaha Hoose"
## [16] "Sool"              "Togdheer"          "Woqooyi Galbeed"
# Plot Somalia with region names

plot(Somalia2, main = "Somalia Regions with Names")
text(Somalia2, labels = Somalia2$NAME_1, col = "black", cex = 0.8)

# Print district names
print(Somalia3$NAME_2)
##  [1] "Baki"          "Boorama"       "Lughaya"       "Zeylac"       
##  [5] "Ceel Barde"    "Rab Dhuure"    "Tiyeeglow"     "Wajid"        
##  [9] "Xudur"         "Mogadisho"     "Bander-Beyla"  "Bosaaso"      
## [13] "Calawla"       "Iskushuban"    "Qandala"       "Qardho"       
## [17] "Baydhabo"      "Buur Xakaba"   "Diinsoor"      "Qansax Dheere"
## [21] "Caabudwaaq"    "Cadaado"       "Ceel Buur"     "Ceel Dheer"   
## [25] "Dhuusamareeb"  "Baar-Dheere"   "Beled Xaawo"   "Ceel Waaq"    
## [29] "Dolow"         "Garbahaaray"   "Luuk"          "Beled Weyn"   
## [33] "Buulo Burdo"   "Jalalaqsi"     "Bu'aale"       "Jilib"        
## [37] "Saakow"        "Afmadow"       "Badhaadhe"     "Jamaame"      
## [41] "Kismaayo"      "Gaalkacayo"    "Goldogob"      "Hobyo"        
## [45] "Jariiban"      "Xarardheere"   "Burtinle"      "Eyl"          
## [49] "Garoowe"       "Badhan"        "Ceel-Afwein"   "Ceerigaabo"   
## [53] "Aadan"         "Balcad"        "Cadale"        "Jawhar"       
## [57] "Afgooye"       "Baraawe"       "Kuntuwaaray"   "Marka"        
## [61] "Qoryooley"     "Sablale"       "Wanla Weyn"    "Caynabo"      
## [65] "Lascaanod"     "Taleex"        "Xudun"         "Burao"        
## [69] "Buuhoodle"     "Oodweyne"      "Sheekh"        "Berbera"      
## [73] "Gabiley"       "Hargeysa"
# Plot Somalia with district names

plot(Somalia3, main = "Somalia Districts with Names")
text(Somalia3, labels = Somalia3$NAME_2, col = "black", cex = 0.5)

2.4 Plotting with Colors

We can add colors to distinguish the regions and districts.

# Color regions
region_colors <- rainbow(length(Somalia2$NAME_1))

plot(Somalia2, col = region_colors, main = "Somalia Regions Colored")
text(Somalia2, labels = Somalia2$NAME_1, col = "black", cex = 0.8)

# Color districts
district_colors <- rainbow(length(Somalia3$NAME_2))

plot(Somalia3, col = district_colors, main = "Somalia Districts Colored")
text(Somalia3, labels = Somalia3$NAME_2, col = "black", cex = 0.5)

2.5 Highlighting specific regions

Let us say we are interested in a specific part of Somalia. Let’s say we want to highlight the areas that are part of the region that is known as Somaliland. Here we will also change the colors to improve visualization.

# Define the regions you want to highlight for Somaliland
selected_regions <- c("Awdal", "Woqooyi Galbeed", "Togdheer", "Sool", "Sanaag")  # Replace with the regions you want

# Select only these regions from the shapefile
selected_Somalia <- Somalia2[Somalia2$NAME_1 %in% selected_regions, ]

# Define a color palette for the regions
region_colors <- rainbow(length(selected_Somalia$NAME_1))

# Plot the selected regions

plot(selected_Somalia, col = region_colors, main = "Somaliland Regions")
text(selected_Somalia, labels = selected_Somalia$NAME_1, col = "black", cex = 0.8)

Now, we want to do the same but this time selecting the Southwest State in Somalia.

# Define the regions you want to highlight for Southwest State
selected_regions1 <- c("Bay", "Bakool", "Shabeellaha Hoose")

# Select only these regions from the shapefile
selected_Somalia1 <- Somalia2[Somalia2$NAME_1 %in% selected_regions1, ]
# Define a color palette for the regions
region_colors1<- rainbow(length(selected_Somalia1$NAME_1))

# Plot the selected regions
plot(selected_Somalia1, col = region_colors1, main = "South West State Regions")
text(selected_Somalia1, labels = selected_Somalia1$NAME_1, col = "black", cex = 0.8)

3. Working with Shapefiles (sf package)

Now, we move to the sf package, which allows us to work directly with shapefiles. Shapefiles are a common way to store geospatial data.

3.1 Loading Shapefiles

Let’s use some shapefiles that we have in a folder. You will need to change the path to the folders in your computer. We will load Admin0 (Whole Country Data), Admin1 (Regional Data) and Admin2 (District Data) and then we will visualize them.

# Load the necessary library
library(sf)
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
# Set the file path to your shapefile - **Change this to your file paths**
filepath0 <- "~/gadm41_SOM_0.shp"
data0 <- st_read(filepath0)
## Reading layer `gadm41_SOM_0' from data source 
##   `C:\Users\Admin\Documents\gadm41_SOM_0.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 40.9785 ymin: -1.647082 xmax: 51.4157 ymax: 11.98931
## Geodetic CRS:  WGS 84
# Plot the shapefile of the whole country

plot(data0, main ="Admin 0: Whole Country")

filepath1 <- "~/gadm41_SOM_1.shp"
data1 <- st_read(filepath1)
## Reading layer `gadm41_SOM_1' from data source 
##   `C:\Users\Admin\Documents\gadm41_SOM_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 18 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 40.9785 ymin: -1.647082 xmax: 51.4157 ymax: 11.98931
## Geodetic CRS:  WGS 84
# Plot the shapefile of the regions

plot(data1, main ="Admin 1: Regions")
## Warning: plotting the first 10 out of 11 attributes; use max.plot = 11 to plot
## all

filepath2 <- "~/gadm41_SOM_2.shp"
data2 <- st_read(filepath2)
## Reading layer `gadm41_SOM_2' from data source 
##   `C:\Users\Admin\Documents\gadm41_SOM_2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 74 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 40.9785 ymin: -1.647082 xmax: 51.4157 ymax: 11.98931
## Geodetic CRS:  WGS 84
# Plot the shapefile of the districts
plot(data2, main = "Admin 2: Districts")
## Warning: plotting the first 10 out of 13 attributes; use max.plot = 13 to plot
## all

Let us load a shapefile that we called admin2.shx and another one that has the shapefile of Somaliland SL_SixRegions.shx and visualize them.

# Import shape files
library(sf)

# Set the file path to your shapefile - **Change this to your file paths**
filepath <- "C:/Users/Admin/Documents/scopus/Shape File/admin2.shx"
admin2 <- st_read(filepath)
## Reading layer `admin2' from data source 
##   `C:\Users\Admin\Documents\scopus\Shape File\admin2.shx' using driver `ESRI Shapefile'
## Simple feature collection with 91 features and 0 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 40.99488 ymin: -1.664897 xmax: 51.41303 ymax: 11.9852
## CRS:           NA
# Visualize the imported Shape File

plot(admin2, main = "Districts")

somaliland <- "C:/Users/Admin/Documents/scopus/SL Shape File/SL_SixRegions.shx"
SL <- st_read(somaliland)
## Reading layer `SL_SixRegions' from data source 
##   `C:\Users\Admin\Documents\scopus\SL Shape File\SL_SixRegions.shx' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 42.68301 ymin: 8 xmax: 49.08111 ymax: 11.51147
## Geodetic CRS:  WGS 84
# check the area of each region
print(SL$Shape_Area)
## [1] 0.8471137 1.1290371 1.3226062 0.2547984 0.9496473 2.1433153
# Visualize the imported Shape File

plot(SL, main = "Somaliland Regions")

area <- SL$Shape_Area

plot(area, main = "Somaliland Region Areas")

4. Interactive Maps with mapview

mapview makes it super easy to create interactive maps that you can zoom, pan, and explore.

# Load the mapview package
library(mapview)

# Map the whole country shapefile
mapview(data0, main = "Whole Somalia Map")
# Map the regional shapefile
mapview(data1, main ="Somalia Regions Map")
# Map the districts shapefile
mapview(data2, main ="Somalia Districts Map")
# Map the Somaliland regions
mapview(SL, main = "Somaliland Regions Map")

5. Mapping with ggplot2

ggplot2 is a powerful package for creating static maps.

# Load the necessary packages
library(ggplot2)
library(viridis)
## Loading required package: viridisLite
# Create a basic map of districts using ggplot2
ggplot(data2) +
  geom_sf(aes(fill = NAME_2)) + # FID is a unique identifier that we will use as a test
   scale_fill_viridis(option = "magma", discrete = TRUE) + # Use discrete viridis scale
    labs(title = "Districts of Somalia Map using ggplot2")

6. Mapping with leaflet

leaflet is a popular package for creating interactive web maps.

# Load the leaflet package
library(leaflet)

# Transform the data to have a CRS of 4326
data1_4326 <- st_transform(data1, 4326)

# Create a color palette using leaflet
pal <- colorNumeric(palette = "YlOrRd", domain = 1:nrow(data1_4326))

# Create the interactive map using leaflet
l <- leaflet(data1_4326) %>%
  addTiles() %>%
  addPolygons(color = "white", fillColor = ~pal(1:nrow(data1_4326)), fillOpacity = 0.8) %>%
  addLegend(pal = pal, values = 1:nrow(data1_4326), opacity = 0.8)

# Add a mini map
l <- l %>% addMiniMap()

# Print the map
print(l)

# Save the map as an HTML file
library(htmlwidgets)
saveWidget(widget = l, file = "leaflet_map.html")

# Save the map as a PNG image (requires webshot or webshot2)
# library(webshot)
# webshot::install_phantomjs() # Install if you don't have it
# webshot(url = "leaflet_map.html", file = "leaflet_map.png")

# Alternatively, use webshot2
# library(webshot2)
# webshot2(url = "leaflet_map.html", file = "leaflet_map.png")

7. Side-by-side maps and Synchronized maps

Here we are going to learn how to compare maps using mapview and how to synchronize them using leafsync.

7.1 Side-by-Side maps with mapview

We will visualize the regions of Somalia data1 and we will also use the admin2 shapefile that we loaded before. We will compare them using mapview.

# load the necessary libraries
library(leaflet.extras2)
library(RColorBrewer)

# Define color palette
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))

# Create a sequence for the breaks
at <- seq(1, max(c(nrow(data1), nrow(admin2))), length.out = 8)

# Create the first map
m1 <- mapview(data1, zcol = "GID_1", map.types = "CartoDB.Positron",
              col.regions = pal, at = at)
## Warning in vectorColRegions(x = x, zcol = zcol, col.regions = col.regions, :
## ignoring 'at' which is only supported for numeric values
# Create the second map
m2 <- mapview(data2, zcol = "GID_2", map.types = "CartoDB.Positron",
              col.regions = pal, at = at)
## Warning in vectorColRegions(x = x, zcol = zcol, col.regions = col.regions, :
## ignoring 'at' which is only supported for numeric values
# Create the side-by-side map
m1 | m2

7.2 Synchronized maps with leafsync

Now we will use leafsync to visualize the same maps, but in this case, the zoom and pan will be synchronized.

# Load the necessary library
library(leafsync)

# Define color palette
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))

# Create a sequence for the breaks
at <- seq(1, max(c(nrow(data1), nrow(data2))), length.out = 8)

# Create the first map
m1 <- mapview(data1, zcol = "GID_1", map.types = "CartoDB.Positron",
              col.regions = pal, at = at)
## Warning in vectorColRegions(x = x, zcol = zcol, col.regions = col.regions, :
## ignoring 'at' which is only supported for numeric values
# Create the second map
m2 <- mapview(data2, zcol = "GID_2", map.types = "CartoDB.Positron",
              col.regions = pal, at = at)
## Warning in vectorColRegions(x = x, zcol = zcol, col.regions = col.regions, :
## ignoring 'at' which is only supported for numeric values
# Synchronize the two maps
m <- leafsync::sync(m1, m2)
m

8. Mapping with tmap

tmap can create static or interactive maps with different styles. We will use the districts map for this section.

# Load the tmap library
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
# Set to plot mode for a static map
tmap_mode("plot")
## tmap mode set to plotting
# Create a map of Somalia districts
tm_shape(data2) +
  tm_polygons("GID_2", title ="Districts of Somalia Using tmap")
## Warning: Number of levels of the variable "GID_2" is 74, which is larger than
## max.categories (which is 30), so levels are combined. Set
## tmap_options(max.categories = 74) in the layer function to show all levels.
## Legend labels were too wide. The labels have been resized to 0.63, 0.60, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.57, 0.60, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

9. Creating Thematic Maps using Survey Data

Thematic maps display statistical data over geographic areas. We can use survey data to create thematic maps, for example, mapping the prevalence of a disease.

9.1 Preparing Survey Data

  • Aggregate Data: Aggregate individual-level survey data to the desired spatial units (e.g., regions, districts).
  • Calculate Rates/Proportions: Calculate the rates or proportions of interest (e.g., disease prevalence, vaccination coverage) for each spatial unit.
  • Merge with Spatial Data: Merge the aggregated survey data with the corresponding spatial data (shapefile) using a common identifier.

9.2 Creating Thematic Maps

  • ggplot2: Use geom_sf() with the fill aesthetic to map the survey data.
  • leaflet: Use addPolygons() with the fillColor aesthetic to map the survey data.
  • tmap: Use tm_polygons() with the fill argument to map the survey data.

9.3 Example: Mapping Unmet Need for Family Planning

Let’s assume we have survey data on unmet need for family planning aggregated at the regional level.

# Load the necessary libraries
library(sf)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
## 
##     intersect, union
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(viridis)
library(leaflet)
library(tmap)

# Load the shapefile for Somalia regions
filepath1 <- "~/gadm41_SOM_1.shp"
data1 <- st_read(filepath1)
## Reading layer `gadm41_SOM_1' from data source 
##   `C:\Users\Admin\Documents\gadm41_SOM_1.shp' using driver `ESRI Shapefile'
## Simple feature collection with 18 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 40.9785 ymin: -1.647082 xmax: 51.4157 ymax: 11.98931
## Geodetic CRS:  WGS 84
# Create simulated survey data
set.seed(123)
survey_data <- data.frame(
  NAME_1 = data1$NAME_1,
  unmet_need = sample(10:40, length(data1$NAME_1), replace = TRUE)
)

# Merge survey data with spatial data
data1_merged <- left_join(data1, survey_data, by = "NAME_1")

# Create a thematic map using ggplot2
ggplot(data1_merged) +
  geom_sf(aes(fill = unmet_need)) +
  scale_fill_viridis(option = "magma") +
  labs(title = "Unmet Need for Family Planning by Region", fill = "Unmet Need (%)")

# Create a thematic map using leaflet
pal <- colorNumeric(palette = "YlOrRd", domain = data1_merged$unmet_need)
leaflet(data1_merged) %>%
  addTiles() %>%
  addPolygons(fillColor = ~pal(unmet_need), color = "white", fillOpacity = 0.7) %>%
  addLegend(pal = pal, values = ~unmet_need, title = "Unmet Need (%)")
# Create a thematic map using tmap
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(data1_merged) +
  tm_polygons("unmet_need", title = "Unmet Need (%)", palette = "YlOrRd")

10. Customizing Map Styles and Layouts

Customizing maps involves adjusting colors, labels, legends, and other elements to enhance clarity and aesthetics.

10.1 Using ggplot2

  • Colors: Use scale_fill_viridis(), scale_fill_brewer(), or custom color palettes.
  • Labels: Use labs() to add titles, subtitles, and axis labels.
  • Themes: Use theme_minimal(), theme_bw(), or custom themes to adjust the overall appearance.

10.2 Using leaflet

  • Tiles: Use different base map providers with addTiles().
  • Popups: Add interactive popups with addPopups().
  • Legends: Customize legends with addLegend().

10.3 Using tmap

  • Styles: Use tm_style() to apply predefined styles.
  • Layouts: Use tm_layout() to adjust map elements.
  • Palettes: Use tm_fill() with different color palettes.

11. Mobility Flows with flowmapblue

The flowmapblue package can be used to visualize mobility data between locations.

# Install flowmapblue from GitHub if you haven't already
# devtools::install_github("FlowmapBlue/flowmapblue.R")

# Load the flowmapblue package
library(flowmapblue)

# Create the locations data frame
locations <- data.frame(
  id = c(1, 2, 3, 4),
  name = c("Mogadishu", "Boorama", "Hargeisa", "Berbera"),
  lat = c(2.037, 9.936, 9.56, 10.44),
  lon = c(45.344, 43.183, 44.065, 45.014)
)

# Create the flows data frame (example data)
flows <- data.frame(
  origin = c(1, 2, 3, 4, 1, 3, 2, 4),
  dest = c(2, 1, 4, 3, 3, 1, 4, 2),
  count = c(150, 120, 80, 90, 100, 70, 60, 50)
)

# Set your Mapbox access token (replace with your actual token)
mapbox_token <- "YOUR_MAPBOX_ACCESS_TOKEN"

# Create the interactive flow map
flowmapblue(locations, flows, mapboxAccessToken = mapbox_token,
            clustering = TRUE, darkMode = TRUE, animation = FALSE)
# If you don't have a mapbox token, you can use this to visualize the flows without a background map
# flowmapblue(locations, flows, mapboxAccessToken = NULL,
#Clustering = TRUE, darkMode = TRUE, animation = FALSE)

12. Maps of Point Data

In addition to areal data, we can also map point data, such as the locations of cities. We will use data from the maps package and focus on cities in Somalia.

# Load the necessary libraries
library(maps)
## 
## Attaching package: 'maps'
## The following object is masked from 'package:viridis':
## 
##     unemp
library(sf)
library(ggplot2)
library(viridis)
library(leaflet)
library(mapview)
library(tmap)

# Load world cities data
d <- world.cities

# Select Somalia
d <- d[which(d$country.etc == "Somalia"), ]

# Transform data to sf object
d <- st_as_sf(d, coords = c("long", "lat"))

# Assign CRS
st_crs(d) <- 4326

# Create variables for population and size
d$vble <- d$pop
d$size <- sqrt(d$vble) / 100

# Map with ggplot2
ggplot(d) +
  geom_sf(aes(col = vble, size = size)) +
  scale_color_viridis() +
  labs(title = "Somali Cities Map using ggplot2")

# Map with leaflet
pal <- colorNumeric(palette = "viridis", domain = d$vble)
print(leaflet(d) %>%
  addTiles() %>%
  addCircles(lng = st_coordinates(d)[, 1],
             lat = st_coordinates(d)[, 2],
             radius = ~sqrt(vble) * 10,
             color = ~pal(vble), popup = ~name) %>%
  addLegend(pal = pal, values = ~vble, position = "bottomright"))

# Map with mapview
mapview(d, zcol = "vble", cex = "size", main ="Somali Cities Map using mapview")
# Map with tmap
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(d) + tm_dots("vble", scale = sqrt(d$vble) / 500,
                      palette = "viridis", title = "Somali Cities Map using tmap")
## Warning in g$scale * (x_legend/maxX)^scaling: longer object length is not a
## multiple of shorter object length

Conclusion

This module has provided a comprehensive introduction to mapping in R using various packages. We covered static and interactive maps, thematic mapping with survey data, map customization, and visualization of point data and mobility flows. You are now equipped with the skills to create effective and informative maps for spatial analysis.

In the next module, we will explore R packages for spatial analysis, building on the mapping skills you have acquired in this module. ```