sp :: proj4string equivalente no framework PROJ6?

Question

Estou confuso sobre a mudança de PROJ4 para PROJ6 nas versões (relativamente) recentes dos pacotes espaciais de R (e gdal para esse assunto). Eu costumava usar sp::proj4string()para obter o código de projeção / argumento CRS de um spobjeto e usar esse código posteriormente em meus scripts. Agora recebo um aviso (que me diz que minha projeção pode ser imprecisa ). Obviamente, estou fazendo algo errado. Qual função devo usar para obter o código de projeção sem um aviso? Eu pessoalmente preferiria o código EPSG se esse nível de simplicidade fosse possível ...

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"

^{Criado em 2020-09-18 pelo pacote reprex (v0.3.0)}

Veja também as referências e comentários nesta questão

Answer 1

Eu também estava bastante confuso sobre como deixar de usar PROJ Strings e o que usar no lugar até assistir a este vídeo de uma palestra realizada no ano passado na FOSS4G.

Aqui estão minhas principais lições:

Por que PROJ e GDAL estão se afastando das strings PROJ

Com o PROJ4, todas as transformações são feitas usando a "abordagem hub". Portanto, tudo será primeiro convertido para WGS84 e só então transformado para a projeção de destino.

Essa abordagem pode levar a erros maiores no processo de transformação do que a transformação direta entre as duas projeções. Na pior das hipóteses, a transformação necessária nem mesmo é possível com essa abordagem. Visto que as strings PROJ podem ser ambíguas, as transformações diretas geralmente não são possíveis.

Um exemplo da palestra:

O antigo datum australiano GDA94 e o novo GDA2020 são ambos baseados no WGS84. No entanto, um está usando o quadro de referência ITRF2014 na Epoch 2020.0 e o outro ITRF92 na Epoch 1994.0. Essas diferenças não podem ser capturadas em uma string PROJ (tanto quanto eu entendi) e, portanto, todas as projeções têm quase exatamente a mesma string PROJ que WGS84. A conversão entre ambos com a abordagem de hub não levaria a nenhuma mudança nas coordenadas. Enquanto na verdade a Austrália moveu 1,8 metros entre 1994 e 2020 e a conversão do datum deve refletir isso.

Com uma transformação direta, erros como esse não aconteceriam.

O que você deve usar no lugar

Em vez disso, você deve usar definições não ambíguas para projeções como códigos EPSG ou definições no formato WKT2.

No entanto, se você não pretende fazer transformações de datum (ou então não se importa com alguns metros ou mais de imprecisão), você ainda pode usar PROJ-Strings, mas esteja ciente de que a melhor prática é usar códigos EPSG ou WKT2.

Como isso é implementado em sp

CRS-objetos em sp agora têm um campo adicional commentno qual uma representação WKT2 do CRS é salva.

Portanto, o fluxo de trabalho de práticas recomendadas spagora se parece com este (tomado e expandido um pouco a partir daqui ):

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

Em suma, sftem funções muito mais convenientes para acessar e manipular crs do que sp. Você pode ler mais sobre como sflidar com isso no mesmo link fornecido acima.

Answer 2

Proj4strings estão sendo eliminados. Graças ao pacote sf e seus desenvolvedores, há um truque simples para converter um proj4string para o formato WKT.

Simplesmente ignore os avisos de datum vindos de rgdal e, na verdade, você pode simplesmente desligá-los porque eles se tornam enfadonhos. Adicione options("rgdal_show_exportToProj4_warnings"="none")ao seu arquivo Rprofile.site ou emita-o antes de adicionar quaisquer bibliotecas.

Por exemplo, com a projeção geográfica WGS84, o uso sf::st_crs("+proj=longlat +datum=WGS84 +no_defs")resultará em

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

Do blog R-spatial R spatial segue o estado dos desenvolvedores sf de desenvolvimento GDAL e PROJ :

O uso das chamadas strings PROJ4 (como + proj = longlat + datum = WGS84) são desencorajados, eles não oferecem mais uma descrição suficiente dos sistemas de referência de coordenadas; uso de + init = epsg: XXXX leva a avisos

Answer 3

Esse aviso não afeta sua saída. Como as outras pessoas sugeriram, o pacote sf pode ser útil para uma codificação EPSG mais limpa e simples. Aqui está o seu código adaptado:

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)