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)
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)
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:
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
y el segundo (convertir, transformar) por
Porque no eres el primero en pensar que el primer comando realmente hace lo que hace el segundo comando, te
sp
da la siguiente advertencia cuando intentas el primer formulario peroa
ya tiene un CRS diferente:(que ya no es del todo correcto, ya que
spTransform
es una funciónsp
que llama a los métodosrgdal
).fuente