Interagir avec PostGIS depuis R

PostGIS est une extentions de base de données PostgreSQL pour traiter les données spatiales. Grâce à PostGIS, votre base de données prend en charge les requêtes géographiques à exécuter directement en SQL. Dans ce post de blog, nous allons nous connecter et interagir avec une base de données PostGIS depuis R, en utilisant {DBI} et {sf}..

Le package {sf} et PostGIS sont potes

Le package {sf} est similaire à une base de données PostGIS de plusieurs points de vue :

• Le package R {sf} a été créé pour travailler au format “simple features” (d’où son nom, “sf”). “Simple features” est une norme ISO définie par l’Open Geospatial Consortium (OGC) pour normaliser le stockage et l’accès aux jeux de données géographiques. PostGIS gère les objets et fonctions de la norme “Simple Features for SQL” de l’OGC.

• La plupart des fonctions {sf} commencent par st_ pour suivre la dénomination des fonctions dans PostGIS. Par exemple, la fonction sf::st_union a son équivalent ST_Union dans Postgis. Si tu le savais pas, bah maintenant tu le sais !

• Les données spatiales avec {sf} sont manipulées à la manière SQL. Les fonctions simples de {sf} permettent de manipuler des objets géographiques vectoriels avec {dplyr}, dont la syntaxe est similaire aux noms de fonctions SQL. Par exemple, dplyr::select, qui peut être appliqué aux objets {sf}, a son équivalent SELECT dans la syntaxe SQL.

Configuration de votre session

Packages

Les principaux packages utilisés dans cet article sont {DBI}, {RPostgres}, {sqlpetr} et {sf}. Si vous voulez installer correctement les packages de manipulation de données géographiques sur un Ubuntu Server, vous pouvez vous relire notre post : installation de R sur ubuntu et astuces pour les packages de cartographie. Pour changer de notre bien-aimée palette Viridis, cet article contient quelques couleurs et palettes ThinkR personnalisées que vous devrez modifier pour reproduire le code. Notre package interne correspondant a été largement inspiré par l’article de @drsimonj : Création de palettes de couleurs pour entreprise avec ggplot2

library(DBI)
library(RPostgres)
library(sqlpetr)
library(sf)
# library(rpostgis)
library(dplyr)
library(dbplyr)
library(rnaturalearth)
library(ggplot2)
# ThinkR internal package for color palette and other brand identity
library(thinkridentity)

Vérification des drivers GDAL

Pour pouvoir connecter la base de données PostGIS avec R, vous aurez besoin d’une version GDAL qui supporte le pilote PostGIS. Vous pouvez le vérifier en utilisant sf::st_drivers().

# GDAL drivers
st_drivers() %>% 
  filter(grepl("Post", name))

##         name          long_name write  copy is_raster is_vector   vsi
## 1 PostgreSQL PostgreSQL/PostGIS  TRUE FALSE     FALSE      TRUE FALSE

Créer un conteneur Docker avec Postgres et PostGIS activés

Chez ThinkR, nous adorons utiliser les conteneurs Docker pour mettre R en production sur des serveurs. Nous avons déjà utilisé Docker pour déployer notre propre dépôt d’archives R, ou une application Shiny. Pour la reproductibilité de cet article de blog, il nous est donc tout naturel d’utiliser Docker. Ici, nous construisons un conteneur avec une base de données Postgres et PostGIS activé.

• Postgis sur DockerHub : https://hub.docker.com/r/mdillon/postgis
• Connexion : superutilisateur : “postgres”, pwd : “mysecretpassword”

Pour utiliser ce conteneur Docker en production sur un serveur, n’oubliez pas de définir un dossier persistant.

docker_cmd <- "run --detach --name some-postgis --publish 5432:5432 --env POSTGRES_PASSWORD=mysecretpassword -d mdillon/postgis"
system2("docker", docker_cmd, stdout = TRUE, stderr = TRUE)

## [1] "81371b924dd0156a2fb2d1b45ddc5d2e932b45c7ccc02926a20db29e4f694d46"

Connexion à la base de données PostGIS depuis R

Vous pouvez vous connecter à n’importe quelle base de données à l’aide du paquet {DBI}, à condition d’avoir les bons pilotes installés.

# With DBI
con <- dbConnect(
  RPostgres::Postgres(),
  host = "localhost",
  dbname = "postgres",
  port = 5432,
  user = "postgres",
  password = "mysecretpassword"
)

Écrire et lire les données mondiales dans PostGIS

Dans les exemples à suivre, nous allons jouer avec une carte du monde du package {rnaturalearth}. La projection “Equal Earth” est utilisée ici pour la représentation des cartes. La projection “Equal Earth” est dans la dernière version de “proj”. Si elle n’est pas disponible pour vous, vous pouvez en choisir une autre dans le code ci-dessous. Mais merci d’éviter la projection Mercator.

Nous pouvons écrire les données dans notre base de données PostGIS directement en utilisant sf::st_write. Les paramètres par défaut overwrite = FALSE, append = FALSE évitent d’écraser vos données si elles sont déjà dans la base de données. {sf} essaye lui-même de déterminer la méthode d’écriture dans la base.

st_write(ne_world, dsn = con, layer = "ne_world",
         overwrite = FALSE, append = FALSE)

Ensuite, nous pouvons lire ces données dans la base avec sf::st_read.

world_sf <- st_read(con, layer = "ne_world")

Utiliser les requêtes SQL avant de charger les données géographiques

Lorsque vous utilisez st_read, le jeu de données est chargé en totalité dans la mémoire de votre ordinateur. Si votre jeu de données spatiales est volumineux, vous voudrez peut-être n’en charger qu’une partie. Les requêtes SQL classiques peuvent être utilisées, par exemple pour extraire l’Afrique uniquement, en utilisant st_read avec le paramètre query.

query <- paste(
  'SELECT "name", "name_long", "geometry"',
  'FROM "ne_world"',
  'WHERE ("continent" = \'Africa\');'
)
africa_sf <- st_read(con, query = query)

Astuces : Avez-vous besoin d’aide pour votre requête SQL ?

PostGIS permet de faire des requêtes SQL avec des fonctions de géomatique que vous pouvez également définir avant de lire votre base de données avec {sf}. Cependant, si vous n’êtes pas totalement à l’aise avec SQL mais que vous êtes un ninja en {dplyr}, alors vous allez pouvoir vous en sortir. En utilisant dplyr::show_query, vous pouvez obtenir la syntaxe SQL de vos opérations {sf}. Prenons un exemple.
Avec {sf}, l’union des pays par continent pourrait s’écrire comme ça :

world_sf %>% 
 group_by(continent) %>% 
  summarise()

En réalité, la version complète du code {sf} serait plutôt :

world_sf %>% 
 group_by(continent) %>% 
  summarise(geometry = st_union(geometry))

Maintenant, connectons-nous à la base de données Postgres avec {dplyr}. Lorsqu’on est connecté à une base de données avec {dplyr}, les requêtes sont exécutées directement en SQL dans la base de données, et non dans la mémoire de l’ordinateur.
Utilisez ensuite show_query pour obtenir la traduction en SQL proposée par {dplyr}.

Notez qu’avec {dbplyr} (>1.4.0), au lieu de tbl, vous pouvez utiliser lazy_frame, ce qui évite une réelle connexion à votre base de données pour réaliser vos tests.

# Read with dplyr directly in database
world_tbl <- tbl(con, "ne_world")
# Create query
world_tbl %>% 
 group_by(continent) %>% 
  summarise(
   geometry = st_union(geometry)) %>% 
 show_query()

## <SQL>
## SELECT "continent", ST_UNION("geometry") AS "geometry"
## FROM "ne_world"
## GROUP BY "continent"

{dplyr} n’a pas été conçu pour traduire des opérations spatiales en SQL, mais comme vous pouvez le voir ci-dessous, la différence avec la requête correcte est très faible. Seule la fonction PostGIS n’est pas correctement écrite (ST_Union).

query <- paste(
  'SELECT "continent", ST_Union(geometry) AS "geometry"',
  'FROM "ne_world"',
  'GROUP BY "continent";'
)
union_sf <- st_read(con, query = query, quiet = TRUE)
# Plot (Choose your own palette)
union_sf %>% 
  st_transform(world_map_crs) %>%
  ggplot() +
  geom_sf(aes(fill = continent), colour = thinkr_cols("dark")) +
  scale_fill_thinkr_d(palette = "discrete2_long", guide = FALSE) +
  theme(panel.background = element_rect(fill = thinkr_cols("primary_light")))

Attention. Vous ne pouvez pas faire confiance aux opérations avec {dplyr}/{dbplyr} sur une base de données géographique car la colonne “geometry” n’est pas importée comme colonne spatiale avec tbl et vous la perdriez. Cependant, vous pouvez utiliser la syntaxe {dplyr} pour explorer votre jeu de données et toutes les colonnes, hors information spatiale. Toutes les opérations seront réalisées à l’intérieur de la base de données, ce qui permet de ne pas surcharger la mémoire.

Et les rasters ?

Si vous suivez ce blog, vous savez que j’aime les rasters. Le paquet {sf} ne traite que des objets spatiaux de type vectoriels. Il ne peut pas gérer les rasters dans une base de données PostGIS. Pour ça, vous pouvez utiliser le package {rpostgis} avec les fonctions pgGetRast() et pgWriteRast(). Ce package fonctionne avec les Rasters du package {raster}. Vous devrez vous connecter en utilisant le package {RPostgreSQL}, qui n’a pas de “contrat de connexion” Rstudio. Notez que vous pouvez avoir deux connexions à votre base de données, l’une avec {sqlpetr} pour conserver le volet “Connection” de Rstudio, et l’autre avec {RPostgreSQL} pour vous connecter avec les objets Raster.

• Connectez-vous à la base de données avec {RPostgreSQL}.

library(rpostgis)
con_raster <- RPostgreSQL::dbConnect("PostgreSQL",
  host = "localhost", dbname = "postgres",
  user = "postgres", password = "mysecretpassword", port = 5432
)
# check if the database has PostGIS
pgPostGIS(con_raster)

## [1] TRUE

• Créer un raster vide et le sauver dans la base PostGIS

# From documentation basic test
r <- raster::raster(
  nrows = 180, ncols = 360, xmn = -180, xmx = 180,
  ymn = -90, ymx = 90, vals = 1
)
# Write Raster in the database
pgWriteRast(con_raster,
  name = "test", raster = r, overwrite = TRUE
)

## [1] TRUE

• Lire le raster dans la base de données

r_read <- pgGetRast(con_raster, name = "test")
print(r_read)

## class       : RasterLayer 
## dimensions  : 180, 360, 64800  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 1, 1  (min, max)

• Ne pas oublier de se déconnecter de la base avant de partir

# Disconnect database
RPostgreSQL::dbDisconnect(con_raster)

## [1] TRUE

Se déconnecter de la base de données et fermer le conteneur Docker

• Déconnecter la connexion à la base de données faite avec {sqlpetr}
• Arrêter et supprimer le conteneur Docker si nécessaire

# Disconnect database
DBI::dbDisconnect(con)
# Stop Docker container
system2("docker", "stop some-postgis")
system2("docker", "rm --force some-postgis")
# Verify docker
system2("docker", "ps -a")

Pour aller plus loin

• ThinkR proposes son expertise pour l’installation d’infrastructures en lien avec R et les produits Rstudio

• ThinkR fait de la formation à la cartographie avec R. Cette formation inclut la présentation de PostGIS.

• Documentation de PostGIS

• [Dépôt de développement du package {sf}

• Introduction à la manipulation des objets spatiaux avec {sf} par Sébastien Rochette (moi)

• “Geocomputation in R” par Robin Lovelace et al.

• Introduction à Docker par Colin Fay

• “Introduction to R, databases and Docker” by John David Smith et al.

Packages utilisés dans cet article

# remotes::install_github("ThinkR-open/attachment")
sort(attachment::att_from_rmd("2019-04-27_postgis.Rmd"))

##  [1] "attachment"     "DBI"            "dbplyr"         "dplyr"         
##  [5] "ggplot2"        "knitr"          "raster"         "rnaturalearth" 
##  [9] "rpostgis"       "RPostgres"      "RPostgreSQL"    "sf"            
## [13] "sqlpetr"        "thinkridentity"