FAIRE DES CARTES AVEC R

# FAIRE DES CARTES AVEC R
## La frontière États-Unis - Mexique
### Nicolas Lambert & Ronan Ysebaert
### 20 novembre 2019

---

---
# Import des géométries utiles
Les géométries sont importées avec les fonctions du package **`rnaturalearth`** (et diverses sources).

```r
countries <- ne_countries(scale = 50, type = "countries", continent = NULL,
                          country = NULL, geounit = NULL, sovereignty = NULL,
                          returnclass = "sf")
```
--
```r
rivers <- ne_download(scale = 50, type = "rivers_lake_centerlines", 
                      category = "physical", returnclass = "sf")
```
--
```r
coastline <- ne_download(scale = 50, type = "coastline",
                         category = "physical", returnclass = "sf")
```
--
```r
ocean <- ne_download(scale = 50, type = "ocean", category = "physical",
                     returnclass = "sf")
```
--
```r
subregions <- st_read(dsn = "data/regions/mex_us_admin.shp" ,
                      options = "ENCODING=UTF-8", stringsAsFactors = FALSE)
```
--
```r
fences <- geojson_sf("data/data.world/border-fence.geojson")
```

---
# Quelques transformations préalables
Puis sont transformées dans la projection cible du modèle cartographique et intersectées avec l'emprise de l'espace d'étude (fonctions du package **`sf`**).

```r
# Choix de la projection
*albers <- "+proj=aea +lat_1=14.5 +lat_2=32.5 +lat_0=24 +lon_0=-105 +x_0=0 +y_0=0
*          +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
```

```r
# Transformation des couches dans la bonne projection
countries_aea <- st_transform(countries,crs = albers) 
```
--
```r
# Création de la bounding box
bbox_aea <- st_as_sfc(st_bbox(c(xmin = -1542784 , xmax = 993341, 
                                ymin = 1417850, ymax = -839750), 
                              crs = albers))
```
--
```r
# Intersection du fond régional avec la bbox
subregions_aea <- st_intersection(x = subregions_aea, st_geometry(bbox_aea))
```

---

# Template 1 (projection Albers)

La fonction *`lay_aea`* crée le modèle cartographique qui sera utilisé dans les représentations ultérieures : ordre d'affichage et styles associés des couches géographiques.

```r

lay_aea <- function(title = ""){
  par(mar = c(0,0,1.2,0))
  plot(st_geometry(ocean_aea), col= "#b8d5e3", border = NA, xlim = bb_aea[c(1,3)],
       ylim = bb_aea[c(2,4)])
  plot(st_geometry(subregions_aea) + c(-10000, -10000), col ="#827e6c50",
       border = NA, add = T)
  plot(st_geometry(subregions_aea), col= "#ede6bb", border = "white",
       cex = 0.5, add=T)
  plot(st_geometry(coastline_aea), col= "#6d9cb3",lwd = 1 ,add= T)
  plot(st_geometry(rivers_aea), col= "#6d9cb3",lwd = 1 ,add= T)
  plot(st_geometry(fences_aea), col= "#3d3c3c",lwd = 3 ,add= T)
  layoutLayer(title = title,
              author =  authors,
              scale = 300, south = TRUE, frame = TRUE,
              col = "#6d9cb3", coltitle = "white")
}

```
--
```r
sizes_aea <- getFigDim(x = bbox_aea, width = 1500,mar = c(0,0,1.2,0), res = 150)
```
--
```r
png("img/fig01.png", width = sizes_aea[1], height = sizes_aea[2], res = 150)
*lay_aea("Template cartographique 1 (projection Albers)")
dev.off()
```

---
background-image: url(img/fig01.png)
background-size: contain

---

# Template 2 (projection orthographique)

Un second modèle "vu du sud" est créé en utilisant une projection orthographique.

```r

* ortho <- "+proj=ortho +lat_0=-35 +lon_0=-104 +x_0=0 +y_0=0 +ellps=
*           WGS84 +units=m +no_defs"

```

![](img/proj-persp5x.png)

--
```r
*lay_ortho("Template cartographique 2 (projection orthographique)")
```
---

---

# Template 2bis (projection orthographique et mur effet 2.5D)

En translatant 20 fois de 5000 m vers le nord le mur frontalier, on renforce
graphiquement l'idée du mur.

```r
   line <- st_geometry(fences_ortho)
   plot(line, col= "#3d3c3c",lwd = 2 ,add= T)

*   for (i in 1:20){
*    line <- line + c(0,5000)
*     plot(line, col= "#565b6380",lwd = 2 ,add= T)  
*   }

plot(line, col= "#3d3c3c",lwd = 2 ,add= T)

```

---
background-image: url(img/fig03.png)
background-size: contain

---

# Évolution du PIB et des moins de 20 ans dans le temps

Plusieurs packages utilisent les API des principaux fournisseurs
de données internationaux, comme **`OECD`**.

C'est un moyen efficace sans sortir de R pour récupérer des données.

```r

# Télécharger les données de la table PDB_LV pour USA, Mexique, pays de l'OCDE
* df <- get_dataset(dataset = "PDB_LV",
*                   filter = list(c("MEX", "USA","OECD"), "T_GDPPOP", "CPC"))

# Transformer la date au format numérique
df$obsTime <- as.numeric(df$obsTime)

```

Le package **`ggthemes`** propose des modèles de mise en forme assez esthétiques.

```r
# Représentation graphique
ggplot(data = df, aes(x = obsTime, y = obsValue, color = LOCATION)) + 
  geom_line(size = 1) +  
  labs(x = NULL, y = "Dollars, prix courant", color = NULL,
       title =  "Évolution comparée du PIB par habitant (Mexique - USA - OCDE)") +
*  theme_hc() +
  scale_color_manual(values=c("#999999", "#E69F00", "#56B4E9"))

```

---
background-image: url(img/fig05.png)
background-size: contain

---
background-image: url(img/fig06.png)
background-size: contain

---

# Des murs dans les données ?

Les fonctions *`getBorders`* et  *`discLayer`* du  package **`cartography`**  permettent d'extraire les frontières entre des polygones et représenter l'importance des discontinuités territoriales.

```r
# Carte choroplèthe
choroLayer(x = subregions_ortho, var = "POP65_POP15", 
           breaks = c(min(subregions_ortho$POP65_POP15, na.rm = T),
                      20,25,35,50,65, max(subregions_ortho$POP65_POP15, na.rm = T)),
           col = carto.pal(pal1 = "green.pal", n1 = 3, pal2 = "red.pal", n2 = 3),
           legend.pos = c(-1700000, 5000000),
           legend.horiz = TRUE, legend.title.cex = 0.7, legend.values.cex = 0.5,
           legend.title.txt = "Rapport entre la population âgée de plus de
                               65 ans\net la population âgée de moins de 15 ans\nen 2015 (%)",
           border = NA, add = TRUE)
```
--
```r
# Extraction des frontières
* subregions.borders <- getBorders(subregions_ortho)
```
--
```r
# Discontiuités
*discontinuities <- discLayer(x = subregions.borders, df = subregions_ortho,
                             var = "POP65_POP15", col = "black", nclass=3,
                             method = "equal", threshold = 0.3, sizemin = 0.5,
                             sizemax = 10, type = "abs",legend.values.rnd = 0,
                             legend.title.txt = "Discontinuités sur l'indice\nde veillissement 2015\n(différences\nabsolues)",
                             legend.pos = c(-1700000, 5300000), legend.title.cex = 0.7, 
                             legend.values.cex = 0.5, add = TRUE)
```

---
background-image: url(img/fig07.png)
background-size: contain

---
background-image: url(img/fig08.png)
background-size: contain

---

# Des (vrais) murs dans les données ?

Ce bout de code permet de translater les frontières proportionnellement à la valeur de leur discontinuité respective. Plus le mur est haut, plus la discontinuité est importante !

```r
# On ne garde que les discontinuités les plus fortes (30%) 
threshold <- 0.3
minvar <- as.numeric(quantile(discontinuities$disc, probs = c(1 - threshold)))
discontinuities <- discontinuities[discontinuities$disc >= minvar,]
```
--
```r
# On déduit de la valeur des discontinuités un nombre d'itérations
discontinuities$height <- round(discontinuities$disc / 2,0)
```
--
```r
# On crée une fonction pour itèrer chaque ligne en augmentant la valeur en Y. 
extrude <- function(id){
 line <- st_geometry(discontinuities[id,])
 plot(line, col = "black",lwd = 2 , add = TRUE)
* nb <- as.numeric(discontinuities[id,"height"])[1]
* for (j in 1:nb){
*   line <- st_geometry(line) + c(0,5000)
*   plot(st_geometry(line), col= "#ebd23490", lwd = 2 ,add = TRUE)  
 }
 plot(line, col= "black", lwd = 2 , add = TRUE)
}

```
--
```r
# On dessine
for (i in 1:length(discontinuities$disc))
{
extrude(i)
}

```

---
background-image: url(img/fig09.png)
background-size: contain

---
background-image: url(img/fig10.png)
background-size: contain

---
background-image: url(img/fig11.png)
background-size: contain

---
background-image: url(img/fig12.png)
background-size: contain

---

# Des points de contrôle...

Le package **`osmdata`** permet d'effectuer des requêtes (pas trop importantes)
sur le contenu de la base OpenStreetMap via l'API Overpass Turbo.

La fonction *`osmdata_sf`* retourne un objet sf en fonction d'une bounding box prédéfinie et d'un couple **clé/valeur**.

```r

# Convertir la bounding box en WGS 84
bbox <- st_transform(bbox_aea, 4326)

# Définir la requête (clé/valeur OSM sur bounding box)
*opqbox <- opq(bbox = bbox , timeout = 5000)
*opquery <- add_osm_feature(opq = opqbox, key = "barrier", value = "border_control")
*feat <- osmdata_sf(opquery)

```
---

# Des points de contrôle...

On extrait les centroides des polygones que l'on combine aux points qui répondent à la requête...

```r

# Extraire les points qui répondent à la requête
featpt <- st_transform(feat$osm_points, albers)
featpt <- featpt[featpt[["barrier"]] %in% "border_control", ]

# Extraire les polygones qui répondent à la requête
featpo <- st_transform(feat$osm_polygons, albers)
st_geometry(featpo) <- st_centroid(st_geometry(featpo))
featpo$osm_id <- row.names(featpo)

# Combiner points et polygones, les intersecter avec la bounding box
featpt <- rbind(featpt[, c("osm_id", "geometry")], featpo[, c("osm_id", "geometry")])
poi_osm <- st_intersection(x = featpt, st_geometry(subregions_aea))

```

... Puis on les comptabilise et on les représente.

```r

*lay_aea(paste0("Au moins ", sum(grid$ncops)," postes frontières..."))
*plot(st_geometry(poi_osm), bg = "red", col = NA, pch = 21, cex = 0.8, add = TRUE)
legtxt <- "Chaque point rouge\nreprésente un poste\nfrontalier recensé\ndans OpenStreetMap."
text(-1400000, y = -100000, legtxt  , cex = 0.9, pos = 4, font = 2)

```

---
background-image: url(img/fig13.png)
background-size: contain

---

# Aux carreaux de contrôle...

Pour y voir plus clair, on crée une grille de 50 km² et on dénombre de postes frontières pour chaque élément de la grille.
Les postes frontières peuvent alors être représentés sous forme de figurés proportionnels...

```r

# Créer une grille sur l'espace d'étude
grid <- st_make_grid(subregions_aea, cellsize = 50000)
grid <- st_sf(grid)

# Compter le nombre de postes de police par points de grille
*grid$ncops <- lengths(st_covers(grid, poi_osm))

# Figurés proportionnels
lay_aea("Les hot-spots du contrôle frontalier - dénombrement")
propSymbolsLayer(grid, var = "ncops", col = "red", symbols = "square", add = T,
                 legend.pos = "left",
                 legend.title.cex = 0.7, legend.values.cex = 0.6,
                 legend.title.txt = "Nombre de postes frontière\n(zones de 50km²)")

```

---
background-image: url(img/fig14.png)
background-size: contain

---

# Aux carreaux de contrôle...

Ou directement sur le carroyage. On s'intéresse ici à la densité des points de contrôle,
exprimés en part du total de l'espace d'étude.

```r

# Calcul de densité
grid$dens <- (grid$ncops / sum(grid$ncops)) * 100
grid <- grid[grid$ncops != 0, ]

lay_aea("Les hot-spots du contrôle frontalier - Part du total")
choroLayer(x = grid, var = "dens", breaks = c(min(grid$dens), 1, 2, 5, 10, max(grid$dens)),
           col = carto.pal(pal1 = "brown.pal", n1 = 5),
           legend.pos = "left",
           legend.title.cex = 0.7, legend.values.cex = 0.6,
           legend.title.txt = "Part des postes frontaliers sur l'espace d'étude (%)",
           border = NA,
           add = TRUE)
plot(st_geometry(fences_aea), col= "#3d3c3c",lwd = 3 ,add= T)

```

---
background-image: url(img/fig15.png)
background-size: contain

---

# Et aux tours de contrôle !

L'usage de l'argument **bar** de la fonction *`propSymbolsLayer`* (package **`cartography`**), associé à une projection orthographique permet de donner une impression 3d intéressante.

L'ajout du mur de séparation permet d'apprécier que ces "tours de contrôle" sont bien situées en territoire américain...

```r

lay_ortho(title = "Sacrées tours de contrôle !")

propSymbolsLayer(x = grid, var = "ncops", col = "darkred",
*                symbols = "bar",
                 inches = 1.3,
                 border = "white", lwd = 1, legend.pos = c(-1700000, 5000000), 
                 legend.title.txt = "Nombre de postes frontières\n(zones de 20 km²)",
                 legend.style = "e")

line <- st_geometry(fences_ortho)
for (i in 1:15){
  line <- st_geometry(line) + c(0,5000)
  plot(st_geometry(line), col= "#565b6380",lwd = 2 ,add= T)  
}
plot(st_geometry(line), col= "#3d3c3c",lwd = 2 ,add= T)

```

---

---

# Morts aux frontières - en chiffres

L'OIM compile la localisation de personnes qui ont trouvé la mort au cours de leur trajectoire migratoire.

![](img/iom.png)

---

# Morts aux frontières - en chiffres

Après avoir importé les données et les avoir nettoyé, on peut réaliser plusieurs barplots pour
visualiser l'importance de la tragédie migratoire dans l'espace...

```r
# Import du fichier brut (OIM)
*iom <- read.csv("data/iom/MissingMigrants-Global-2019-10-29T14-11-50.csv", stringsAsFactors = F)
```
--
```r
# Un peu de nettoyage

iom <- iom[(iom$Location.Coordinates)!="",]
latlon <- matrix(as.numeric(unlist(strsplit(iom$Location.Coordinates, split = ", "))), 
                 ncol = 2, byrow = T)
colnames(latlon) <- c("lat", 'lon')
iom <- cbind(iom, latlon)

iom <- iom[,c("Web.ID","Reported.Year","Total.Dead.and.Missing","Number.of.Survivors",
              "Region.of.Incident","lat","lon")]
colnames(iom) <- c("id","year","deads","survivors","region","latitude","longitude")
iom$deads <- as.numeric(iom$deads)
iom <- iom[!is.na(iom$deads),]
iom$latitude <- as.numeric(iom$latitude)
iom$longitude <- as.numeric(iom$longitude)

```

---

# Morts aux frontières - en chiffres

```r
# Agréger par zone géographique
*med <- aggregate(iom_sf$deads,list(iom_sf$region), sum, simplify = TRUE )
barplot(med$nb, ylab = "Nombre de personnes", names.arg = med$region, las = 2, 
        border="#991313",col = cols, cex.names = 0.8, cex.axis = 0.8)

```
--
<img src="img/fig17.png" width="90%"></img>

---

# Morts aux frontières - en chiffres
```r
# Agréger par années
*med <- aggregate(iom_sf$deads,list(iom_sf$year), sum, simplify = TRUE )
barplot(med$nb, xlab=paste0("Total sur la période: ",total,"\n(*) Du 1er janvier au 29 octobre 2019"), 
        ylab="Nombre de personnes", names.arg=med$year,
        border="#991313",col=c("red","red","red","red","red","#ffbaba"),
        cex.names = 0.8, cex.axis = 0.8)
        
```
--
<img src="img/fig18.png" width="80%"></img>

---

# Morts aux frontières - carte de localisation

On transforme le dataframe en objet sf pour obtenir la localisation et on l'intersecte avec l'emprise de l'espace d'étude.

```r
iom_sf <- st_as_sf(iom, coords = c("longitude", "latitude"), crs = 4326, agr = "constant")
iom_ortho <- st_transform(iom_sf,crs = ortho)
iom_ortho <- st_intersection(x = iom_ortho, st_geometry(bbox_ortho))
```

On crée la carte

--
```r
# Affichage du template
lay_ortho("Migrants morts et portés disparus à la frontière USA-Mexique, 2014 - 2019")
```
--
```r
# Affichage des points
*plot(st_geometry(iom_ortho), pch=20, col= "#eb3850", cex = 0.5, add=T)
```
--
```r
# Légende
legtxt <- "Sur cette carte, chaque point\ncorrespond
           à un évenement\nayant donné la mort d'au\nmoins une
           personne\nsur la période 2014 - 2019"
text(-1700000, y = 5200000, legtxt  , cex = 0.9, pos = 4, font = 2) 
```

---
background-image: url(img/fig19.png)
background-size: contain

---

# Morts aux frontières - Localisation et nombre de personnes

Dans la base de données de l'OIM, un événement peut regrouper plusieurs décès...

```r

lay_ortho("Migrants morts et portés disparus à la frontière USA-Mexique, 2014 - 2019")

*propSymbolsLayer(x = iom_ortho, var = "deads",
*                symbols = "circle", col =  "#eb3850",
*                legend.pos = "left", border = "#ede6bb", lwd = 0.5,
*                legend.title.txt = "Nombre de morts\net portés
*                                    disparus\nsur la période\n2014 - 2019",
*                legend.style = "e")

plot(st_geometry(fences_ortho), col= "#3d3c3c",lwd = 3 ,add= T)

```

---

---

# Morts aux frontières - Cartogramme de Dorling

La fonction *`cartogram_dorling`* du package **`cartogram`** permet de d'éviter le recouvrement graphique des cercles.
La fonction *`st_jitter`* du package **`sf`** optimise leur localisation.

```r

# Localisation des points précise plus importante pour gros événements
iom_ortho$m_weight <- 1
iom_ortho$m_weight[iom_ortho$deads > 1] <- 0.5
iom_ortho$m_weight[iom_ortho$deads >= 25] <- 0
```
--

```r
# Cartogramme Dorling / Jitter (pour ajouter un peu de bruit dans la localisation)
*deathsdor <- cartogram_dorling(x = st_jitter(iom_ortho),
                               weight = "deads", 
                               m_weight = iom_ortho$m_weight, k = .4)
```
--
```r
# Réalisation de la carte
lay_ortho("Migrants morts et portés disparus à la frontière USA-Mexique, 2014 - 2019")
plot(st_geometry(deathsdor), pch=20, col= "#eb3850", border ="#ede6bb", cex = 0.1, add=T)
plot(st_geometry(fences_ortho), col= "#3d3c3c", lwd = 3 ,add= T)

```

---

---

# Morts aux frontières - Cartogramme de Dorling (avec désagrégation)

Cette boucle permet de désagréger les événements comprenant plusieurs décès de telle
sorte qu'un point représente un décès.

```r
iom_unique <- all[all$deads == 1,]
iom_multi <-  all[all$deads > 1,]

*for (i in 1:dim(iom_multi)[1]){
*  nb <- as.numeric(iom_multi[i,"deads"])[1]
*  tmp <- iom_multi[i,]
*  tmp$deads <- 1
*  for (j in 1:nb){ iom_unique <- rbind(iom_unique,tmp)}
*}

deathsdor2 <- cartogram_dorling(x = st_jitter(iom_unique), weight = "deads", k = .003)

```

---
background-image: url(img/fig22.png)
background-size: contain

---

# Morts aux frontières - Carte animée

Création d'un vecteur avec toutes les dates

```r

m <- c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")
y <- c("2014","2015","2016","2017","2018","2019")

for (i in 1:length(y)){
  d <- paste0(m," ",y[i])
  if (i == 1) { all <- d } else {all <- c(all,d)}
}

```
--
Nombre de morts par mois
```r
bymonth <- aggregate(iom_ortho$deads,list(iom_ortho$date), sum, simplify = TRUE )
all <- merge (x = all, y = bymonth, 
              by.x = "date",
              by.y = "date",
              all.x = T)
```
--
Création d'un objet sf
```r
iom_ortho <- merge (x = iom_ortho, y = all[,c("date","id")], 
                    by.x = "date",
                    by.y = "date",
                    all.x = T)
```

---
# Morts aux frontières - Carte animée

Nombre de morts au cours du mois (en jaune)
```r
*layer1 <- iom_ortho[iom_ortho$id == i,]
propSymbolsLayer(x = layer1, var = "deads",
                 symbols = "circle", col =  "#eb385040",
                 inches = inches, fixmax = fixmax,
                 legend.pos = NA, border = NULL, lwd = 0.5)
```
--
Nombre de morts depuis le 1er janvier 2014 (en rouge)
```r
*layer2 <- iom_ortho[iom_ortho$id <= i,]
propSymbolsLayer(x = layer2, var = "deads",
                 symbols = "circle", col =  "#ffe100",
                 inches = inches, fixmax = fixmax,
                 legend.pos = NA, border = "black", lwd = 1)
  
```
--
On crée une boucle pour créer une carte par mois dans le template orthographique (largement personnalisé graphiquement)
```r
for (i in 1:length(all$id)){
...
}
```
---
background-image: url(img/animate.gif)
background-size: contain

---

# Morts aux frontières - Carte interactive

Avec le package **`leaflet`**.

```r
pins <- makeIcon(
  iconUrl = "data/pin.svg",
  iconWidth = 30, iconHeight = 30,
  iconAnchorX = 15, iconAnchorY = 15
)
```
![](data/pin.svg)
--

Création des infobulles (en html)
```r
iom$label <- paste0(
  "<h1>",iom$cause,"</h1><h3>year: </b>",iom$year,"
   <br/>location: ",iom$location,"</h3><i>Source: ",iom$source,"</i>"
)
```

--
Création de la carte leaflet (cluster)
```r
m <- leaflet(iom) %>%
  addProviderTiles(providers$Esri.WorldStreetMap) %>%
  setView(lng = -104, lat = 30, zoom = 06) %>%
  addMarkers(~lon, ~lat, popup = ~label, clusterOptions = markerClusterOptions(), icon = pins ) %>%
  addScaleBar(position = "bottomleft") %>%
  addPolylines(data = fences, color = "black", weight = 7, opacity = 1)
```

---

---
class: clear, bottom

`NICOLAS LAMBERT`<br/>
[**nicolas.lambert@cnrs.fr**](mailto:nicolas.lambert@cnrs.fr)

`RONAN YSEBAERT`<br/>
[**ronan.ysebaert@cnrs.fr**](mailto:nronan.ysebaert@cnrs.fr)

`CODE SOURCE`<br/>
[**github.com/riatelab/mexusaborder**](https://github.com/riatelab/mexusaborder)

`PRESENTATION`<br/>
[**riatelab.github.io/mexusaborder**](https://riatelab.github.io/mexusaborder)