Apakah mengubah ke proyeksi baru, lalu kembali, memengaruhi akurasi data?

13

Saya memiliki kelas fitur (kabupaten Carolina Selatan, jadi area geografis yang cukup besar) di NAD83 SC State Plane. Itu perlu ditransformasikan ke proyeksi kedua (NAD83 UTM 17), kemudian diubah kembali ke aslinya. Saya akan menggunakan alat Proyek Esri untuk mencapai ini.

Dapatkah transformasi ganda ini menyebabkan pergeseran lokasi koordinat poligon, dan seberapa jauh - sentimeter, meter, kilometer?

Erica
sumber
Karena: resolusi transformasi, perbedaan resolusi sistem koordinat, dan resolusi dan toleransi penyimpanan geometri. Masing-masing dari "variabel" ini berbeda. Jadi, Anda perlu membaca dokumentasi untuk masing-masing.
GISI
... dan jika Anda menggunakan ArcGIS, ada ratusan persamaan transformasi yang berpotensi dalam urutan transformasi resolusi tertinggi untuk domain spasial data Anda.
GISI
1
Hasil yang biasa dari A -> B -> A 'adalah A ~ = A', tetapi menambahkan transformasi datum ke dalam campuran dapat benar-benar merusak segalanya jika Anda melakukannya dengan salah. Banyak tergantung pada bagaimana referensi koordinat didefinisikan (dan oleh karena itu pemotongan dalam unit peta dari setiap sistem koordinat).
Vince

Jawaban:

19

Saya tidak tahu mesin proyeksi mana yang digunakan ArcGis, tetapi pertanyaan yang sangat menarik juga untuk proj.4. Jadi saya mencobanya untuk menguji mesin proyeksi proj.4 dalam lingkungan GNU-R. Saya menggunakan sudut NAD 83 - UTM 17 dan EPSG 26917 dan memproyeksikan ulang 10000 dan 1000000 kali secara berulang dan menghitung perbedaannya dengan nilai awal.

Inilah hasilnya:

Tampaknya kesalahan "proyeksi ulang" berada dalam kisaran sentimeter untuk 10.000 putaran.

"LON/LAT differences after  10000  loops"
           DLON          DLAT
1 -2.441464e-07 -1.341807e-07
2  2.441129e-07 -1.341807e-07
3  1.852679e-07 -1.691737e-08
4 -1.853157e-07 -1.691819e-08

"X/Y differences after  10000  loops"
            DX           DY
1 -0.025169783 -0.014338141
2  0.025166375 -0.014338208
3  0.002419045 -0.002016762
4 -0.002419690 -0.002016889

Dan berkembang menjadi kesalahan dalam rentang meter jika Anda menjalankan loop 1000000 kali.

"LON/LAT differences after  1000000  loops"
           DLON          DLAT
1 -2.441464e-05 -1.341845e-05
2  2.441128e-05 -1.341846e-05
3  1.852621e-05 -1.691837e-06
4 -1.853105e-05 -1.691828e-06

"X/Y differences after  1000000  loops"
          DX         DY
1 -2.5172288 -1.4339977
2  2.5168869 -1.4340064
3  0.2419201 -0.2017070
4 -0.2419859 -0.2017094

Ini skripnya.

# load the package
require('proj4')

# the LON/LAT frame of NAD83 UTM 17 
lon = c(-84.00, -78.00, -84.00, -78.00 ) 
lat = c( 24.00,  24.00,  83.00,  83.00)

# build the projection conform object
ll0 = matrix(c(lon,lat),nrow=4,ncol=2)
xy0 = project(ll0,"+init=epsg:26917",ellps.default='GRS80')

# make a copy
ll1 = ll0
xy1 = xy0

# number of iterations
num = 1000000

# reproject the stuff num times
for(i in 1:num) {
 # project forward  
 xy1 = project(ll1,"+init=epsg:26917", ellps.default='GRS80')
 # project backward
 ll1 = project(xy1,"+init=epsg:26917", inverse=T, ellps.default='GRS80')
}

# build difference table ll
dll = as.data.frame(ll1-ll0)
names(dll) = c('DLON','DLAT')
# print results LON/LAT
print(paste("LON/LAT differences after ", num," loops"))
print(dll)

# build difference table xy
dxy = as.data.frame(xy1-xy0)
names(dxy) = c('DX','DY')
# print results X/Y
print(paste("X/Y differences after ", num," loops"))
print(dxy)

Tes lebih lanjut dalam lingkungan statistik harus mudah. Penjelasan skrip dan kode untuk lingkungan linux tersedia di github.com/bigopensky .

huckfinn
sumber
Ini bahkan lebih menyeluruh daripada yang saya harapkan, dan sangat membesarkan hati. Terima kasih untuk tes ini, dan untuk skrip contoh untuk mereplikasi dengan data saya sendiri!
Erica
Bisakah Anda memasukkan apa yang Anda maksud dengan sudut UTM NAD83? Jika mereka berada di ujung zona (misalnya, lintang tinggi), menggunakan titik di AS kemungkinan akan memberikan hasil yang lebih baik.
mkennedy
Saya berasumsi, bahwa batas WGS84 yang dikirim dengan EPSG 26917 di spatialreference.org/ref/epsg/nad83-utm-zone-17n .. WGS84 Bounds: -84.0000, 24.0000, -78.0000, 83.0000adalah wilayah yang tepat untuk diminati. Apakah saya melakukan kesalahan?
huckfinn
@ Huckfinn Duh, saya seharusnya melihat nilai-nilai dalam kode! Maaf atas pertanyaan bodoh ini. Nilai hebat untuk jawaban umum tentang UTM.
mkennedy
7

Esri memiliki mesin proyeksi sendiri.

Sebagian besar metode proyeksi dan transformasi geografis / datum berperilaku baik ketika digunakan dalam bidang yang sesuai minat. Jika Anda terlalu jauh di luar zona UTM, Mercator melintang tidak selalu 'terbalik' (konversikan ke lintang-bujur). Proyeksi yang digunakan untuk seluruh dunia mungkin memiliki beberapa masalah di atau sekitar kutub atau +/- 180 meridian atau 'anti-meridian' (meridian yang berlawanan dengan pusat sistem referensi koordinat yang diproyeksikan).

Saya berlari 4 poin yang jatuh di luar Carolina Selatan melalui mesin proyeksi Esri. Untuk stress test 1k atau 10k atau 1M poin, saya harus memberi kode sesuatu karena tes serupa saya yang ada hanya melakukan 'round-trip' - diproyeksikan ke geografis untuk diproyeksikan. 32133 adalah NAD 1983 State Plane South Carolina (meter). 26917 adalah NAD 1983 UTM zona 17 Utara.

C:\Users\melita>inverse 32133
382000 20000
      -83.40806392522212        31.98974518135408
382000 383000
      -83.50098893136905        35.26180827475587
839100 20000
      -78.57184097446545        31.98934439195045
839100 383000
      -78.47814111839074        35.26139222680582

C:\Users\melita>forward 26917
  -83.40806392522212        31.98974518135408
       272490.5730967618        3541832.738731374
  -83.50098893136905        35.26180827475587
       272485.6257057797         3904944.98998655
  -78.57184097446545        31.98934439195045
       729409.4734382738        3541830.781689366
  -78.47814111839074        35.26139222680582
       729414.4926270114        3904946.919009762

C:\Users\melita>inverse 26917
 272490.5730967618        3541832.738731374
      -83.40806392522212        31.98974518135408
  272485.6257057797         3904944.98998655
      -83.50098893136905        35.26180827475587
  729409.4734382738        3541830.781689366
      -78.57184097446545        31.98934439195045
  729414.4926270114        3904946.919009762
      -78.47814111839074        35.26139222680582
^Z

C:\Users\melita>forward 32133
  -83.40806392522212        31.98974518135408
                382000.0                  20000.0
  -83.50098893136905        35.26180827475587
                382000.0                 383000.0
  -78.57184097446545        31.98934439195045
                839100.0        19999.99999999814
  -78.47814111839074        35.26139222680582
                839100.0        382999.9999999981

Jadi Anda dapat melihat kami memiliki dua poin yang kembali pada 10e-09.

Penanganan di ArcGIS dipersulit oleh fakta bahwa ada referensi spasial. Referensi spasial mencakup sistem koordinat ditambah beberapa nilai penyimpanan dan analisis. Secara default, sistem koordinat yang menggunakan meter disimpan dengan presisi sepersepuluh milimeter, 0,0001.

Pengungkapan: Saya bekerja untuk Esri.

mkennedy
sumber
5

Saya pikir ini adalah kasus di mana Anda perlu menguji alur kerja yang diusulkan terhadap beberapa fitur titik uji, yang mudah untuk menambahkan bidang koordinat XY.

Bandingkan nilai-nilai XY poin awal Anda dengan yang telah Anda proyeksikan / ubah (betapapun berkali-kali), dan Anda akan menghitung perbedaannya.

PolyGeo
sumber
1
Setuju. Perlu diketahui juga bahwa ArcGIS menampilkan 6 tempat desimal dari Double datatype secara default di Table View. Anda bisa mengubah properti bidang untuk menampilkan 12 tempat desimal, dalam tampilan tabel. Nilai xy geografis biasanya 9 angka desimal atau lebih, presisi ganda.
klewis