La transformation en une nouvelle projection, puis en arrière, affecte-t-elle la précision des données?

J'ai une classe d'entités (comtés de Caroline du Sud, donc une zone géographique assez étendue) dans NAD83 SC State Plane. Il doit être transformé en une seconde projection (NAD83 UTM 17), puis retransformé à l'original. Je vais utiliser l' outil Projet d'Esri pour y parvenir.

Cette double transformation peut-elle entraîner un changement dans l'emplacement des coordonnées des polygones et de combien - centimètres, mètres, kilomètres?

Je ne sais pas quel moteur de projection ArcGis utilise, mais une question très intéressante aussi pour proj.4. J'essaie donc de tester le moteur de projection proj.4 dans l'environnement GNU-R. J'utilise les coins NAD 83 - UTM 17 et EPSG 26917 et le reprojete 10000 et 1000000 fois de manière récursive et calcule la différence avec les valeurs de départ.

Voici les résultats:

Il semble que l'erreur de "reprojection" se situe dans une plage d'un centimètre pour 10000 boucles.

"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

Et passez à une erreur dans une plage de mètres si vous exécutez la boucle 1000000 fois.

"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

Voici le 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)

D'autres tests dans un environnement statistique devraient être faciles. Les scripts et l'explication du code pour un environnement linux sont disponibles sur github.com/bigopensky .

Esri a son propre moteur de projection.

La plupart des projections et des méthodes de transformations géographiques / datum se comportent bien lorsqu'elles sont utilisées dans une zone d'intérêt appropriée. Si vous vous éloignez trop loin d'une zone UTM, le Mercator transverse n'inverse pas toujours exactement (conversion en latitude-longitude). Les projections utilisées pour le monde entier peuvent avoir des problèmes au niveau ou autour des pôles ou du méridien +/- 180 ou de l'anti-méridien (le méridien qui est opposé au centre du système de référence de coordonnées projeté).

J'ai couru 4 points qui tombent en dehors de la Caroline du Sud grâce au moteur de projection Esri. Pour un test de stress de 1k ou 10k ou 1M points, je vais devoir coder quelque chose car mon test similaire existant ne fait qu'un `` aller-retour '' - projeté à géographique à projeté. 32133 est NAD 1983 State Plane South Carolina (mètres). 26917 est 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

Vous pouvez donc voir que nous avions deux points qui sont revenus à 10e-09.

La gestion dans ArcGIS est compliquée par le fait qu'il existe une référence spatiale. La référence spatiale comprend le système de coordonnées ainsi que certaines valeurs de stockage et d'analyse. Par défaut, les systèmes de coordonnées qui utilisent des mètres sont stockés avec une précision d'un dixième de millimètre, 0,0001.

Divulgation: je travaille pour Esri.

Je pense que c'est un cas où vous devez tester votre flux de travail proposé par rapport à certaines fonctionnalités de point de test, auxquelles il est facile d'ajouter des champs de coordonnées XY.

Comparez les valeurs XY de vos points initiaux avec celles que vous avez projetées / transformées (cependant plusieurs fois), et vous aurez quantifié la différence.