cartographie, projections, sf, rnaturalearth, ggplot2, R
Recentrer, zoomer, éplucher : trois projections, trois bonnes raisons de tout envoyer valser.
Cet article s’appuie sur mes travaux publiés dans Geophysical Research Letters :
Panaïotis, T., Wilson, J., & Cael, B. B. (2025). A machine learning-based dissolved organic carbon climatology. Geophysical Research Letters, 52(7), e2024GL112792. https://doi.org/10.1029/2024GL112792.
Cet article est né d’un cas bien réel. Pour un papier publié dans Geophysical Research Letters sur la distribution mondiale du carbone organique dissous (le DOC, pour dissolved organic carbon) dans l’océan, mon encadrant m’avait demandé des cartes en projection Robinson centrées sur le Pacifique plutôt que sur l’Atlantique, une habitude courante en océanographie pour ne pas couper le Pacifique en deux au milieu de la carte. Sur le papier, ça aurait dû me prendre 10 minutes : changer la projection et son centre semblait anodin, un simple argument en plus dans la chaîne PROJ. Dans la pratique, ça a viré au cauchemar : les continents se chevauchaient et les données de DOC traçaient des lignes parasites en plein milieu de l’océan. Revenons sur une journée entière de prise de tête.
Reproduisons ce cas avec les vraies données du papier, publiques sur le dépôt GitHub du projet. Commençons par préparer de quoi tracer nos cartes.
library(tidyverse)
library(sf)
library(rnaturalearth)
world <- ne_countries(scale = "medium", returnclass = "sf")
doc <- read_csv("https://raw.githubusercontent.com/ThelmaPana/DOClimato/main/output/annual_climatologies.csv") |>
select(lon, lat, doc_avg = surf_doc_avg)
head(doc)# A tibble: 6 × 3
lon lat doc_avg
<dbl> <dbl> <dbl>
1 -180. -76.5 50.0
2 -180. -75.5 47.9
3 -180. -74.5 47.8
4 -180. -73.5 45.7
5 -180. -72.5 46.1
6 -180. -71.5 46.1doc contient une grille régulière au degré près de carbone organique dissous prédit en surface de l’océan, en µmol par kg. Rien d’exotique pour l’instant, une carte en projection lat/lon classique s’en sort très bien.
sf ou pas ?
doc est un tableau classique, pas un objet sf : chaque ligne est un simple triplet lon/lat/doc_avg, sans colonne geometry. geom_tile() s’en accommode très bien, il trace un rectangle par combinaison unique de x et y, sans rien savoir des projections. C’est justement le rôle de default_crs dans coord_sf() : lui dire que ces nombres bruts sont des coordonnées lon/lat, à reprojeter comme le reste.
C’est justement cette demande de recentrage sur le Pacifique qui m’a fait sortir de coord_quickmap() et coord_map(), que j’utilisais encore à l’époque du papier bien qu’elles soient déjà dépréciées. Pour cet article de blog, j’ai donc été tentée de tout convertir doc en sf (avec st_as_sf(), en reconstruisant un polygone par case à partir de lon/lat et de la résolution de la grille) pour rester cohérente avec world. J’ai finalement gardé la version non sf avec geom_tile(), plus proche de mon code d’origine.
ggplot(doc) +
geom_tile(aes(x = lon, y = lat, fill = doc_avg, colour = doc_avg)) +
geom_sf(data = world, fill = "gray80", colour = NA) +
scale_fill_viridis_c(option = "F", trans = "log1p") +
scale_colour_viridis_c(option = "F", trans = "log1p", guide = "none") +
coord_sf(crs = "EPSG:4326") + # système de coordonnées lon/lat
labs(fill = "DOC (μmol kg<sup>-1</sup>)") +
theme_void() +
theme(legend.title = element_markdown())
Une palette adaptée
viridis, ici sa variante « rocket » (option = "F"), est une palette séquentielle : pensée pour une variable continue sans centre significatif, comme une concentration. Elle est perceptuellement uniforme et lisible pour les daltoniens, contrairement à une palette arc-en-ciel. Le trans = "log1p" étire l’échelle sur les valeurs basses, largement majoritaires, sans quoi les quelques zones à forte concentration écraseraient toutes les nuances. Plus de détails dans l’article sur les palettes de couleurs.
Passons d’abord en Robinson tout court, celle de l’article précédent, sans toucher à son centre par défaut.
robin <- "+proj=robin"
ggplot(doc) +
geom_tile(aes(x = lon, y = lat, fill = doc_avg, colour = doc_avg)) +
geom_sf(data = world, fill = "gray80", colour = NA) +
scale_fill_viridis_c(option = "F", trans = "log1p") +
scale_colour_viridis_c(option = "F", trans = "log1p", guide = "none") +
coord_sf(crs = robin, default_crs = sf::st_crs(4326)) +
labs(fill = "DOC (μmol kg<sup>-1</sup>)") +
theme_void() +
theme(legend.title = element_markdown())
crs et default_crs, à quoi ça sert ?
Dans coord_sf(), crs fixe la projection cible, celle de tous les objets sf déjà tracés (ici world). default_crs, lui, s’applique uniquement aux données qui n’ont pas de système de coordonnées propre, comme les colonnes lon et lat brutes de doc utilisées par geom_tile(). Sans lui, ggplot ne saurait pas dans quel système ces nombres sont exprimés, ni comment les reprojeter pour qu’ils tombent au bon endroit sur la carte.
Aucune surprise : la couture1 de la carte reste là où elle a toujours été, à 180° de longitude, exactement là où world et la grille de DOC sont déjà découpés. Rien ne casse.
Recentrons maintenant sur le Pacifique : pour cela on ajoute lon_0=180 à notre chaîne PROJ, puis on regarde nos espoirs s’effondrer.
robin_pacific <- "+proj=robin +lon_0=180"
ggplot(doc) +
geom_tile(aes(x = lon, y = lat, fill = doc_avg, colour = doc_avg)) +
geom_sf(data = world, fill = "gray80", colour = NA) +
scale_fill_viridis_c(option = "F", trans = "log1p") +
scale_colour_viridis_c(option = "F", trans = "log1p", guide = "none") +
coord_sf(crs = robin_pacific, default_crs = sf::st_crs(4326)) +
labs(fill = "DOC (μmol kg<sup>-1</sup>)") +
theme_void() +
theme(legend.title = element_markdown())
Deux problèmes distincts, mais une seule et même cause. Passer le centre de 0 à 180 revient aussi à faire passer la couture de 180 à 0 : les deux ont simplement échangé leur place. Problème, ni les pays ni la grille de DOC n’ont été pensés pour être coupés à cette nouvelle longitude. Toute géométrie qui la chevauche encore se retrouve avec une partie de ses points tout à gauche de la carte reprojetée et l’autre tout à droite, alors qu’ils devraient rester côte à côte :
lon = -0.5 (qui couvrent donc de -1 à 0), une par latitude océanique : leurs deux bords se retrouvent eux aussi aux extrémités opposées de la carte, et geom_tile() relie ces points en autant de bandes qui traversent tout l’océan.Une donnée gridée n’est pas plus à l’abri qu’un polygone du problème de couture : toute case, tout pixel qui chevauche la nouvelle limite doit être découpé ou filtré avant la reprojection, jamais laissé à coord_sf().
On coupe les données à l’endroit exact de la nouvelle couture, avant de les reprojeter, exactement comme pour le fond de carte de l’article précédent. Ici, la couture tombe à la longitude 0 : on découpe world avec une lame trèèèèèès fine à cette longitude.
# Lame fine à la longitude 0, la nouvelle couture une fois recentré sur le Pacifique
seam <- st_polygon(list(rbind(
c(-0.0001, 90), c(0, 90), c(0, -90), c(-0.0001, -90), c(-0.0001, 90)
))) |>
st_sfc(crs = 4326)
world_robin <- world |>
st_difference(seam) |>
st_transform(crs = robin_pacific)world_robin est désormais déjà dans le bon système de coordonnées : plus besoin de préciser crs = dans coord_sf(), seulement default_crs pour que la grille de DOC, elle, toujours en lat/lon brut, soit reprojetée au même endroit.
ggplot(doc) +
geom_tile(aes(x = lon, y = lat, fill = doc_avg, colour = doc_avg)) +
geom_sf(data = world_robin, fill = "gray80", colour = NA) +
scale_fill_viridis_c(option = "F", trans = "log1p") +
scale_colour_viridis_c(option = "F", trans = "log1p", guide = "none") +
coord_sf(default_crs = sf::st_crs(4326)) +
labs(fill = "DOC (μmol kg<sup>-1</sup>)") +
theme_void() +
theme(legend.title = element_markdown())
Le problème est résolu pour les continents ! En revanche, pour le DOC océanique, ce n’est pas encore tout à fait cela.
Toutes les cases à lon = -0.5 chevauchent la nouvelle couture, quelle que soit leur latitude. Plutôt que de les découper proprement, ce qui demanderait de re-grider les données, le plus simple est de les retirer.
doc_clean <- doc |>
filter(lon != -0.5)On retrouve ici le réflexe du fond de carte de l’article précédent : une étendue continue (le rectangle globe d’un côté, la grille de DOC de l’autre) doit explicitement suivre la nouvelle couture, contrairement à une géométrie isolée comme un seul pays. La différence tient à la nature de l’objet : globe est un vrai polygone sf, qu’on découpe proprement avec st_difference() ; la grille de DOC reste un tableau classique tracé via geom_tile(), sans géométrie à découper, juste une case à exclure.
Continents et grille corrigés, il ne reste plus qu’à assembler le tout.
ggplot(doc_clean) +
geom_tile(aes(x = lon, y = lat, fill = doc_avg, colour = doc_avg)) +
geom_sf(data = world_robin, fill = "gray80", colour = NA) +
scale_fill_viridis_c(option = "F", trans = "log1p") +
scale_colour_viridis_c(option = "F", trans = "log1p", guide = "none") +
coord_sf(default_crs = sf::st_crs(4326)) +
labs(fill = "DOC (μmol kg<sup>-1</sup>)") +
theme_void() +
theme(legend.title = element_markdown())
C’est à peu de choses près la figure 1a de l’article : la vraie version y ajoute l’incertitude des prédictions et l’amplitude saisonnière, juste en dessous. Continents nets, plus aucune ligne parasite : le recentrage sur le Pacifique n’a plus rien d’un cauchemar.
Le recentrage n’est pas le seul terrain à pièges. Deux autres projections exotiques valent le détour : une vue satellite qui ampute des pays entiers, et une carte en écorce d’orange qui laisse des continents déborder de leurs lobes. Les symptômes n’ont rien à voir l’un avec l’autre, mais gratter un peu fait apparaître la même maladie : la reprojection naïve ne sait pas gérer une géométrie qui déborde du domaine que la projection peut représenter. Il faut la découper soi-même avant de la transformer, jamais après.
Le fond qui épouse la projection
Dans l’article précédent, on a construit globe, un rectangle en lat/lon qui, une fois projeté, épouse la forme réelle de la projection (ellipse, cercle…) plutôt que de rester un rectangle. On le réutilise tel quel pour Goode. La vue satellite, qui n’affiche qu’un seul hémisphère à la fois, demande sa propre variante, détaillée plus bas.
On réutilise world, déjà chargé en section précédente, et on recrée globe pour ces deux projections.
globe <- st_sfc(
st_polygon(list(cbind(
c(seq(-180, 180, 1), rep(180, 181), seq(180, -180, -1), rep(-180, 181)),
c(rep(-90, 361), seq(-90, 90, 1), rep(90, 361), seq(90, -90, -1))
))),
crs = "EPSG:4326"
)La projection orthographique simule une vue de la Terre depuis l’espace : elle n’affiche donc qu’un seul hémisphère à la fois. L’argument +lat_0 fixe la latitude du point central de la vue. Reprojetons world et globe sans traitement particulier.
ggplot(world) +
geom_sf(data = globe, fill = "#E4EFF6", color = NA) +
geom_sf(fill = "#6BAED6", color = "white", linewidth = 0.1) +
coord_sf(crs = "+proj=ortho +lat_0=45 +lon_0=10", expand = FALSE) +
theme_void() +
theme(legend.title = element_markdown())
Plusieurs pays en contact avec le bord de l’hémisphère visible ont purement disparu, comme les États-Unis ou l’Afrique du Sud : la projection orthographique ne sait pas proprement découper un pays qui chevauche exactement l’horizon. Les océans, eux, ne sont pas représentés du tout, pas de disque bleu clair pour délimiter les contours de la Terre : globe, pensé pour une projection qui couvre le monde entier sans interruption, ne produit ici qu’une géométrie vide côté face cachée.
On découpe avant de reprojeter, comme pour le Pacifique, plutôt que de compter sur coord_sf(). Ici on travaille directement avec sf : on construit d’abord le cercle de l’hémisphère visible, directement dans le CRS de la projection cible, son rayon correspond au rayon de la Terre pour une projection orthographique. C’est ce même cercle qui sert de fond de carte.
ortho_crs <- "+proj=ortho +lat_0=45 +lon_0=10"
# Cercle servant de fond de carte, construit directement dans le CRS cible
circle <- st_point(c(0, 0)) |>
st_buffer(dist = 6371000) |> # rayon de la Terre en m
st_sfc(crs = ortho_crs)
# Cercle servant à découper les pays : petite marge de sécurité (1000 m),
# repassé en lat/lon pour découper world qui y est encore exprimé
circle_clip <- st_point(c(0, 0)) |>
st_buffer(dist = 6371000 - 1000) |>
st_sfc(crs = ortho_crs) |>
st_transform(crs = "EPSG:4326")
world_visible <- st_intersection(st_make_valid(world), circle_clip)circle, déjà dans le bon CRS, sert directement de fond ; world_visible remplace world dans le geom_sf() suivant. st_make_valid() corrige les éventuelles géométries invalides avant la découpe, une précaution utile dès qu’on enchaîne avec st_intersection().
ggplot() +
geom_sf(data = circle, fill = "#E4EFF6", color = NA) +
geom_sf(data = world_visible, fill = "#6BAED6", color = "white", linewidth = 0.1) +
coord_sf(crs = ortho_crs, expand = FALSE) +
theme_void() +
theme(legend.title = element_markdown())
Pourquoi ne pas simplement réutiliser globe ?
globe fonctionne pour les projections qui affichent le monde entier sans interruption (Equal Earth, Robinson…), où chaque point du rectangle a un équivalent valide une fois projeté. Sur une projection à hémisphère unique, la moitié du rectangle correspond à la face cachée du globe, sans coordonnées valides dans cette projection : le résultat est une géométrie vide ou dégénérée. circle, lui, est construit directement dans les limites de ce qui est visible, ce problème ne se pose donc pas.
La projection interrompue de Goode préserve exactement les surfaces, comme Equal Earth, mais découpe le globe en lobes le long des océans, un peu comme une orange épluchée à plat. Reprojetons world et globe sans traitement particulier, comme pour l’orthographique.
ggplot() +
geom_sf(data = globe, fill = "#E4EFF6", color = NA) +
geom_sf(data = world, fill = "#6BAED6", color = "white", linewidth = 0.1) +
coord_sf(crs = "+proj=igh", expand = FALSE) +
theme_void() +
theme(legend.title = element_markdown())
Le Groenland et l’l’Antarctique débordent au-dessus des lobes, tandis que le fond de carte continue là où il devrait s’arrêter : globe, un rectangle continu, ne peut pas suivre les coupures de Goode.
Deux choses à corriger ici : le fond de carte d’abord, puis les débordements. On commence par un fond construit de façon directement ajustée aux méridiens de coupure.
goode_crs <- "+proj=igh"
# Contour construit directement aux méridiens de coupure de Goode : il
# épouse ainsi exactement la forme des lobes. On dessine le contour d'un
# seul trait, dans l'ordre : on descend le méridien 180° (hémisphère
# sud et nord réunis à l'est), puis on remonte et redescend à chaque
# coupure de l'hémisphère sud (océan Indien ~80°E, Atlantique Sud ~20°O,
# Pacifique Sud ~100°O), on longe le méridien -180° à l’ouest pour rejoindre
# l'hémisphère nord, puis on referme avec la seule coupure au nord
# (Atlantique Nord ~40°O) avant de revenir au point de départ.
lats <- c(90:-90, -90:0, 0:-90, -90:0, 0:-90, -90:0, 0:-90, -90:90, 90:0, 0:90, 90)
longs <- c(
rep(180, 181),
rep(c(80.01, 79.99), each = 91),
rep(c(-19.99, -20.01), each = 91),
rep(c(-99.99, -100.01), each = 91),
rep(-180, 181),
rep(c(-40.01, -39.99), each = 91),
180
)
goode_outline <- list(cbind(longs, lats)) |>
st_polygon() |>
st_sfc(crs = "EPSG:4326") |>
st_transform(crs = goode_crs)ggplot() +
geom_sf(data = goode_outline, fill = "#E4EFF6", color = NA) +
geom_sf(data = world, fill = "#6BAED6", color = "white", linewidth = 0.1) +
coord_sf(crs = goode_crs, expand = FALSE) +
theme_void() +
theme(legend.title = element_markdown())
Les lobes sont nets, mais le Groenland et l’Antarctique chevauchent toujours les coupures. Plutôt que de découper ces deux pays précisément à cet endroit, la solution la plus simple consiste à masquer les débordements après coup.
# Masque, de la couleur de fond de la page, posé par-dessus tout ce qui
# déborde des lobes (Groenland, Antarctique) pour le cacher plutôt que
# d'essayer de découper ces pays eux-mêmes
goode_xlim <- st_bbox(goode_outline)[c("xmin", "xmax")] * 1.1
goode_ylim <- st_bbox(goode_outline)[c("ymin", "ymax")] * 1.1
# Rectangle englobant, un peu plus large que le contour (marge de 10 %)
# pour couvrir tout débordement possible
goode_encl_rect <- list(cbind(
c(goode_xlim[1], goode_xlim[2], goode_xlim[2], goode_xlim[1], goode_xlim[1]),
c(goode_ylim[1], goode_ylim[1], goode_ylim[2], goode_ylim[2], goode_ylim[1])
)) |>
st_polygon() |>
st_sfc(crs = goode_crs)
# Masque = tout ce qui est dans le rectangle mais hors du contour des lobes
goode_mask <- st_difference(goode_encl_rect, goode_outline)ggplot() +
geom_sf(data = goode_outline, fill = "#E4EFF6", color = NA) +
geom_sf(data = world, fill = "#6BAED6", color = "white", linewidth = 0.1) +
geom_sf(data = goode_mask, fill = "white", color = NA) +
coord_sf(crs = goode_crs, expand = FALSE) +
theme_void() +
theme(legend.title = element_markdown())
Pourquoi un masque plutôt que découper le Groenland et l’Antarctique ?
Un masque de la couleur de fond est plus simple que de découper le Groenland et l’Antarctique à l’endroit exact des coupures, comme on l’a fait pour l’hémisphère visible de l’orthographique : il cache le débordement sans toucher aux données, et s’adapte tout seul si on change un jour de fond de carte ou de résolution.
Trois projections, trois symptômes différents en apparence, mais un même mécanisme en toile de fond.
| Projection | Symptôme | Cause | Solution |
|---|---|---|---|
| Robinson recentré (Pacifique) | Continents qui se chevauchent, lignes parasites dans la grille de données | La nouvelle couture (fixée par lon_0) tombe sur des géométries jamais pensées pour être coupées à cet endroit |
Découper polygones et grille à la nouvelle couture avant de reprojeter, jamais après |
| Orthographique (vue satellite) | Pays en bordure d’hémisphère amputés, fond de carte manquant | coord_sf() ne sait pas gérer une géométrie qui déborde de l’hémisphère visible |
Construire le cercle de l’hémisphère visible et découper les données avec, avant de reprojeter |
| Goode homolosine (écorce d’orange) | Continents qui débordent des lobes, fond de carte qui continue hors des coupures | Un fond rectangulaire générique ne peut pas suivre les coupures de la projection | Fond de carte calé sur les méridiens de coupure, complété d’un masque pour les débordements résiduels |
Le réflexe à garder n’est donc pas une recette différente par projection, mais une question à se poser à chaque fois : où se trouve la nouvelle limite ou la nouvelle couture, et mes données la traversent-elles ? Si oui, elles doivent être découpées ou filtrées avant la reprojection, jamais confiées telles quelles à coord_sf().
En pratique, ça tient en trois réflexes, toujours dans le même ordre :
Ce qui m’a pris une journée entière la première fois (voir le tout début de cet article) tient maintenant en quelques lignes. Ça n’aurait pas rendu la découverte initiale moins pénible sur le moment, mais ça aide à se souvenir que le problème est presque toujours le même, juste déguisé différemment.
Si Mercator, Equal Earth et Robinson suffisent déjà à la plupart de vos usages, ce qui est le cas la majorité du temps, direction l’article précédent pour les bases.
Votre carte refuse de coopérer ? On peut sûrement arranger ça.
la longitude où la carte se referme sur elle-même, située à l’opposé du méridien central↩︎