Cet atelier se veut un survol des outils disponibles en R pour l’analyse de données géoréférencées. Ce type de données apparaît de plus en plus fréquemment dans divers domaines (ex.: photos aériennes, images satellite, données du recensement, lieux rattachés aux messages sur les réseaux sociaux, etc.). Il existe deux grandes catégories de données géoréférencées:

  • les données matricielles représentent des variables définies à chaque point d’une grille couvrant tout l’espace représenté (comme une image satellite);
  • les données vectorielles associent des variables à des objets géométriques placés à des endroits précis (comme la position des villes et des chemins sur une carte routière).

Au départ, l’utilisation de commandes de programmation pour manipuler des données géographiques peut sembler moins intuitif que l’interface graphique de logiciels spécialisés (ex.: ArcGIS). Voici quelques avantages d’une analyse programmée:

  • Il est facile de répéter l’analyse pour de nouvelles données en ré-exécutant le programme.
  • Il est facile pour d’autres chercheurs de reproduire la méthodologie s’ils ont accès au même langage de programmation.
  • Dans le cas spécifique de R, on peut extraire des variables spatiales et les combiner à d’autres analyses statistiques avec un seul programme.

Objectifs

  • Se familiariser avec les principaux packages permettant le traitement et la visualisation simple de données vectorielles (sf) et matricielles (stars) en R.
  • Effectuer des transformations de données courantes à l’aide des fonctions de ces packages.
  • Connaître certains packages permettant des visualisations plus complexes: ggplot2 pour des cartes statiques et mapview pour des cartes interactives.

Note sur les packages

L’ensemble des packages disponibles pour l’analyse spatiale en R évolue rapidement. Il y a quelques années, les packages sp and raster étaient les principaux outils pour l’analyse des données vectorielles et matricielles, respectivement. sf et stars font partie d’une initiative récente pour faciliter le traitement des données spatiales dans R (https://www.r-spatial.org/). Le package sf représente les tableaux de données spatiaux selon un format standard basé sur les bases de données géospatiales et s’intègre bien avec des packages populaires pour la manipulation et la visualisation des données (notamment dplyr et ggplot2). Le package stars est plus récent, mais présente déjà certains avantages par rapport à raster. Une version précédente de cet ateiler utilisait raster plutôt que stars.

Explorer un jeu de données vectoriel

Tous les jeux de données de cet atelier se trouvent dans le répertoire data. Le jeu de données mrc contient les coordonnées des municipalités régionales de comté (MRC) du Québec en format ESRI shapefile. Notez que l’information pour chaque jeu de données est contenue dans plusieurs fichiers, qui portent le même nom mais diffèrent par leur extension (mrc.dbf, mrc.prj, mrc.shp and mrc.shx).

Pour charger un jeu de données dans R, nous appelons la fonction st_read (toutes les fonctions du package sf ont le préfixe st_, pour spatiotemporel). Le premier argument de st_read est le chemin vers le fichier de données. L’argument stringsAsFactors = FALSE évite que R ne convertisse les variables textuelles en facteur dans le tableau de données résultant.

library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

mrc <- st_read("data/mrc.shp", stringsAsFactors = FALSE)
## Reading layer `mrc' from data source `Z:\cours\rgeo-csee2019\data\mrc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 104 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -79.7625 ymin: 44.99136 xmax: -56.93495 ymax: 62.58217
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs

Le texte de sortie indique certaines propriétés du jeu de données chargé, incluant le type de géométrie (MULTIPOLYGON), les limites spatiales du jeu de données (bbox) et le système de coordonnées utilisé. Ce dernier est décrit sous deux formats: un code numérique epsg et une chaîne de caractère proj4string. Nous traiterons davantage des systèmes de coordonnées plus loin. Notons pour l’instant que “+proj=longlat” dans la proj4string signifie que les données sont en degrés de longitude et latitude.

Les limites spatiales et le système de coordonnées (CRS) peuvent être extraits séparément à l’aide des fonctions st_bbox et st_crs:

st_bbox(mrc)
##      xmin      ymin      xmax      ymax 
## -79.76250  44.99136 -56.93495  62.58217
st_crs(mrc)
## Coordinate Reference System:
##   EPSG: 4269 
##   proj4string: "+proj=longlat +datum=NAD83 +no_defs"

Regardons un aperçu du jeu de données:

class(mrc)
## [1] "sf"         "data.frame"
head(mrc)
## Simple feature collection with 6 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -79.51776 ymin: 45.44911 xmax: -63.31941 ymax: 62.58217
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
##                           mrc_name reg_id              reg_name pop2016
## 1                          Abitibi     08 Abitibi-Temiscamingue   24707
## 2                    Abitibi-Ouest     08 Abitibi-Temiscamingue   20582
## 3                            Acton     16            Monteregie   15623
## 4 Administration regionale Kativik     10        Nord-du-Quebec   13343
## 5                  Antoine-Labelle     15           Laurentides   35410
## 6                       Argenteuil     15           Laurentides   32477
##                         geometry
## 1 MULTIPOLYGON (((-78.64359 4...
## 2 MULTIPOLYGON (((-79.51776 4...
## 3 MULTIPOLYGON (((-72.70415 4...
## 4 MULTIPOLYGON (((-77.82293 5...
## 5 MULTIPOLYGON (((-75.52133 4...
## 6 MULTIPOLYGON (((-74.71004 4...

Un objet de classe sf est en fait un data.frame spécialisé, où chaque ligne contient des champs de données qui sont associés à un élément géométrique, décrit dans la colonne geometry. Les types géométriques les plus courants sont:

  • POINT: Coordonnées (x, y) d’un point.
  • LINESTRING: Séquence de points reliés par des segments de droite.
  • POLYGON: Séquence de points formant un polygone simple fermé.
  • MULTIPOINT, MULTILINESTRING ou MULTIPOLYGON: Jeu de données où chaque élément peut être composé de plusieurs points ou lignes ou polygones.

La fonction plot appliquée à un objet sf crée une carte de chaque attribut du jeu de données.

plot(mrc)

Pour n’afficher qu’un attribut, il faut sélectionner la colonne correspondante. Pour n’afficher que la carte sans données, on peut choisir la colonne geometry. Le paramètre axes = TRUE indique au programme d’afficher les axes de coordonnées.

plot(mrc[, "geometry"], axes = TRUE)

Vous pouvez choisir une partie des rangées ou des colonnes d’un objet sf comme s’il s’agissait d’un data.frame régulier.

# Choisir la 5e rangée
mrc[5, ]
## Simple feature collection with 1 feature and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -76.14932 ymin: 45.95376 xmax: -74.34359 ymax: 47.76333
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
##          mrc_name reg_id    reg_name pop2016
## 5 Antoine-Labelle     15 Laurentides   35410
##                         geometry
## 5 MULTIPOLYGON (((-75.52133 4...
# Choisir le nom et la population de la MRC pour les MRC avec une population de 200 000 et plus
mrc[mrc$pop2016 > 200000, c("mrc_name", "pop2016")] 
## Simple feature collection with 5 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -75.90834 ymin: 45.37075 xmax: -71.13162 ymax: 46.98186
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
##     mrc_name pop2016                       geometry
## 23  Gatineau  278198 MULTIPOLYGON (((-75.89532 4...
## 47     Laval  425225 MULTIPOLYGON (((-73.74488 4...
## 69 Longueuil  417497 MULTIPOLYGON (((-73.44474 4...
## 82  Montreal 1959836 MULTIPOLYGON (((-73.49801 4...
## 89    Quebec  572755 MULTIPOLYGON (((-71.55293 4...

Notez que la colonne contenant l’information spatiale (geometry) demeure toujours dans le tableau, même si elle n’est pas choisie. Pour éliminer cette colonne et convertir l’objet sf en tableau de données régulier, vous pouvez utiliser la fonction st_drop_geometry.

Exercice 1

Choisissez les MRC de la région du Bas-St-Laurent (reg_id: 01) et de la Gaspésie (reg_id: 11), puis affichez leur population de 2016 sur une carte. Indice: L’opérateur %in% permet de vérifier si une variable prend une valeur parmi une liste. Par exemple, x %in% c(1, 3) est égal à TRUE si x = 1 ou x = 3.

Solution

Intégration avec le package dplyr

Les fonctions de manipulation de données du package dplyr fonctionnent aussi sur des objets sf. Par exemple, nous pouvons réécrire l’exemple ci-dessus (choisir le nom et la population des MRC avec une population de plus de 200 000) avec filter et select.

library(dplyr)
## 
## 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

filter(mrc, pop2016 > 200000) %>%
    select(mrc_name, pop2016)
## Simple feature collection with 5 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -75.90834 ymin: 45.37075 xmax: -71.13162 ymax: 46.98186
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
##    mrc_name pop2016                       geometry
## 1  Gatineau  278198 MULTIPOLYGON (((-75.89532 4...
## 2     Laval  425225 MULTIPOLYGON (((-73.74488 4...
## 3 Longueuil  417497 MULTIPOLYGON (((-73.44474 4...
## 4  Montreal 1959836 MULTIPOLYGON (((-73.49801 4...
## 5    Quebec  572755 MULTIPOLYGON (((-71.55293 4...

Lorsque nous calculons une statistique sommaire par groupe, les éléments géométriques sont aussi soudés pour chaque groupe. Pour illustrer ce point, regroupons les MRC et leur population par région.

regions <- group_by(mrc, reg_name) %>%
    summarize(pop2016 = sum(pop2016))
head(regions)
## Simple feature collection with 6 features and 2 fields
## geometry type:  GEOMETRY
## dimension:      XY
## bbox:           xmin: -79.57933 ymin: 45.58943 xmax: -56.93495 ymax: 55.00006
## epsg (SRID):    4269
## proj4string:    +proj=longlat +datum=NAD83 +no_defs
## # A tibble: 6 x 3
##   reg_name        pop2016                                          geometry
##   <chr>             <dbl>                                   <GEOMETRY [°]>
## 1 Abitibi-Temisc~  147282 POLYGON ((-77.93147 47.27012, -77.93155 47.26624~
## 2 Bas-Saint-Laur~  197806 POLYGON ((-68.38319 47.88014, -68.38314 47.83307~
## 3 Capitale-Natio~  733898 POLYGON ((-70.69151 47.03311, -70.69453 47.02324~
## 4 Centre-du-Queb~  243354 POLYGON ((-72.09283 45.79699, -72.11938 45.77338~
## 5 Chaudiere-Appa~  421993 POLYGON ((-70.2784 46.05675, -70.2784 46.05608, ~
## 6 Cote-Nord         92712 MULTIPOLYGON (((-66.25978 55.00001, -66.25001 55~
plot(regions["pop2016"])