Conversión del sistema de coordenadas geográficas en R

14

Tengo puntos en el sistema de coordenadas geográficas y quería convertirlos a cuadrícula suiza (CH1903 +).

Data de muestra:

 id       lon      lat
  2 7.173500 45.86880
  3 7.172540 45.86887
  4 7.171636 45.86924
  5 7.180180 45.87158
  6 7.178070 45.87014
  7 7.177229 45.86923  
  8 7.175240 45.86808
  9 7.181409 45.87177
  10 7.179299 45.87020

Resultados Respetados:

id       E              N      
2     2579408.2431  1079721.1499
3     2579333.7158  1079729.1852
4     2579263.6502  1079770.1125
5     2579928.0358  1080028.4605
6     2579763.6471  1079868.9218
7     2579698.0756  1079767.9762
8     2579543.1019  1079640.6512
9     2580023.6226  1080049.2672
10    2579859.1889  1079875.2740 
Topdombili
fuente
3
@ Aaron, soy el mismo chico. No tengo la respuesta adecuada allí. ¿Puedes ayudarme? Estoy realmente sorprendido de lo exigente que eres.
Topdombili
1
@ Arriba Esto no es exigente, es la política de StackExchange La publicación cruzada crea todo tipo de inconsistencias y problemas. (También puede haber notado que publicar en el foro incorrecto también puede obtener respuestas poco útiles). Elimine la versión SO que publicó.
whuber
@Topdombili, solo quería señalar que, según la respuesta de Whuber, los valores de entrada están en WGS84 y están experimentando una transformación de datos más una proyección a la cuadrícula CH1903 + / LV95.
mkennedy
@mkennedy Gracias por esa observación. Fui negligente al no explicar que asumí que las coordenadas originales (lat, lon) son WGS 84 (esa suposición está oculta en un comentario en el código). Si no es así, cambie el valor de la proj4string(d)consecuencia. Mi atención se dirigió principalmente a los falsos parámetros este y norte, x0y y0, debido a que algunas referencias populares en la Web (como en el primer comentario en el código) han bajado sus dígitos más significativos, desplazando todos los puntos unos pocos miles de kilómetros. :-).
whuber
1
@whuber, ouch re: las referencias! Vi su comentario sobre la entrada establecida en WGS 84, pero quería asegurarme de que Topdombili lo viera.
mkennedy

Respuestas:

18

Use el paquete RGDAL . Hay un problema de qué CRS usar; RGDAL no reconoce el código EPSG. Debe proporcionar los parámetros explícitamente, como se muestra aquí. (Aparentemente, estas son aproximaciones, pero se supone que son bastante buenas. Parecen estar dentro de aproximadamente 0.1 m de los valores previstos).

# References:
# http://lists.maptools.org/pipermail/proj/2001-September/000248.html (has typos)
# http://www.remotesensing.org/geotiff/proj_list/swiss_oblique_cylindrical.html
#
# Input coordinates.
#
x <- c(7.173500, 7.172540, 7.171636, 7.180180, 7.178070, 7.177229, 7.175240, 7.181409, 7.179299)
y <- c(45.86880, 45.86887, 45.86924, 45.87158, 45.87014, 45.86923, 45.86808, 45.87177, 45.87020)
#
# Define the coordinate systems.
#
library(rgdal)
d <- data.frame(lon=x, lat=y)
coordinates(d) <- c("lon", "lat")
proj4string(d) <- CRS("+init=epsg:4326") # WGS 84
CRS.new <- CRS("+proj=somerc +lat_0=46.9524056 +lon_0=7.43958333 +ellps=bessel +x_0=2600000 +y_0=1200000 +towgs84=674.374,15.056,405.346 +units=m +k_0=1 +no_defs")
# (@mdsumner points out that
#    CRS.new <- CRS("+init=epsg:2056")
# will work, and indeed it does. See http://spatialreference.org/ref/epsg/2056/proj4/.)
d.ch1903 <- spTransform(d, CRS.new)
#
# Plot the results.
#
par(mfrow=c(1,3))
plot.default(x,y, main="Raw data", cex.axis=.95)
plot(d, axes=TRUE, main="Original lat-lon", cex.axis=.95)
plot(d.ch1903, axes=TRUE, main="Projected", cex.axis=.95)
unclass(d.ch1903)

Parcelas

        lon        lat  

[1,] 2.579.408,24 1.079.721,15
[2,] 2.579.333,69 1.079.729,18
[3,] 2.579.263,65 1.079.770,55
[4,] 2.579.928,04 1.080.028,46
[5,] 2.579.763,65 1.079.868,92
[6,] 2.579.698,00 1.079.767,97
[7,] 2.579.543,10 1.079.640,65
[8,] 2.580.023,55 1.080.049,26
[9 ,] 2579859.11 1079875.27

whuber
fuente
Quería preguntarle que el resultado de la precisión de conversión podría ser menor mientras no haya valores decimales, mientras que las coordenadas disponibles en los resultados respetados están en 10 grados más cercanos, es decir, 2 dígitos decimales.
Topdombili
1
No entiendo tu comentario. ¿Pregunta qué tan precisas son las coordenadas proyectadas cuando los valores originales (lat, lon) se especifican con precisión limitada? Si es así, puede encontrar esta respuesta útil.
whuber
1
rgdal reconoce EPSG: 2056, FWIW
mdsumner
@md Gracias! Encontré una referencia que dice que esto es EPSG: 9814, pero RGDAL no lo reconoció.
whuber