Mapping in R Webinar

WCA R-EiE WG

Agenda

  • Introduction (10 min)
  • Overview of R software (15 min)
  • Overview of R packages for mapping (15 min)
  • Reading and manipulating spatial/geographic data (15 min))
  • Creating maps (static, interactive)

– Static maps with packages:ggplot2, tmap, mapsf (45 min)

– Interactive maps with packages:Leaflet, tmap, mapview (45 min)

  • Questions & Answers session (30 min)
  • Conclusion (5 min)

Introduction

Overview on GIS

  • GIS

Geographic Information Systems (GIS) store, analyze, and visualize data for geographic positions on Earth’s surface. The four major features of GIS are listed below.

-Create geographic data.

-Manage spatial data in a database.

-Analyze and find patterns.

-Map Visualization.

Introduction

Overview on GIS

  • Well-Known Text (WKT)

Well-known text (WKT) is a text markup language for representing vector geometry objects. The common type of geometric objects are marked as shown in the table below.

Introduction

Overview on GIS

  • Shapefile

The shapefile format is a universal geospatial vector data format.The shapefile format can illustrate the vector features: points, lines, and polygons.Each file is briefly introduced as the following:

Files Features
.shp shape format; the feature geometry
.dbf attribute format; attributes for each shape, stored as two-dimensional table
.prj projection description, using a well-known text representation of coordinate reference systems
.shx shape index format; a positional index of the feature geometry to allow seeking forwards and backwards quickly

Introduction

Overview on GIS

  • Coordinate Reference System (CRS)

-Geometric data is not geospatial unless it is accompanied by CRS information, which allows GIS to display and operate the data accurately.It includes two major components, datum and projection.

-Datum is a model of the shape of the earth. It has angular units (degrees) and defines the starting point (0,0), and hence the coordinate can represent a specific spot on the earth.

-Projection is a mathematical transformation of the angular measurements on a round earth to a flat surface. The units associated with a given projection are usually linear (feet, meters, etc.).

Introduction

Overview on GIS

  • Coordinate Reference System (CRS)
Type Features Example
Geographic Coordinate Systems A global or spherical coordinate system such as latitude–longitude. EPSG:4326
(World Geodetic System 1984)
Projected Coordinate System Based on a map projection such as transverse Mercator, Albers equal area, or Robinson, which project the spherical surface onto a two-dimensional coordinate plane. EPSG:3826
(TWD97 / TM2 zone 121, Taiwan Datum 1997)

Introduction

The problem for cartographers

Many maps produced in humanitarian context are made with a set of software products that slices the cartographic process.

A solution

R Simplifies the map making process in an unified workflow thanks to its spatial ecosystem !

Overview of R software

R is a popular programming language used for statistical computing and graphical presentation.

R is derived form the S language developed at ATT Bell Laboratories.

R was originally developed for teaching at the University of Auckland, New Zealand, by Ross Ihaka and Robert Gentleman.

-It is a great resource for data analysis, data visualization, data science and machine learning

-It provides many statistical techniques

-It is easy to draw graphs in R,

-It works on different platforms

-RStudio/Posit is an integrated development environment (IDE) for R.

-It has many packages (libraries of functions)

-It is open-source and free

-It has a large community support

Overview of R packages for mapping (15 min)

Spatial operations

  • before 2003: spatial, sgeostat, splancs, akima, geoR, spatstat, spdep, maptools.

  • 2003: rgdal (Roger Bivand, Tim Keitt & Barry Rowlingson), interface between R and GDAL and PROJ4

  • 2005: sp (Edzer Pebesma & Roger Bivand), classes methods for spatial objects, quickly adopted.

  • 2008: sp surpport in ggplot2 (Hadley Wickham)

  • 2010: rgeos (Roger Bivand & Colin Rundel), interface between R & GEOS.

  • 2010: raster (Robert J. Hijmans), support for raster data

  • 2016: sf (Edzer Pebesma), supersedes sp, rgdal & rgeos

  • 2018: stars (Edzer Pebesma), supersedes raster.

  • 2020: terra (Robert J. Hijmans) also supersedes raster.

Aperçu des packages de R pour la cartographie (15 min)

Thematic mapping

  • 2014: tmap (Martijn Tennekes)

  • 2017: ggplot2 + ggspatial (Dewey Dunnington)

  • 2021: mapsf (T. Giraud) supersedes cartography.

Overview of R packages for mapping (15 min)

Interactive mapping

  • 2015: leaflet (Joe Cheng et al.), relies on the leaflet javascript library.

  • 2015: mapview (Tim Appelhans et al.), relies on the leaflet package.

  • 2018: mapdeck (David Cooley), relies on Mapbox GL & Deck.gl libraries.

Reading and manipulating spatial/geographic data

sf, corner stone of the R spatial ecosystem

Interface between R and extensively used libraries in GIS:

  • GDAL - Geospatial Data Abstraction Library
  • PROJ - Coordinate Transformation Software
  • GEOS - Geometry Engine - Open Source

Pebesma and Bivand (2023)

Reading and manipulating spatial/geographic data

The sf object does not only contain the geometry information, but sets of features with attributes. The basic component of its geometry data is sfg, which is recorded in WKT format.

A simple feature must contain feature attributes and a simple feature geometries (sfg) object, which defines the location of the tuple. The list-column with the simple feature geometries (sfg) for each tuple is then called simple feature columns (sfc).

expand for full code 
library(sf)
ls("package:sf")
  [1] "%>%"                          "as_Spatial"                  
  [3] "dbDataType"                   "dbWriteTable"                
  [5] "gdal_addo"                    "gdal_create"                 
  [7] "gdal_crs"                     "gdal_extract"                
  [9] "gdal_inv_geotransform"        "gdal_metadata"               
 [11] "gdal_polygonize"              "gdal_rasterize"              
 [13] "gdal_read"                    "gdal_read_mdim"              
 [15] "gdal_subdatasets"             "gdal_utils"                  
 [17] "gdal_write"                   "gdal_write_mdim"             
 [19] "get_key_pos"                  "NA_agr_"                     
 [21] "NA_bbox_"                     "NA_crs_"                     
 [23] "NA_m_range_"                  "NA_z_range_"                 
 [25] "pivot_wider.sf"               "plot_sf"                     
 [27] "rawToHex"                     "read_sf"                     
 [29] "sf.colors"                    "sf_add_proj_units"           
 [31] "sf_extSoftVersion"            "sf_proj_info"                
 [33] "sf_proj_network"              "sf_proj_pipelines"           
 [35] "sf_proj_search_paths"         "sf_project"                  
 [37] "sf_use_s2"                    "st_agr"                      
 [39] "st_agr<-"                     "st_area"                     
 [41] "st_as_binary"                 "st_as_grob"                  
 [43] "st_as_s2"                     "st_as_sf"                    
 [45] "st_as_sfc"                    "st_as_text"                  
 [47] "st_axis_order"                "st_bbox"                     
 [49] "st_bind_cols"                 "st_boundary"                 
 [51] "st_break_antimeridian"        "st_buffer"                   
 [53] "st_can_transform"             "st_cast"                     
 [55] "st_centroid"                  "st_collection_extract"       
 [57] "st_combine"                   "st_concave_hull"             
 [59] "st_contains"                  "st_contains_properly"        
 [61] "st_convex_hull"               "st_coordinates"              
 [63] "st_covered_by"                "st_covers"                   
 [65] "st_crop"                      "st_crosses"                  
 [67] "st_crs"                       "st_crs<-"                    
 [69] "st_delete"                    "st_difference"               
 [71] "st_dimension"                 "st_disjoint"                 
 [73] "st_distance"                  "st_drivers"                  
 [75] "st_drop_geometry"             "st_equals"                   
 [77] "st_equals_exact"              "st_filter"                   
 [79] "st_geometry"                  "st_geometry_type"            
 [81] "st_geometry<-"                "st_geometrycollection"       
 [83] "st_graticule"                 "st_inscribed_circle"         
 [85] "st_interpolate_aw"            "st_intersection"             
 [87] "st_intersects"                "st_is"                       
 [89] "st_is_empty"                  "st_is_longlat"               
 [91] "st_is_simple"                 "st_is_valid"                 
 [93] "st_is_within_distance"        "st_jitter"                   
 [95] "st_join"                      "st_layers"                   
 [97] "st_length"                    "st_line_merge"               
 [99] "st_line_sample"               "st_linestring"               
[101] "st_m_range"                   "st_make_grid"                
[103] "st_make_valid"                "st_minimum_rotated_rectangle"
[105] "st_multilinestring"           "st_multipoint"               
[107] "st_multipolygon"              "st_nearest_feature"          
[109] "st_nearest_points"            "st_node"                     
[111] "st_normalize"                 "st_overlaps"                 
[113] "st_perimeter"                 "st_point"                    
[115] "st_point_on_surface"          "st_polygon"                  
[117] "st_polygonize"                "st_precision"                
[119] "st_precision<-"               "st_read"                     
[121] "st_read_db"                   "st_relate"                   
[123] "st_reverse"                   "st_sample"                   
[125] "st_segmentize"                "st_set_agr"                  
[127] "st_set_crs"                   "st_set_geometry"             
[129] "st_set_precision"             "st_sf"                       
[131] "st_sfc"                       "st_shift_longitude"          
[133] "st_simplify"                  "st_snap"                     
[135] "st_sym_difference"            "st_touches"                  
[137] "st_transform"                 "st_triangulate"              
[139] "st_triangulate_constrained"   "st_union"                    
[141] "st_viewport"                  "st_voronoi"                  
[143] "st_within"                    "st_wrap_dateline"            
[145] "st_write"                     "st_write_db"                 
[147] "st_z_range"                   "st_zm"                       
[149] "vec_cast.sfc"                 "vec_ptype2.sfc"              
[151] "write_sf"

Geometric checking

st_intersects: touch or overlap

st_disjoint: !intersects

st_touches: touch

st_crosses: cross (do not touch)

st_within: within

st_contains: contains

Geometric operations

st_read/st_write: reads a sf object / writes a sf object,

st_buffer: compute a buffer around this geometry/each geometry

st_transform: transforms data to a new CRS,

st_intersection: intersects sf objects,

st_union: combines several sf objects into one,

st_as_sf: converts to a sf object.

library(pacman)
pacman::p_load(
  readxl,         # import/export excel files
  tidyverse,      # includes ggplot2 and other
  rio,            # import/export
  here,           # file locator
  stringr,        # working with characters   
  scales,         # transform numbers
  ggrepel,        # smartly-placed labels
  gghighlight,    # highlight one part of plot
  RColorBrewer,   # color scales
  viridis,        # color scales
  sf,             # simple features for...
  tmap,           # tmap for thematic and interactive map 
  mapsf,          # mapsf for thematic map
  leaflet,        # leaflet for interactive map
  mapview,        # mapview for interactive map
  janitor,        # clean variables names
  ggspatial,      # add map's scale and north arrow
  gganimate,      # graphics animation for ggplot2
  leafsync,      #interactive map synchronisation
  knitr          # Getting nice table
  )
hno <- read_excel(here("data","ner_hno_pin_data_2024.xlsx"))  %>% clean_names() %>% select(id:population_totale,cluster,pin) %>%  filter(cluster=="Education") %>% mutate(admin_2_p_code=str_remove(admin_2_p_code,"R"))
ner_shp_adm2<-st_read(here("data","ner_adm_ignn_20230720_ab_shp","NER_admbnda_adm2_IGNN_20230720.shp"))
Reading layer `NER_admbnda_adm2_IGNN_20230720' from data source 
  `C:\Users\Souleymane SAKANDE\OneDrive - Norwegian Refugee Council\Desktop\Webinar folder\mapping\cartography_r\data\ner_adm_ignn_20230720_ab_shp\NER_admbnda_adm2_IGNN_20230720.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 67 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 0.1662505 ymin: 11.69697 xmax: 15.99698 ymax: 23.52503
Geodetic CRS:  WGS 84
Code 
ner_shp_adm2_hno <-left_join(ner_shp_adm2,hno,by=c("ADM2_PCODE"="admin_2_p_code"))
# ner_shp_adm2<-read_sf()
# ner_gjson_adm1<-st_write(ner_shp_adm1,"./data/NER_admin1.geojson")
# ner_gjson_adm2<-write_sf()
head(ner_shp_adm2,3)
Simple feature collection with 3 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 3.861768 ymin: 13.70453 xmax: 8.087891 ymax: 16.61792
Geodetic CRS:  WGS 84
  ADM2_FR ADM2_PCODE ADM2_REF ADM1_FR ADM1_PCODE     ADM0_FR ADM0_PCODE
1  Abalak   NE005001     <NA>  Tahoua      NE005 Niger (the)         NE
2  Tillia   NE005012     <NA>  Tahoua      NE005 Niger (the)         NE
3  Mayahi   NE004008     <NA>  Maradi      NE004 Niger (the)         NE
        date    validOn validTo Shape_Leng Shape_Area AREA_SQKM
1 2023-06-21 2023-07-25    <NA>   5.540334  1.1761078 13962.561
2 2023-06-21 2023-07-25    <NA>   5.209354  1.4855289 17603.184
3 2023-06-21 2023-07-25    <NA>   3.455926  0.5502204  6573.335
                        geometry
1 MULTIPOLYGON (((6.580872 15...
2 MULTIPOLYGON (((3.861768 15...
3 MULTIPOLYGON (((7.637695 13...
plot(ner_shp_adm2)

plot(ner_shp_adm2$geometry)

plot(st_geometry(ner_shp_adm2))

ner_shp_adm0_u<-st_union(ner_shp_adm2)
plot(st_geometry(ner_shp_adm0_u),border="blue",lwd=2)

ner_shp_adm1<-ner_shp_adm2 %>% group_by(ADM1_FR) %>% count()
plot(st_geometry(ner_shp_adm1),border="orange",lwd=2)

expand for full code 
ner_shp_adm2_u<-st_union(ner_shp_adm2)
ner_shp_adm2_100km<-st_buffer(ner_shp_adm2_u,100000)

#Plot
plot(st_geometry(ner_shp_adm2_100km),border="red",col="grey")
plot(st_geometry(ner_shp_adm2),add=T,col="lightblue")
plot(st_geometry(ner_shp_adm2_u),border="blue",add=T,lwd=2)

expand for full code 
plot(st_geometry(ner_shp_adm2))

expand for full code 
ner_shp_adm1_aga<-ner_shp_adm1 %>% filter(ADM1_FR=="Agadez")
plot(st_geometry(ner_shp_adm1_aga),border="red")

expand for full code 
ner_shp_adm2_aga<-st_intersection(ner_shp_adm2,ner_shp_adm1_aga)
plot(st_geometry(ner_shp_adm1_aga),border="red")
plot(st_geometry(ner_shp_adm2_aga),add=T,col="orange",border="blue")

hno_sf= st_as_sf(hno, coords = c("admin_2_long", "admin_2_lat"), crs = 4326)
plot(st_geometry(hno_sf),col="red",pch = 3)

The st_transform() function is used to convert the existing spatial object’s coordinates into another projection. For example, let’s transform our hno_sf object to the Equal Earth projection (EPSG 8857).

Note

hno_sf_trans = st_transform(hno_sf, 8857)

ner_shp_adm2_ctd<-st_centroid(ner_shp_adm2)
plot(st_geometry(ner_shp_adm2), col = sf.colors(12, categorical = TRUE))
plot(st_geometry(ner_shp_adm2_ctd),col="red",pch = 19,add=T)

Static maps

  • ggplot2
bfa_shp_adm2<-st_read(here("data","bfa_adm_igb_20200323_shp","bfa_admbnda_adm2_igb_20200323.shp"))
Reading layer `bfa_admbnda_adm2_igb_20200323' from data source 
  `C:\Users\Souleymane SAKANDE\OneDrive - Norwegian Refugee Council\Desktop\Webinar folder\mapping\cartography_r\data\bfa_adm_igb_20200323_shp\bfa_admbnda_adm2_igb_20200323.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 45 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -5.511255 ymin: 9.415337 xmax: 2.407427 ymax: 15.08311
Geodetic CRS:  WGS 84
acled <- read_excel(here("data","bfa_acled_21feb2024.xlsx"),sheet="Data") %>% clean_names() %>% 
  group_by(admin2,admin2_pcode,year) %>% summarize(events=sum(events),fatalities=sum(fatalities)) %>% 
  filter(year %in% c("2020","2021","2022","2023"))
bfa_shp_adm2_acled<-bfa_shp_adm2 %>% left_join(acled,by=c("ADM2_PCODE"="admin2_pcode"))
ggplot() +
  geom_sf(data = bfa_shp_adm2_acled, aes(fill = events))

bfa_shp_adm2_acled$eventsClass<- cut (bfa_shp_adm2_acled$events, breaks =c(0,100,300,500,750,1000,Inf),
labels=c('< 100', '100-300', '300-500', '500-750','750-1000','> 1000'))

bfa_m<-ggplot()+
geom_sf(aes(fill=eventsClass), color='grey10',data=bfa_shp_adm2_acled) + 
  scale_fill_brewer (palette="Reds",name='Nb. of events',guide=guide_legend (direction='horizontal',title.position='top',
title.hjust=.5,
label.hjust=.5,
label.position='right',
keywidth=unit(5, units = "mm"),
keyheight=unit(3, units = "mm"),
nrow=1,
byrow=TRUE
))+
labs (title="Burkina Faso's security context",
subtitle='Acled events',
caption=c('Source: ACLED, Feb 2024'))+
theme_void()+
theme (title=element_text(face='bold'),
legend.position='bottom',
plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)
) 

bfa_m + annotation_scale(location="br") +
annotation_north_arrow(location="tl") 

bfa_m + facet_wrap(~year,nrow=1)

bfa_m_anim<-bfa_m + transition_manual(year) + 
  labs(subtitle = "Year:{current_frame}")
animate(bfa_m_anim, nframes =4, duration = 10)

#geom_sf_text(aes(label = NAME), col="white", size=1.5)
points<- st_centroid(bfa_shp_adm2_acled, of_largest_polygon = TRUE)

ggplot() +
  geom_sf(data=bfa_shp_adm2_acled) +
  geom_sf(data =points,pch = 21,aes(size = fatalities, fill = fatalities),col = "grey20") + 
  scale_size(guide = guide_legend(direction = "horizontal",nrow = 1,label.position="bottom")) +
  labs(size="") +
  scale_fill_gradientn(colours = hcl.colors(4, "RdBu",rev = TRUE,alpha = 0.9)) +
  guides(fill = guide_legend(title = ""))  +
  theme_void() +
  theme(legend.position = "top")

Static maps

  • tmap
mli_shp_adm2 <- read_sf(here("data","mli_admb_1m_gov_20211220_shp","mli_admbnda_adm2_1m_gov_20211220.shp"))
mli_Q3_3W<-read_excel(here("data","mali-3w_troisieme-trimestre2023.xlsx")) %>% slice(-1) %>% clean_names()
mli_Q3_3w<-mli_Q3_3W %>% select(organisation_accronyme:cluster) %>% 
group_by(admin2pcod) %>% mutate(n_org_adm2=n_distinct(nom_organisation)) %>% ungroup() %>% 
group_by(admin2pcod) %>% mutate(n_cluster_adm2=n_distinct(cluster)) %>% ungroup() %>% 
group_by(admin2pcod,cluster) %>% mutate(n_org_adm2_clus=n_distinct(nom_organisation)) %>% 
ungroup()
mli_shp_adm2_3w<-mli_shp_adm2 %>% left_join(mli_Q3_3w,by=c("ADM2_PCODE"="admin2pcod"))
qtm(mli_shp_adm2_3w,fill="n_org_adm2")

mli_m<-tm_shape(mli_shp_adm2_3w) +
  tm_polygons(col = "n_org_adm2",style="quantile",pal="-viridis",title="Nb.Org") +
  tm_layout(main.title="Mali operationnal presence",main.title.position ="center",frame=F)
mli_m2<-mli_m +  
  tm_logo(here("data","gec_mali_logo_rgb_0.png"), height = 1) +
  tm_scale_bar(position = c("left", "bottom"), width = 0.1) +
  tm_compass(position = c("right", "top"), size = 1) +
  tm_credits ("Source: Mali's HDX page", align = "right") + 
  tm_text("cercle", size ="AREA",print.tiny = TRUE)
mli_m2

mli_m + tm_facets(by = "ADM1_FR", ncol = 2) +
  tm_layout(legend.outside.size = 0.2)

mli_m3<-mli_m + tm_facets(along = "ADM1_FR",free.coords = F) +
  tm_layout(legend.position = c("left", "bottom")) 
tmap_animation(mli_m3, filename = "mli_adm1.gif", delay = 50, width = 700, height = 400)
#tm_layout(panel.labels = c("label1", "label2")) 

tm_shape(mli_shp_adm2_3w) +
  tm_fill(c("n_org_adm2", "n_cluster_adm2"), palette = "Oranges", title = c("Nb.Org", "Nb.Cluster"),style="quantile",n=4) +
  tm_layout(frame = F)

tm_shape(mli_shp_adm2_3w) + tm_fill("n_cluster_adm2", palette="Blues", title = "Nb.Cluster") +
  tm_borders(alpha=.4) +
  tm_bubbles(size = "n_org_adm2", col = "n_org_adm2",
                                      palette = "YlOrRd", style = "quantile",
                                      title.col = "Nb.Org",
                                      title.size= "Nb.Org",
                                      border.col = "black", border.lwd = 0.1,
                                      border.alpha = 0.1) +
  tm_layout(legend.outside=T,legend.text.size = 0.8, legend.title.size = 1.1, frame = FALSE)

Static maps

  • mapsf

Note

mf_map(x=objet_sf,var =“variable”,type=“map type”,…)

WCA R-EiE Relief web page