Beeinträchtigt die Umwandlung in eine neue Projektion die Datengenauigkeit?

13

Ich habe eine Feature-Class (Grafschaften in South Carolina, also ein ziemlich großes geografisches Gebiet) in NAD83 SC State Plane. Es muss in eine zweite Projektion umgewandelt werden (NAD83 UTM 17) und dann wieder in das Original umgewandelt werden. Ich werde das Projekt-Tool von Esri verwenden , um dies zu erreichen.

Kann diese doppelte Transformation zu einer Verschiebung der Position der Polygonkoordinaten und um wie viel Zentimeter, Meter, Kilometer führen?

Erica
quelle
Aus folgenden Gründen: Transformationsauflösung, Unterschiede in der Koordinatensystemauflösung sowie Auflösung und Toleranz des Geometriespeichers. Jede dieser "Variablen" ist anders. Sie müssen also die Dokumentation für jeden lesen.
GISI
... und wenn Sie ArcGIS verwenden, werden möglicherweise Hunderte von Transformationsgleichungen in der Reihenfolge der Transformationen mit der höchsten Auflösung für die räumliche Domäne Ihrer Daten aufgelistet.
GISI
1
Das übliche Ergebnis von A -> B -> A 'ist A ~ = A', aber die Hinzufügung einer Datumstransformation zum Mix kann die Sache wirklich durcheinander bringen, wenn Sie es falsch machen. Viel hängt davon ab, wie die Koordinatenreferenzen definiert sind (und daher die Kürzung in den Karteneinheiten jedes Koordinatensystems).
Vince

Antworten:

19

Ich weiß nicht, welche Projektions-Engine ArcGis verwendet, aber eine sehr interessante Frage auch für Proj.4. Deshalb versuche ich, die Projektions-Engine proj.4 in der GNU-R-Umgebung zu testen. Ich benutze die Ecken NAD 83 - UTM 17 und EPSG 26917 und projiziere sie 10000 und 1000000 mal rekursiv und berechne die Differenz zu den Startwerten.

Hier sind die Ergebnisse:

Es scheint, dass der "Reprojektions" -Fehler bei 10000 Schleifen innerhalb eines Zentimeterbereichs liegt.

"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

Wenn Sie die Schleife 1000000-mal ausführen, kann dies zu einem Fehler in einem Messbereich führen.

"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

Hier ist das Drehbuch.

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

Weitere Tests in einer statistischen Umgebung sollten einfach sein. Die Skripte und Code-Erklärungen für eine Linux-Umgebung finden Sie unter github.com/bigopensky .

huckfinn
quelle
Dies ist noch gründlicher als ich gehofft hatte und sehr ermutigend. Vielen Dank für den Test und für das Beispielskript, um es mit meinen eigenen Daten zu replizieren!
Erica
Können Sie einschließen, was Sie mit den NAD83 UTM-Ecken meinen? Wenn sie sich am äußersten Rand der Zone befinden (z. B. im hohen Breitengrad), wird die Verwendung von Punkten in den USA wahrscheinlich zu noch besseren Ergebnissen führen.
Mkennedy
Ich gehe davon aus, dass die mit EPSG 26917 ausgelieferten WGS84-Grenzen unter spaciousreference.org/ref/epsg/nad83-utm-zone-17n .. WGS84 Bounds: -84.0000, 24.0000, -78.0000, 83.0000die richtige Region von Interesse sind. Habe ich einen Fehler gemacht?
Huckfinn
@huckfinn Duh, ich hätte die Werte im Code sehen sollen! Entschuldigung für die blöde Frage. Tolle Werte für eine allgemeine Antwort zu UTM.
Mkennedy
7

Esri hat eine eigene Projektionsmaschine.

Die meisten Projektionen und Methoden zur geografischen / Datumsumwandlung verhalten sich in einem geeigneten Bereich von Interesse gut. Wenn Sie sich zu weit außerhalb einer UTM-Zone befinden, wird Mercator in Querrichtung nicht immer genau umgekehrt (in Breiten- und Längengrad konvertiert). Projektionen, die für die ganze Welt verwendet werden, können Probleme an oder um die Pole oder am +/- 180-Meridian oder am 'Anti-Meridian' (dem Meridian, der sich gegenüber dem Zentrum des projizierten Koordinatenreferenzsystems befindet) haben.

Ich lief 4 Punkte, die außerhalb von South Carolina durch die Esri-Projektionsmaschine fallen. Für einen Stresstest von 1k oder 10k oder 1M Punkten muss ich etwas codieren, da mein bestehender ähnlicher Test nur einen "Roundtrip" durchführt - von geografisch nach projiziert. 32133 ist NAD 1983 State Plane South Carolina (Meter). 26917 ist NAD 1983 UTM Zone 17 Nord.

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

Sie sehen also, wir hatten zwei Punkte, die zwischen 10 und 09 zurückkamen.

Die Handhabung in ArcGIS wird durch einen räumlichen Bezug erschwert. Der Raumbezug enthält das Koordinatensystem sowie einige Speicher- und Analysewerte. Standardmäßig werden Koordinatensysteme, die Meter verwenden, mit einer Genauigkeit von einem Zehntel Millimeter (0,0001) gespeichert.

Offenlegung: Ich arbeite für Esri.

mkennedy
quelle
5

Ich denke, dies ist ein Fall, in dem Sie Ihren vorgeschlagenen Workflow mit einigen Testpunkt-Features testen müssen, denen Sie leicht XY-Koordinatenfelder hinzufügen können.

Vergleichen Sie die XY-Werte Ihrer Anfangspunkte mit denen, die Sie projiziert / transformiert haben (wie oft auch immer), und Sie werden den Unterschied quantifiziert haben.

PolyGeo
quelle
1
Zustimmen. Beachten Sie außerdem, dass ArcGIS in der Tabellenansicht standardmäßig 6 Dezimalstellen eines Double-Datentyps anzeigt. Sie können die Eigenschaften des Felds so ändern, dass in der Tabellenansicht 12 Dezimalstellen angezeigt werden. Geografische xy-Werte sind in der Regel 9 Dezimalstellen oder so, doppelte Genauigkeit.
Klewis