Alternativas do spTransform?

8

Digamos que tenhamos um shapefile com uma certa projeção.

s<-readOGR(dsn=".",layer="Spain")

Também temos os aeroportos da Espanha como pontos em outra projeção.

a<-readOGR(dsn=".",layer="airports")

Se pretendemos colocar os pontos no shapefile da Espanha, temos que organizar as coordenadas para que sejam iguais. Normalmente, é feito assim:

 a<-spTransform(a,CRS(proj4string(s))

Mas é o mesmo com isso?

proj4string(a)<-proj4string(s)

Se sim, então por que não é a maneira padrão de fazer isso, pois é mais simples e precisa recorrer ao spTransform?

gsa
fonte

Respostas:

6

A atribuição de um CRS diferente não altera a projeção dos dados espaciais subjacentes - o CRS é uma parte interna do objeto espacial que indica a R como interpretar as coordenadas espaciais.

library(rgdal)

# Load Tanzania in UTM 36
tz.36 <- readOGR(dsn = ".", layer = "tz_36")
summary(tz.36) # Show the bounding coordinates:

Coordinates:
         min     max
x   94000.58 1315978
y 8699697.87 9889701
Is projected: TRUE 
proj4string :
[+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]
....

Agora transforme a forma em uma zona UTM vizinha:

tz.37 <- spTransform(tz.36, CRS("+init=epsg:32737"))
summary(tz.37)

Coordinates:
        min       max
x -576091.7  657248.1
y 8700995.0 9888925.5
Is projected: TRUE 
proj4string :
[+init=epsg:32737 +proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]

As coordenadas da geometria foram alteradas - corretamente - para se ajustarem à zona UTM vizinha. E se simplesmente redesignássemos o CRS sem transformação?

# Make a copy of the original object.
tz.37.b <- tz.36 
# Assign the CRS
proj4string(tz.37.b) <- CRS("+init=epsg:32737") 

Warning message:
In `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) :
  A new CRS was assigned to an object with an existing CRS:
+proj=utm +zone=36 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
without reprojecting.
For reprojection, use function spTransform in package rgdal

Fomos avisados ​​... então, como são as coordenadas?

summary(tz.37.b)

    Coordinates:
         min     max
x   94000.58 1315978
y 8699697.87 9889701
Is projected: TRUE 
proj4string :
[+init=epsg:32737 +proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]

Os números (coordenadas) são os originais da zona UTM 36, mas a forma será mapeada para a zona UTM incorreta agora e aparecerá no lugar errado.

Editar

Usando o método exato sugerido pelo OP:

# Make a new copy of the original UTM 36 shape:
tz.37.c <- tz.36
# Now assign the proj4string using OP's suggestion:
proj4string(tz.37.c) <- proj4string(tz.37)
summary(tz.37.c)

Coordinates:
         min     max
x   94000.58 1315978
y 8699697.87 9889701
Is projected: TRUE 
proj4string :
[+init=epsg:32737 +proj=utm +zone=37 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]

A extensão é igual ao objeto utm36 original, mas o proj4string agora é utm37. (Curiosamente, não houve aviso neste momento.) Nota: o resultado é EXATAMENTE o mesmo que ter atribuído o CRS com o código EPSG. Testar:

identical(tz.37.b, tz.37.c)
TRUE

Vamos verificar como são as formas reais, usando uma camada UTM37 real das regiões da Tanzânia.

# Load the region shape
tz.regions.37 <- readOGR(".", "tz_regions_37_simp")
# Plot it
plot(tz.regions.37, lwd = 1, border = 'red')
# Now add the CRS-assigned (not transformed!) object:
plot(tz.37.c, add = T, border = 'blue', lwd = 2)

insira a descrição da imagem aqui

Eles não se alinham, então esse método não está funcionando! E o objeto transformado ?

plot(tz.37, add = T, border = 'darkgreen', lwd = 2)

insira a descrição da imagem aqui

O objeto transformado aparece no lugar certo (fundo verde escuro).

NB Se você estiver transformando entre zonas UTM, mas o mesmo dado (WGS84 aqui), as diferenças serão dramáticas, como no meu exemplo - mas se você estiver transformando entre dados diferentes, as diferenças poderão ser muito mais sutis. Por exemplo, abaixo estão duas formas de Zanzibar, ambas UTM 37, mas uma é WGS84 e a outra é ARC1960:

insira a descrição da imagem aqui

Simbamangu
fonte
Tente, em vez de proj4string (tz.37.b) <- CRS ("+ init = epsg: 32737"), fazer como fiz no post pro4string (tz.37.b) <- proj4string (um shapefile com o projeto u são depois) e ver se funciona porque eu fiz como eu mencionei na descrição e ele apareceu na posição correta.
GSA
Portanto, se um shp não possui CRS, qual é o primeiro passo antes de aplicarmos o spTransform ?, dê os crs desejados e, em seguida, aplique o spTransform e tudo ficará bem?
GSA
Definitivamente, você deve descobrir o CRS correto e aplicá-lo primeiro, ou o spTransform não terá uma base para reprojetar o objeto!
Simbamangu 18/10/2015
4

Não: essas duas formas não são as mesmas, e isso ocorre por design. O motivo é que, em alguns casos, você deseja atribuir um CRS a um objeto que não possui um CRS ou alterar um CRS sem alterar as coordenadas e, em outros casos, deseja converter ou transformar coordenadas de um CRS em outro CRS. O primeiro (atribuindo, substituindo) é obtido por

proj4string(a)<-proj4string(s)

e o segundo (converter, transformar) por

a<-spTransform(a,CRS(proj4string(s))

Porque você não é o primeiro a pensar que o primeiro comando realmente faz o que o segundo comando, spemite o seguinte aviso quando você tenta o primeiro formulário, mas ajá possui um CRS diferente:

Warning message:
In ReplProj4string(obj, CRS(value)) :
  A new CRS was assigned to an object with an existing CRS:
+init=epsg:28992 +proj=sterea [etc]
without reprojecting.
For reprojection, use function spTransform in package rgdal

(que não é mais totalmente correta, pois spTransformé uma função spque chama métodos em rgdal).

Edzer Pebesma
fonte