¿La transformación a una nueva proyección, luego de regreso, afecta la precisión de los datos?

13

Tengo una clase de entidad (condados de Carolina del Sur, por lo que es un área geográfica bastante grande) en NAD83 SC State Plane. Necesita ser transformado a una segunda proyección (NAD83 UTM 17), luego transformado nuevamente al original. Usaré la herramienta Proyecto de Esri para lograr esto.

¿Puede esta doble transformación causar un cambio en la ubicación de las coordenadas de los polígonos y en cuánto - centímetros, metros, kilómetros?

arcgis-desktop coordinate-system united-states Erica
fuente
Debido a: resolución de transformación, diferencias de resolución del sistema de coordenadas y resolución y tolerancia de almacenamiento de geometría. Cada una de estas "variables" son diferentes. Por lo tanto, debe leer la documentación de cada uno.
GISI
... y si está utilizando ArcGIS, hay potencialmente cientos de ecuaciones de transformación enumeradas en el orden de las transformaciones de mayor resolución para el dominio espacial de sus datos.
GISI
1
El resultado habitual de A -> B -> A 'es A ~ = A', pero agregar una transformación de datos a la mezcla realmente podría estropear las cosas si lo haces mal. Mucho depende de cómo se definen las referencias de coordenadas (y, por lo tanto, el truncamiento en las unidades de mapa de cada sistema de coordenadas).
Vince

Respuestas:

19

No sé qué motor de proyección usa ArcGis, pero una pregunta muy interesante también para el proyecto 4. Así que intento probar el motor de proyección proj.4 dentro del entorno GNU-R. Utilizo las esquinas NAD 83 - UTM 17 y EPSG 26917 y lo vuelvo a proyectar 10000 y 1000000 veces de forma recursiva y calculo la diferencia con los valores iniciales.

Aquí están los resultados:

Parece que el error de "reproyección" está dentro de un rango de centímetros para 10000 bucles.

"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

Y crezca hasta un error en un rango de medidor si ejecuta el bucle 1000000 veces.

"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

Aquí está el guión.

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

Otras pruebas dentro de un entorno estadístico deberían ser fáciles. Los scripts y la explicación del código para un entorno Linux están disponibles en github.com/bigopensky .

huckfinn
fuente
Esto es aún más exhaustivo de lo que esperaba y muy alentador. ¡Gracias por la prueba y por el script de ejemplo para replicarlo con mis propios datos!
Erica
¿Puedes incluir lo que quieres decir con las esquinas NAD83 UTM? Si están en los extremos de la zona (alta latitud, por ejemplo), el uso de puntos dentro de los EE. UU. Probablemente dará mejores resultados.
mkennedy
Supongo que los límites de WGS84 enviados con EPSG 26917 en spatialreference.org/ref/epsg/nad83-utm-zone-17n .. WGS84 Bounds: -84.0000, 24.0000, -78.0000, 83.0000es la región correcta de interés. ¿Cometí un error?
huckfinn
@huckfinn Duh, ¡debería haber visto los valores en el código! Perdón por la estúpida pregunta. Excelentes valores para una respuesta generalizada sobre UTM.
mkennedy
7

Esri tiene su propio motor de proyección.

La mayoría de las proyecciones y los métodos de transformación geográfica / datum se comportan bien cuando se usan en un área de interés apropiada. Si se aleja demasiado de una zona UTM, el Mercator transversal no siempre es 'inverso' (se convierte en latitud-longitud) exactamente. Las proyecciones utilizadas para todo el mundo pueden tener algunos problemas en o alrededor de los polos o el meridiano +/- 180 o el 'anti-meridiano' (el meridiano que está opuesto al centro del sistema de referencia de coordenadas proyectadas).

Corrí 4 puntos que caen fuera de Carolina del Sur a través del motor de proyección Esri. Para una prueba de esfuerzo de 1k o 10k o 1M puntos, tendré que codificar algo ya que mi prueba similar existente solo hace un 'viaje de ida y vuelta', proyectado de geográfico a proyectado. 32133 es NAD 1983 State Plane South Carolina (metros). 26917 es NAD 1983 UTM zona 17 norte.

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

Como pueden ver, tuvimos dos puntos que volvieron a las 10e-09.

El manejo en ArcGIS es complicado por el hecho de que hay una referencia espacial. La referencia espacial incluye el sistema de coordenadas más algunos valores de almacenamiento y análisis. Por defecto, los sistemas de coordenadas que usan medidores se almacenan con una precisión de una décima de milímetro, 0.0001.

Divulgación: trabajo para Esri.

mkennedy
fuente
5

Creo que este es un caso en el que necesita probar su flujo de trabajo propuesto contra algunas características de puntos de prueba, a las que es fácil agregar campos de coordenadas XY.

Compare los valores XY de sus puntos iniciales con los que ha proyectado / transformado (sin importar cuántas veces), y habrá cuantificado la diferencia.

PolyGeo
fuente
1
De acuerdo. Además, tenga en cuenta que ArcGIS muestra 6 lugares decimales de un tipo de datos Doble de forma predeterminada en la Vista de tabla. Puede cambiar las propiedades del campo para mostrar 12 decimales, en la vista de tabla. Los valores geográficos xy son típicamente 9 lugares decimales más o menos, doble precisión.
klewis