sp :: proj4string équivalent dans le framework PROJ6?

Question

Je suis confus par le passage de PROJ4 à PROJ6 dans les versions (relativement) récentes des packages spatiaux de R (et gdal d'ailleurs). J'avais l'habitude d' sp::proj4string()obtenir le code de projection / l'argument CRS d'un spobjet et d'utiliser ce code plus loin dans mes scripts. Maintenant, je reçois un avertissement (qui me dit que ma projection pourrait être inexacte ). Je fais évidemment quelque chose de mal. Quelle fonction dois-je utiliser pour obtenir le code de projection sans avertissement? Personnellement, je préférerais le code EPSG si ce niveau de simplicité était possible ...

library(sp)

x <- SpatialPolygons(
            list(Polygons(
                 list(Polygon(cbind(c(1,2,3,4,1), c(3,8,10,12,3)))), 
                 ID = 1)), 
            proj4string = CRS(SRS_string = "EPSG:4326"))

proj4string(x)
#> Warning in proj4string(x): CRS object has comment, which is lost in output
#> [1] "+proj=longlat +datum=WGS84 +no_defs"

^{Créé le 18/09/2020 par le package reprex (v0.3.0)}

Voir aussi les références et commentaires dans cette question

Answer 1

J'étais également assez confus quant à l'abandon de PROJ Strings et à ce qu'il fallait utiliser à la place jusqu'à ce que je regarde cette vidéo d'une conférence tenue l'année dernière à FOSS4G.

Voici mes principaux points à retenir:

Pourquoi PROJ et GDAL s'éloignent des chaînes PROJ

Avec PROJ4, chaque transformation est effectuée en utilisant «l'approche hub». Ainsi, tout sera d'abord converti en WGS84 et ensuite seulement transformé en projection cible.

Cette approche peut conduire à des erreurs plus importantes dans le processus de transformation que la transformation directe entre les deux projections. Dans le pire des cas, la transformation nécessaire n'est même pas possible avec cette approche. Étant donné que les chaînes PROJ peuvent être ambiguës, les transformations directes ne sont souvent pas possibles.

Un exemple de la conférence:

L'ancienne datum australienne GDA94 et la nouvelle GDA2020 sont toutes deux basées sur WGS84. Cependant, l'un utilise le cadre de référence ITRF2014 à Epoch 2020.0 et l'autre ITRF92 à Epoch 1994.0. Ces différences ne peuvent pas être capturées dans une chaîne PROJ (pour autant que je sache) et donc toutes les projections ont presque exactement la même chaîne PROJ que WGS84. La conversion entre les deux avec l'approche hub n'entraînerait aucun changement dans les coordonnées. Alors qu'en fait, l'Australie a déplacé 1,8 mètre entre 1994 et 2020 et la conversion du datum devrait en tenir compte.

Avec une transformation directe, des erreurs comme celles-là ne se produiraient pas.

Ce que vous devriez utiliser à la place

Au lieu de cela, vous devriez idéalement utiliser des définitions non ambiguës pour des projections telles que les codes EPSG ou des définitions au format WKT2.

Cependant, si vous n'avez pas l'intention de faire des transformations de données (ou que vous ne vous souciez pas de quelques mètres ou plus d'inexactitude), vous pouvez toujours utiliser PROJ-Strings, mais sachez que la meilleure pratique consiste à utiliser les codes EPSG ou WKT2.

Comment est-ce implémenté dans sp

CRS-objets dans sp ont maintenant un champ supplémentaire commentdans lequel une représentation WKT2 du CRS est enregistrée.

Ainsi, le flux de travail des meilleures pratiques avec spmaintenant ressemble à ceci (pris et développé un peu à partir d' ici ):

# Define the CRS using an EPSG Code
x <- CRS(SRS_string='EPSG:4326')

# Display the stored CRS using comment()
cat(comment(x), "\n")

# Store the wkt in a variable
wkt <- comment(x)

# Use this to assign the CRS of another sp-object
y <- CRS(SRS_string = wkt)

Dans l'ensemble, il sfa des fonctions beaucoup plus pratiques pour accéder et gérer les crs que sp. Vous pouvez en savoir plus sur la façon dont le sfgère dans le même lien ci-dessus.

Answer 2

Proj4strings est en cours de suppression. Grâce au package sf et à ses développeurs, il existe une astuce simple pour convertir une chaîne proj4 au format WKT.

Ignorez simplement les avertissements de datum provenant de rgdal et, en fait, vous pouvez simplement les désactiver car ils deviennent ennuyeux. Ajoutez options("rgdal_show_exportToProj4_warnings"="none")à votre fichier Rprofile.site ou émettez-le avant d'ajouter des bibliothèques.

Par exemple, avec la projection géographique WGS84, l'utilisation sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")se traduira par

Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +no_defs 
  wkt:
GEOGCRS["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]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]]

À partir du blog R-spatial, R spatial suit l'état des développeurs GDAL et PROJ development sf:

L'utilisation de chaînes dites PROJ4 (comme + proj = longlat + datum = WGS84) est déconseillée, elles n'offrent plus une description suffisante des systèmes de référence de coordonnées; l'utilisation de + init = epsg: XXXX entraîne des avertissements

Answer 3

Cet avertissement n'affecte pas du tout votre sortie. Comme les autres l'ont suggéré, le package sf peut être utile pour un codage EPSG plus propre et simple. Voici votre code adapté:

library(sf)
library(mapview)
coords <- list(rbind(c(1,3), c(2,8), c(3,10), c(4,12), c(1,3)))
x2 <- st_sf(st_geometry(st_polygon(x = coords)), crs="EPSG:4326")
x2$ID <- 1
mapview(x2)