sp :: proj4string-Äquivalent im PROJ6-Framework?

Question

Ich bin verwirrt über die Verschiebung von PROJ4 zu PROJ6 in den (relativ) neueren Versionen der räumlichen Pakete von R (und gdal). Früher habe ich verwendet, sp::proj4string()um das Projektionscode- / CRS-Argument von einem spObjekt abzurufen und diesen Code in meinen Skripten weiter zu verwenden. Jetzt erhalte ich eine Warnung (die mir sagt, dass meine Projektion möglicherweise ungenau ist ). Ich mache offensichtlich etwas falsch. Mit welcher Funktion kann ich den Projektionscode ohne Warnung abrufen? Ich persönlich würde den EPSG-Code bevorzugen, wenn diese Einfachheit möglich wäre ...

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"

^{Erstellt am 2020-09-18 durch das reprex-Paket (v0.3.0)}

Siehe auch die Referenzen und Kommentare in dieser Frage

Answer 1

Ich war auch ziemlich verwirrt über die Abkehr von PROJ Strings und was ich stattdessen verwenden sollte, bis ich dieses Video eines Vortrags sah, der letztes Jahr bei FOSS4G stattfand.

Hier sind meine wichtigsten Erkenntnisse:

Warum PROJ und GDAL sich von PROJ-Strings entfernen

Mit PROJ4 erfolgt jede Transformation nach dem "Hub-Ansatz". Also wird alles zuerst in WGS84 konvertiert und erst dann in die Zielprojektion transformiert.

Dieser Ansatz kann zu größeren Fehlern im Transformationsprozess führen als die direkte Transformation zwischen den beiden Projektionen. Im schlimmsten Fall ist die erforderliche Transformation mit diesem Ansatz nicht einmal möglich. Da PROJ-Strings mehrdeutig sein können, sind direkte Transformationen oft nicht möglich.

Ein Beispiel aus dem Vortrag:

Das alte australische Datum GDA94 und das neue GDA2020 basieren beide auf WGS84. Einer verwendet jedoch den Referenzrahmen ITRF2014 in Epoche 2020.0 und der andere ITRF92 in Epoche 1994.0. Diese Unterschiede können (soweit ich verstanden habe) nicht in einem PROJ-String erfasst werden und daher haben alle Projektionen fast genau den gleichen PROJ-String wie WGS84. Das Konvertieren zwischen beiden mit dem Hub-Ansatz würde zu keiner Änderung der Koordinaten führen. Während sich Australien zwischen 1994 und 2020 tatsächlich um 1,8 Meter bewegt hat, sollte die Bezugsumrechnung dies widerspiegeln.

Mit einer direkten Transformation würden solche Fehler nicht passieren.

Was Sie stattdessen verwenden sollten

Stattdessen sollten Sie idealerweise nicht mehrdeutige Definitionen für Projektionen wie EPSG-Codes oder Definitionen im WKT2-Format verwenden.

Wenn Sie jedoch nicht beabsichtigen, Bezugstransformationen durchzuführen (oder sich dann nicht um einige Meter oder mehr Ungenauigkeit kümmern), können Sie weiterhin PROJ-Strings verwenden. Beachten Sie jedoch, dass die Verwendung von EPSG-Codes oder WKT2 die beste Vorgehensweise ist.

Wie wird dies in sp

CRS-Objekte in sp haben jetzt ein zusätzliches Feld, commentin dem eine WKT2-Darstellung des CRS gespeichert wird.

Der Best-Practice-Workflow mit spjetzt sieht also so aus (von hier aus ein wenig übernommen und erweitert ):

# 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)

Alles in allem sfhat es viel bequemere Funktionen für den Zugriff auf und die Handhabung von crs als sp. Weitere Informationen zum Umgang damit finden sfSie unter dem oben angegebenen Link.

Answer 2

Proj4strings werden auslaufen. Dank des sf-Pakets und seiner Entwickler gibt es einen einfachen Trick zum Konvertieren eines proj4string in das WKT-Format.

Ignorieren Sie einfach die von rgdal kommenden Datumswarnungen, und Sie können sie einfach ausschalten, da sie lästig werden. options("rgdal_show_exportToProj4_warnings"="none")Fügen Sie sie entweder zu Ihrer Rprofile.site-Datei hinzu oder geben Sie sie aus, bevor Sie Bibliotheken hinzufügen.

Zum Beispiel mit der geographischen WGS84 - Projektion, indem sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")führt in

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]]]]

Aus der R-Raum Blog R räumlichen folgt GDAL und PROJ Entwicklung sf Entwickler Zustand:

Von der Verwendung sogenannter PROJ4-Strings (wie + proj = longlat + datum = WGS84) wird abgeraten, sie bieten keine ausreichende Beschreibung von Koordinatenreferenzsystemen mehr. Die Verwendung von + init = epsg: XXXX führt zu Warnungen

Answer 3

Diese Warnung wirkt sich überhaupt nicht auf Ihre Ausgabe aus. Wie die anderen Leute vorgeschlagen haben, kann das sf-Paket für eine sauberere und einfache EPSG-Codierung nützlich sein. Hier ist Ihr Code angepasst:

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)