¿La transformación a una nueva proyección, luego la vuelta, 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. Debe transformarse en una segunda proyección (NAD83 UTM 17) y luego volver a la original. Estaré usando la herramienta de Esri para lograr esto.

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

    
pregunta Erica 10.02.2016 - 02:13

3 respuestas

19

No sé qué motor de proyección usa ArcGis, pero es una pregunta muy interesante también para proj.4. Así que le doy una oportunidad para probar el motor de proyección proj.4 dentro del entorno GNU-R. Utilizo las esquinas NAD 83 - UTM 17 y EPSG 26917 y las reproduzco 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 aumente a 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 script.

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

Las pruebas adicionales 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 .

    
respondido por el huckfinn 10.02.2016 - 03:53
7

Esri tiene su propio motor de proyección.

La mayoría de las proyecciones y los métodos de transformación de datos / geográficos se comportan bien cuando se utilizan en un área de interés apropiada. Si se aleja demasiado de una zona UTM, Mercator transversal no siempre es "inverso" (se convierte a 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 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 simplemente hace un 'viaje de ida y vuelta' - proyectado a 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

Para ver, tuvimos dos puntos que regresaron a las 10e-09.

El manejo en ArcGIS se complica 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: yo trabajo para Esri.

    
respondido por el mkennedy 10.02.2016 - 19:21
5

Creo que este es un caso en el que necesitas probar el 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 (muchas veces), y habrá cuantificado la diferencia.

    
respondido por el PolyGeo 10.02.2016 - 03:01

Lea otras preguntas en las etiquetas