Konversi sistem koordinat geografis di R

14

Saya memiliki poin dalam sistem koordinat geografis dan saya ingin mengubahnya menjadi swiss grid (CH1903 +).

Contoh data:

 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

Hasil yang Dihormati:

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
convert coordinate-system r Topdombili
sumber
3
@ Harun, saya orang yang sama. Saya belum mendapatkan jawaban yang tepat di sana. Bisakah Anda membantu saya? Saya benar-benar kagum betapa Anda pilih-pilih.
Topdombili
1
@Top Ini bukan pilih-pilih, ini kebijakan StackExchange. Posting silang menciptakan semua jenis ketidakkonsistenan dan masalah. (Anda mungkin juga telah memperhatikan bahwa memposting di forum yang salah dapat memperoleh beberapa jawaban yang kurang berguna.) Silakan hapus versi SO yang Anda posting.
whuber
@ Topdombili, saya hanya ingin menunjukkan bahwa menurut jawaban whuber, nilai input ada di WGS84 dan sedang menjalani transformasi datum plus proyeksi ke kisi CH1903 + / LV95.
mkennedy
@kenkeny Terima kasih atas pengamatan itu. Saya lalai dalam tidak menjelaskan bahwa saya berasumsi koordinat asli (lat, lon) adalah WGS 84 (anggapan itu terkubur dalam komentar dalam kode). Jika tidak, ubah nilainya proj4string(d)sesuai. Perhatian saya terutama tertuju pada parameter easting dan northing yang salah, x0dan y0, karena beberapa referensi populer di Web (seperti dalam komentar pertama dalam kode) telah menjatuhkan angka mereka yang paling signifikan, sehingga menggusur semua titik dengan beberapa ribu kilometer :-).
whuber
1
@whuber, aduh re: referensi! Saya memang melihat komentar Anda tentang set input ke WGS 84, tetapi ingin memastikan Topdombili melihatnya.
mkennedy

Jawaban:

18

Gunakan paket RGDAL . Ada masalah CRS yang digunakan; RGDAL tidak mengenali kode EPSG. Anda harus memberikan parameter secara eksplisit, seperti yang ditunjukkan di sini. (Rupanya ini adalah perkiraan, tetapi mereka seharusnya cukup bagus. Mereka tampaknya berada dalam jarak sekitar 0,1 m dari nilai yang dimaksudkan.)

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

Plot

        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
sumber
Saya ingin bertanya hasil akurasi konversi mungkin lebih sedikit sementara ada tanpa nilai desimal sementara koordinat yang tersedia dalam hasil yang dihormati adalah 10 derajat terdekat, maksud saya 2 digit desimal.
Topdombili
1
Saya tidak mengerti komentar Anda. Apakah Anda bertanya tentang seberapa tepat koordinat yang diproyeksikan ketika nilai asli (lat, lon) ditentukan dengan presisi terbatas? Jika demikian, Anda mungkin menemukan jawaban ini bermanfaat.
whuber
1
rgdal mengenali EPSG: 2056, FWIW
mdsumner
@md, terima kasih! Saya telah menemukan referensi yang menyatakan ini adalah EPSG: 9814, tetapi RGDAL tidak mengenalinya.
whuber