Przekształcanie układu współrzędnych geograficznych w R.

Mam punkty w układzie współrzędnych geograficznych i chciałem przekonwertować je na szwajcarską siatkę (CH1903 +).

Przykładowe dane:

 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

Szacowane wyniki:

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
źródło

@Aaron, jestem tym samym facetem. Nie mam tam właściwej odpowiedzi. Czy możesz mi pomóc Jestem naprawdę zaskoczony, jak bardzo jesteś wybredna.

— Topdombili

@Top To nie jest wybredność, to polityka StackExchange. Cross-posting stwarza wszelkiego rodzaju niespójności i problemy. (Być może zauważyłeś również, że publikowanie na niewłaściwym forum może również dać mniej niż przydatne odpowiedzi). Usuń opublikowaną wersję SO.

— whuber

@Tomdombili, chciałem tylko zauważyć, że zgodnie z odpowiedzią Whubera wartości wejściowe są na WGS84 i przechodzą transformację odniesienia plus rzut na siatkę CH1903 + / LV95.

— mkennedy

@mkennedy Dziękuję za tę obserwację. Byłem niedbały w nie wyjaśniając, że mam założyć oryginał (szer) współrzędnych WGS 84 (to założenie jest pochowany w komentarzu w kodzie). Jeśli nie są, zmień odpowiednio wartość proj4string(d). Moją uwagę zwróciły przede wszystkim fałszywe parametry wschodu i północy, x0a y0ponieważ niektóre popularne odniesienia w sieci (takie jak w pierwszym komentarzu w kodzie) porzuciły swoje najbardziej znaczące cyfry, tym samym przesuwając wszystkie punkty o kilka tysięcy kilometrów :-).

— whuber

@ whuber, ouch re: referencje! Widziałem twój komentarz na temat zestawu danych wejściowych do WGS 84, ale chciałem się upewnić, że Topdombili go zobaczy.

— mkennedy

Użyj pakietu RGDAL . Jest problem, którego CRS użyć; RGDAL nie rozpoznaje kodu EPSG. Musisz podać parametry jawnie, jak pokazano tutaj. (Najwyraźniej są to przybliżenia, ale powinny być całkiem dobre. Wydaje się, że mieszczą się w granicach 0,1 m zamierzonych wartości.)

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

Działki

        lon        lat  
[1] 2579408,24 1079721,15
[2] 2579333,69 1079729,18
[3] 2579263,65 1079770,55
[4] 2579928,04 1080028,46
[5] 2579763,65 1079868,92
[6] 2579698,00 1079767,97
[7] 2579543,10 1079640,65
[8] 2580023,55 1080049,26
[9 ,] 2579859.11 1079875.27

— Whuber
źródło

Chciałem zapytać, czy wynik dokładności konwersji może być mniejszy, gdy są bez wartości dziesiętnych, podczas gdy dostępne współrzędne w szacowanych wynikach są z dokładnością do 10 stopni, miałem na myśli 2 cyfry dziesiętne.

— Topdombili

Nie rozumiem twojego komentarza. Czy pytasz o dokładność rzutowanych współrzędnych, gdy oryginalne wartości (szerokość, długość) są określone z ograniczoną precyzją? Jeśli tak, możesz uznać tę odpowiedź za przydatną.

— whuber

rgdal rozpoznaje EPSG: 2056, FWIW

— mdsumner

@md Dziękujemy! Znalazłem odniesienie stwierdzające, że jest to EPSG: 9814, ale RGDAL tego nie rozpoznał.

— whuber