1 votos

¿Por qué ha fallado spTransform aquí?

Estoy tratando de convertir un mapa que encontré aquí en lat/lon compatible con otro conjunto de datos que creé a partir de las lat/lons de Google Maps.

Desde el directorio con ese .shp archivo, corro:

library(rgdal)
SG = readOGR('../data', 'MP14_PLNG_AREA_NO_SEA_PL')

Un subconjunto de los datos que estoy tratando de co-trazar es:

dim = c(5L, 2L)
# polygon coordinate matrices
poly_mats = list(structure(c(103.843, 103.843, 103.854, 103.854, 103.843, 
                             1.318, 1.324, 1.324, 1.318, 1.318), .Dim = dim), 
                 structure(c(103.843, 103.843, 103.854, 103.854, 103.843,
                             1.324, 1.329, 1.329, 1.324, 1.324), .Dim = dim), 
                 structure(c(103.854, 103.854, 103.865, 103.865, 103.854,
                             1.318, 1.324, 1.324, 1.318, 1.318), .Dim = dim),
                 structure(c(103.854, 103.854, 103.865, 103.865, 103.854, 
                             1.324, 1.329, 1.329, 1.324, 1.324), .Dim = dim))
# nest properly to create SP object
polySP = SpatialPolygons(lapply(seq_along(poly_mats), function(ii) {
  Polygons(list(Polygon(poly_mats[[ii]])), ID = ii)
  # assign lat/lon CRS
}), proj4string = CRS('+init=epsg:3857'))

Ambos objetos se ven bien cuando se trazan por separado:

png('separate.png')
par(mfrow = c(1, 2))
plot(SG)
plot(polySP)
dev.off()

polygons separately

No es que deba importar, pero lo cierto es que los objetos están actualmente en diferentes SRI:

identical(proj4string(SG), proj4string(polySP))
# [1] FALSE

Así que deberíamos ser capaces de transformar SG así:

SG = spTransform(SG, proj4string(polySP))

Pero esto falla:

png('together_no.png')
plot(SG)
plot(polySP, col = 'red', add = TRUE)
dev.off()

not together

Parece que spTransform ha fallado, ya que las coordenadas en la salida no tienen sentido como lat/lons:

coordinates(SG)[1:5, ]
#       [,1]     [,2]
# 0 11559649 153646.0
# 1 11569258 147405.4
# 2 11559462 150848.3
# 3 11544079 146374.9
# 4 11549925 150925.5

¿Qué está fallando/cómo se puede rectificar?

4voto

Jay Bazuzi Puntos 194

Has creado tus polígonos con epsg:3857 :

proj4string = CRS('+init=epsg:3857'))

que es el Web Mercator de Google, no el lat-long. Quiere decir que epsg:4326 , WGS84 lat-long, muy probablemente.

Una lectura más detallada:

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X