Chapitre 8 Les reprojections

Dans ce chapitre, nous allons utiliser les librairies suivantes.

library(sf)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(lwgeom)

Pour rappel, il existe deux types de CRS, les CRS géographiques (longitude/lattitude avec pour unité de compte des degrés) et les CRS projetés (avec un datum et une unité en mètre par exemple).

La plupart des fonctions de sf présupposent s’appliquer sur un CRS projeté, car les fonctions de GEOS sur lesquelles elles se basent le font aussi.

8.1 Un premier exemple de reprojection

Prenons les coordonnées de Nantes en WGS 84:

nantes<-data.frame(lon=-1.553621,lat=47.218371) %>% 
    st_as_sf(coords = c("lon", "lat")) %>% 
  st_set_crs(4326)

On peut visualiser nos données

mapview(nantes)

st_is_longlat() est une fonction de sf qui permet de faire un test sur la famille de CRS à laquelle on a à faire.

st_is_longlat(nantes)
[1] TRUE

Essayons de créer un buffer de 1km autour de Nantes

nantes_buffer= st_buffer(nantes, dist = 1)
Warning in st_buffer.sfc(st_geometry(x), dist, nQuadSegs, endCapStyle = endCapStyle, : st_buffer does not
correctly buffer longitude/latitude data
mapview(list(nantes,nantes_buffer))

Le message est clair et indique qu’un buffer ne marchera pas correctement sur des données en projection longitude / lattitude.

Tentons une reprojection. La fonction permettant une reprojection est st_transform(). On va ici passer en lambert 93 nos données.

nantes_proj<-st_transform(nantes, 2154)

Le CRS lambert 93 est bien un CRS projeté :

st_is_longlat(nantes_proj)
[1] FALSE
nantes_proj_buffer= st_buffer(nantes_proj, dist = 1)
mapview(list(nantes_proj,nantes_proj_buffer))

8.2 Quand reprojeter ?

Quelques cas usuels qui peuvent vous amener à reprojeter vos données :

  • la manipulation de données fournies dans des CRS différents

  • l’usage du package leaflet impose des données spécifiées en WGS 84

  • le besoin de visualiser vos données suivant la conversion de certaines propriétés des objets à la surface de la terre.

  • l’usage de fonctions demandant à utiliser des CRS projetés (comme st_buffer() ci-dessus)

Un exemple d’usage : la distance de Rennes à Nantes

Prenons les coordonnées WGS 84 de Rennes

rennes<-data.frame(lon=-1.6777926,lat=48.117266) %>% 
    st_as_sf(coords = c("lon", "lat")) %>% 
  st_set_crs(4326)
mapview(rennes)

Tentons de calculer la distance de Rennes à Nantes. Avec la données en Lambert 93, la fonction st_distance() renvoie un message d’erreur.

st_distance(rennes,nantes_proj)
Error in st_distance(rennes, nantes_proj): st_crs(x) == st_crs(y) is not TRUE

Avec la données en WGS 84, la fonction st_distance() renvoie bien le résultat.

st_distance(rennes,nantes)
Units: [m]
         [,1]
[1,] 100376.7

8.3 Quel CRS utiliser ?

A cette question, il y a rarement une bonne réponse.

En ce qui concerne les CRS géométriques, le plus simple est d’utiliser le WGS 84, qui est de loin le plus populaire, avec lequel beaucoup de données sont fournies.

En ce qui concerne les CRS projetés, utiliser par défaut le lambert 93, le CRS officiel français, pour les données nationales fait sens.

Ensuite votre choix va dépendre des propriétés que vous souhaitez conserver.

8.4 Comment projeter ?

8.4.1 Projeter des vecteurs

Reprojeter des données vecteur se fait à l’aide de la fonction st_transform() que nous avons vu en utilisant le code epsg que nous voulons.

8.4.2 Modifier la projection d’une carte.

Parfois on souhaite pouvoir aller plus loin dans les reprojections, en adaptant le centre de la projection, pour cela on peut utiliser un proj4string ad hoc.

Pour cela, on va modifier l’argument +proj de notre crs avec st_transform.

Tentons par exemple de reprojeter notre carte du globe en utilisant la projection azimutale équivalente de Lambert, centrée sur Pékin.

data("World")
World_pekin<-st_transform(World,crs = "+proj=laea +x_0=0 +y_0=0 +lon_0=116 +lat_0=40")

Le paramètre +proj=laea permet de redéfinir la projection, les paramètres +lon_0 et lat_0 permettent de définir le centre de la projection. x_0 et y_0 définissent le centre du plan pour les coordonnées.

Qu’est ce qui a changé entre nos deux cartes ?

st_crs(World)
Coordinate Reference System:
  User input: +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Eckert IV"],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
st_crs(World_pekin)
Coordinate Reference System:
  User input: +proj=laea +x_0=0 +y_0=0 +lon_0=116 +lat_0=40 
  wkt:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Lambert Azimuthal Equal Area",
            ID["EPSG",9820]],
        PARAMETER["Latitude of natural origin",40,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",116,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]