sp :: proj4string equivalente nel framework PROJ6?

Question

Sono confuso circa il passaggio da PROJ4 a PROJ6 nelle versioni (relativamente) recenti dei pacchetti spaziali di R (e gdal per quella materia). Usavo sp::proj4string()per ottenere il codice di proiezione / argomento CRS da un spoggetto e utilizzare ulteriormente quel codice nei miei script. Ora ricevo un avviso (che mi dice che la mia proiezione potrebbe essere imprecisa ). Ovviamente sto facendo qualcosa di sbagliato. Quale funzione devo utilizzare per ottenere il codice di proiezione senza preavviso? Personalmente preferirei il codice EPSG se quel livello di semplicità fosse possibile ...

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"

^{Creato il 18/09/2020 dal pacchetto reprex (v0.3.0)}

Vedi anche i riferimenti e i commenti in questa domanda

Answer 1

Ero anche abbastanza confuso riguardo all'allontanamento da PROJ Strings e cosa usare invece fino a quando ho visto questo video di un discorso tenuto l'anno scorso a FOSS4G.

Ecco le mie principali conclusioni:

Perché PROJ e GDAL si stanno allontanando dalle stringhe PROJ

Con PROJ4 ogni trasformazione viene eseguita utilizzando l '"approccio hub". Quindi tutto verrà prima convertito in WGS84 e solo successivamente trasformato nella proiezione target.

Questo approccio può portare a errori più grandi nel processo di trasformazione rispetto alla trasformazione diretta tra le due proiezioni. Nel peggiore dei casi la trasformazione necessaria non è nemmeno possibile con questo approccio. Poiché le stringhe PROJ possono essere ambigue, le trasformazioni dirette spesso non sono possibili.

Un esempio dal discorso:

Il vecchio dato australiano GDA94 e il nuovo GDA2020 sono entrambi basati su WGS84. Tuttavia, uno utilizza il frame di riferimento ITRF2014 a Epoch 2020.0 e l'altro ITRF92 a Epoch 1994.0. Queste differenze non possono essere catturate in una stringa PROJ (per quanto ho capito) e quindi tutte le proiezioni hanno la stessa stringa PROJ quasi identica di WGS84. La conversione tra entrambi con l'approccio hub non porterebbe a nessun cambiamento nelle coordinate. Mentre in realtà l'Australia si è spostata di 1,8 metri tra il 1994 e il 2020 e la conversione del dato dovrebbe riflettere questo.

Con una trasformazione diretta errori del genere non sarebbero accaduti.

Cosa dovresti usare invece

Invece dovresti idealmente usare definizioni non ambigue per proiezioni come codici EPSG o definizioni in formato WKT2.

Tuttavia, se non si intende effettuare trasformazioni di dati (o quindi non si preoccupano di pochi metri o più di imprecisione) è comunque possibile utilizzare le stringhe PROJ, ma si tenga presente che la migliore pratica è utilizzare i codici EPSG o WKT2.

Come viene implementato in sp

CRS-Gli oggetti in sp hanno ora un campo aggiuntivo commentin cui viene salvata una rappresentazione WKT2 del SR.

Quindi il flusso di lavoro delle migliori pratiche con spora è simile a questo (preso e ampliato un po 'da qui ):

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

Tutto sommato sfha funzioni molto più convenienti per l'accesso e la gestione di crs rispetto a sp. Puoi leggere di più su come lo sfgestisce nello stesso link sopra.

Answer 2

Proj4strings è in fase di eliminazione. Grazie al pacchetto sf e ai suoi sviluppatori c'è un semplice trucco per convertire una proj4string nel formato WKT.

Basta ignorare gli avvisi di dato provenienti da rgdal e, infatti, puoi semplicemente disattivarli perché diventano fastidiosi. Aggiungilo options("rgdal_show_exportToProj4_warnings"="none")al tuo file Rprofile.site o emettilo prima di aggiungere qualsiasi libreria.

Ad esempio con la proiezione geografica WGS84, l'utilizzo sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")risulterà 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]]]]

Dal blog R-spatial R spatial segue GDAL e lo sviluppo di PROJ gli sviluppatori dichiarano:

L'uso delle cosiddette stringhe PROJ4 (come + proj = longlat + datum = WGS84) è sconsigliato, poiché non offrono più una descrizione sufficiente dei sistemi di riferimento delle coordinate; l'uso di + init = epsg: XXXX porta ad avvertimenti

Answer 3

Quell'avvertimento non influisce affatto sul tuo output. Come hanno suggerito le altre persone, il pacchetto sf può essere utile per una codifica EPSG semplice e pulita. Ecco il tuo codice adattato:

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)