odpowiednik sp :: proj4string we frameworku PROJ6?

Question

Jestem zdezorientowany przejściem z PROJ4 do PROJ6 w (stosunkowo) najnowszych wersjach pakietów przestrzennych R (i gdal o to chodzi). Kiedyś używałem go sp::proj4string()do pobierania kodu projekcji / argumentu CRS z spobiektu i używania tego kodu dalej w moich skryptach. Teraz otrzymuję ostrzeżenie (które mówi mi, że moja prognoza może być niedokładna ). Oczywiście robię coś złego. Której funkcji należy użyć, aby uzyskać kod projekcji bez ostrzeżenia? Osobiście wolałbym kod EPSG, gdyby taki poziom prostoty był możliwy ...

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"

^{Utworzono 18.09.2020 przez pakiet reprex (v0.3.0)}

Zobacz także odniesienia i komentarze w tym pytaniu

Answer 1

Byłem również dość zdezorientowany odejściem od PROJ Strings i tym, czego zamiast tego użyć, dopóki nie obejrzałem tego wideo z wykładu, który odbył się w zeszłym roku na FOSS4G.

Oto moje główne wnioski:

Dlaczego PROJ i GDAL odchodzą od łańcuchów PROJ

W przypadku PROJ4 każda transformacja jest wykonywana przy użyciu podejścia „koncentratora”. Więc wszystko zostanie najpierw przekonwertowane na WGS84, a dopiero potem przekształcone w projekcję docelową.

Takie podejście może prowadzić do większych błędów w procesie transformacji niż bezpośrednia transformacja między dwiema projekcjami. W najgorszym przypadku potrzebna transformacja nie jest nawet możliwa przy takim podejściu. Ponieważ łańcuchy PROJ mogą być niejednoznaczne, bezpośrednie transformacje często nie są możliwe.

Przykład z wykładu:

Stare australijskie dane GDA94 i nowe GDA2020 są oparte na WGS84. Jednak jeden używa ramki odniesienia ITRF2014 w Epoce 2020.0, a drugi ITRF92 w Epoce 1994.0. Tych różnic nie można uchwycić w ciągu PROJ (o ile zrozumiałem), a zatem wszystkie projekcje mają prawie dokładnie taki sam ciąg PROJ jak WGS84. Konwersja między obydwoma z podejściem typu hub prowadziłaby do braku zmiany współrzędnych. Chociaż w rzeczywistości Australia przesunęła się o 1,8 metra między 1994 a 2020 r., A konwersja odniesienia powinna to odzwierciedlać.

Z bezpośrednią transformacją takie błędy nie miałyby miejsca.

Co powinieneś zamiast tego użyć

Zamiast tego najlepiej byłoby używać jednoznacznych definicji dla odwzorowań, takich jak kody EPSG lub definicje w formacie WKT2.

Jeśli jednak nie zamierzasz wykonywać transformacji układu odniesienia (lub nie przejmujesz się kilkoma metrami lub więcej niedokładnościami), nadal możesz używać łańcuchów PROJ, ale pamiętaj, że najlepszą praktyką jest użycie kodów EPSG lub WKT2.

Jak jest to realizowane w sp

CRS-obiekty w sp mają teraz dodatkowe pole, commentw którym zapisywana jest reprezentacja CRS WKT2.

Tak więc najlepszy przepływ pracy z spteraz wygląda następująco (wzięty i rozwinięty trochę stąd ):

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

W sumie sfma o wiele wygodniejsze funkcje dostępu i obsługi crs niż sp. Możesz przeczytać więcej o tym, jak sfsobie z tym radzi w tym samym linku podanym powyżej.

Answer 2

Proj4strings są wycofywane. Dzięki pakietowi sf i jego twórcom istnieje prosta sztuczka do konwersji proj4string do formatu WKT.

Po prostu zignoruj ​​ostrzeżenia dotyczące danych pochodzące z rgdal i faktycznie możesz je po prostu wyłączyć, ponieważ stają się uciążliwe. Dodaj options("rgdal_show_exportToProj4_warnings"="none")do pliku Rprofile.site lub wydaj go przed dodaniem jakichkolwiek bibliotek.

Na przykład w przypadku odwzorowania geograficznego WGS84 użycie sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")da w wyniku

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

Z bloga R-spatial R spatial podąża za deweloperami sf GDAL i PROJ development.

Używanie tak zwanych łańcuchów PROJ4 (takich jak + proj = longlat + datum = WGS84) jest odradzane, ponieważ nie oferują one już wystarczającego opisu układów odniesienia za pomocą współrzędnych; użycie + init = epsg: XXXX prowadzi do ostrzeżeń

Answer 3

To ostrzeżenie w ogóle nie wpływa na wynik. Jak sugerowali inni ludzie, pakiet sf może być przydatny do czystszego i prostego kodowania EPSG. Oto dostosowany kod:

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)