Alternativas de spTransform?

8

Digamos que tenemos un archivo de forma con cierta proyección.

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

También tenemos los aeropuertos de España como puntos en otra proyección.

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

Si nuestro objetivo es colocar los puntos en el shapefile de España, tenemos que organizar las coordenadas para que sean las mismas. Por lo general, se hace así:

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

¿Pero es lo mismo con eso?

proj4string(a)<-proj4string(s)

En caso afirmativo, ¿por qué no es la forma estándar de hacerlo ya que es más simple y tiene que recurrir a spTransform?

gsa
fuente

Respuestas:

6

Asignar un CRS diferente no cambia la proyección de los datos espaciales subyacentes: el CRS es una parte interna del objeto espacial que le dice a R cómo interpretar las coordenadas espaciales.

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

Ahora transforme la forma en una zona UTM vecina:

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]

Las coordenadas de la geometría han cambiado, correctamente, para adaptarse a la zona UTM vecina. ¿Qué pasa si simplemente reasignamos el CRS sin transformación?

# 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

Nos han advertido ... entonces, ¿cómo son las 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]

Los números (coordenadas) son los originales de la zona UTM 36, pero la forma se asignará a la zona UTM incorrecta ahora y aparecerá en el lugar incorrecto.

Editar

Usando el método exacto sugerido por el 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]

La extensión es igual al objeto utm36 original, pero proj4string ahora es utm37. (Curiosamente, no hubo advertencia esta vez). Nota: el resultado es EXACTAMENTE el mismo que haber asignado el CRS con el código EPSG. Probar:

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

Veamos cómo se ven las formas reales, usando una capa UTM37 real de las regiones de Tanzania.

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

ingrese la descripción de la imagen aquí

No se alinean, por lo que este método no funciona. ¿Qué pasa con el objeto transformado ?

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

ingrese la descripción de la imagen aquí

El objeto transformado aparece en el lugar correcto (fondo verde oscuro).

Nota: si se está transformando entre zonas UTM pero el mismo dato (WGS84 aquí), las diferencias serán dramáticas, como en mi ejemplo, pero si se está transformando entre diferentes datos, las diferencias pueden ser mucho más sutiles. Por ejemplo, a continuación hay dos formas de Zanzíbar, ambas UTM 37, pero una es WGS84 y la otra es ARC1960:

ingrese la descripción de la imagen aquí

Simbamangu
fuente
Intente en lugar de proj4string (tz.37.b) <- CRS ("+ init = epsg: 32737") para hacer lo que hice en el post pro4string (tz.37.b) <- proj4string (un archivo de forma con el proj u son después) y ver si funciona porque lo hice como mencioné en la descripción y apareció en la posición correcta.
gsa
Entonces, si un shp no tiene CRS, ¿cuál es el primer paso antes de aplicar spTransform? ¿Darle los crs que queramos y luego aplicar spTransform y todo está bien?
gsa
Definitivamente, debe averiguar el CRS correcto y aplicarlo primero, ¡o spTransform no tendrá una base sobre la cual volver a proyectar el objeto!
Simbamangu
4

No: estas dos formas no son lo mismo, y esto es por diseño. La razón es que en algunos casos desea asignar un CRS a un objeto que no tiene un CRS o cambiar un CRS sin cambiar las coordenadas, y en otros casos desea convertir o transformar coordenadas de un CRS en otro CRS. El primero (asignar, reemplazar) se obtiene mediante

proj4string(a)<-proj4string(s)

y el segundo (convertir, transformar) por

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

Porque no eres el primero en pensar que el primer comando realmente hace lo que hace el segundo comando, te spda la siguiente advertencia cuando intentas el primer formulario pero aya tiene un 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 ya no es del todo correcto, ya que spTransformes una función spque llama a los métodos rgdal).

Edzer Pebesma
fuente