Czy przejście na nową projekcję, a następnie powrót, wpływa na dokładność danych?

13

Mam klasę obiektów (hrabstwa Karoliny Południowej, czyli dość duży obszar geograficzny) w samolocie stanowym NAD83 SC. Musi zostać przekształcony w drugą projekcję (NAD83 UTM 17), a następnie przekształcony z powrotem w oryginał. W tym celu użyję narzędzia Projekt Esri .

Czy ta podwójna transformacja może spowodować przesunięcie położenia współrzędnych wielokątów i o ile - centymetrów, metrów, kilometrów?

Erica
źródło
Ze względu na: rozdzielczość transformacji, różnice rozdzielczości układu współrzędnych oraz rozdzielczość i tolerancję przechowywania geometrii. Każda z tych „zmiennych” jest inna. Musisz więc przeczytać dokumentację dla każdego z nich.
GISI
... a jeśli korzystasz z ArcGIS, potencjalnie setki równań transformacji są wymienione w kolejności transformacji o najwyższej rozdzielczości dla domeny przestrzennej danych.
GISI
1
Zwykle wynik z A -> B -> A 'to A ~ = A', ale dodanie transformacji układu odniesienia do miksu może naprawdę popsuć sytuację, jeśli zrobisz to źle. Wiele zależy od tego, jak zdefiniowane są odniesienia do współrzędnych (a zatem obcięcia w jednostkach mapy każdego układu współrzędnych).
Vince

Odpowiedzi:

19

Nie wiem, którego silnika projekcyjnego używa ArcGis, ale bardzo interesujące pytanie dotyczy również proj.4. Dlatego próbuję przetestować silnik projekcyjny proj.4 w środowisku GNU-R. Używam narożników NAD 83 - UTM 17 i EPSG 26917 i rekursywnie ponownie 10000 i 1000000 razy i obliczam różnicę do wartości początkowych.

Oto wyniki:

Wygląda na to, że błąd „ponownej projekcji” mieści się w zakresie centymetrów dla 10000 pętli.

"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

I dojdzie do błędu w zakresie metrów, jeśli uruchomisz pętlę 1000000 razy.

"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

Oto skrypt.

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

Dalsze testy w środowisku statystyk powinny być łatwe. Skrypty i objaśnienia kodu dla środowiska linux są dostępne na stronie github.com/bigopensky .

huckfinn
źródło
Jest to nawet dokładniejsze niż się spodziewałem i bardzo zachęcające. Dziękuję za test i przykładowy skrypt do replikacji go z własnymi danymi!
Erica
Czy możesz podać, co masz na myśli, mówiąc o narożnikach NAD83 UTM? Jeśli znajdują się na krańcach strefy (na przykład na dużej szerokości geograficznej), użycie punktów w USA prawdopodobnie da jeszcze lepsze wyniki.
mkennedy
Zakładam, że granice WGS84 dostarczane z EPSG 26917 na spatialreference.org/ref/epsg/nad83-utm-zone-17n .. WGS84 Bounds: -84.0000, 24.0000, -78.0000, 83.0000są właściwym obszarem zainteresowania. Czy popełniłem błąd?
huckfinn
@huckfinn Duh, powinienem był zobaczyć wartości w kodzie! Przepraszam za głupie pytanie. Świetne wartości dla uogólnionej odpowiedzi na temat UTM.
mkennedy
7

Esri ma własny silnik projekcyjny.

Większość metod prognozowania i transformacji geograficznej / danych jest dobrze zachowana, gdy jest stosowana w odpowiednim obszarze zainteresowania. Jeśli znajdziesz się zbyt daleko poza strefą UTM, poprzeczny Mercator nie zawsze dokładnie „odwraca” (przelicz na długość i szerokość geograficzną) dokładnie. W rzutach używanych dla całego świata mogą występować pewne problemy na biegunach lub wokół biegunów lub południka +/- 180 lub „anty-południka” (południk, który znajduje się naprzeciwko środka rzutowanego układu odniesienia za pomocą współrzędnych).

Przebiegłem 4 punkty, które wypadają poza Karoliną Południową za pośrednictwem silnika projekcyjnego Esri. Aby wykonać test warunków skrajnych na 1k, 10k lub 1M punktów, będę musiał coś zakodować, ponieważ mój istniejący podobny test po prostu wykonuje „podróż w obie strony” - od projekcji geograficznej do projekcji. 32133 to NAD 1983 State Plane South Carolina (metry). 26917 to NAD 1983 UTM zone 17 North.

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

Widać więc, że mieliśmy dwa punkty, które wróciły o 10e-09.

Obsługa w ArcGIS komplikuje fakt, że istnieje odniesienie przestrzenne. Odniesienie przestrzenne obejmuje układ współrzędnych oraz pewne wartości przechowywania i analizy. Domyślnie układy współrzędnych wykorzystujące mierniki są przechowywane z dokładnością do jednej dziesiątej milimetra, 0,0001.

Ujawnienie: Pracuję dla Esri.

Mkennedy
źródło
5

Myślę, że jest to przypadek, w którym trzeba przetestować proponowany przepływ pracy pod kątem niektórych funkcji punktów testowych, do których można łatwo dodać pola współrzędnych XY.

Porównaj wartości XY swoich początkowych punktów z tymi, które rzutowałeś / przekształcałeś (jakkolwiek wiele razy), a oszacujesz różnicę.

PolyGeo
źródło
1
Zgodzić się. Należy również pamiętać, że ArcGIS wyświetla domyślnie 6 miejsc dziesiętnych typu danych Double w widoku tabeli. Możesz zmienić właściwości pola, aby wyświetlać 12 miejsc po przecinku w widoku tabeli. Geograficzne wartości xy wynoszą zwykle około 9 miejsc po przecinku, podwójna precyzja.
klewis