Does transforming to a new projection, then back, affect data accuracy?

Question

I have a feature class (South Carolina counties, so a fairly large geographic area) in NAD83 SC State Plane. It needs to be transformed to a second projection (NAD83 UTM 17), then transformed back to the original. I'll be using Esri's Project tool to accomplish this.

Can this dual transformation cause a shift in the location of the polygons coordinates, and by how much -- centimeters, meters, kilometers?

Answer 1

I don't know which projection engine ArcGis uses, but a very interesting question also for proj.4. So I give it a try to test the proj.4 projection engine within the GNU-R environment. I use the NAD 83 - UTM 17 corners and EPSG 26917 and reproject it 10000 and 1000000 times recursivly and calculate the difference to the start values.

Here are the results:

It seems that the "reprojection" error in is within a centimeter range for 10000 loops.

"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

And grow to an error in a meter range if you run the loop 1000000 times.

"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

Here is the 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)

Further tests within a statistcs environment should be easy. The scripts and code explanation for a linux environment are available at github.com/bigopensky.

Answer 2

Esri has its own projection engine.

Most projections and geographic/datum transformations methods are well-behaved when used in an appropriate area of interest. If you get too far outside of a UTM zone, transverse Mercator doesn't always 'inverse' (convert to latitude-longitude) exactly. Projections used for the whole world may have some issues at or around the poles or the +/-180 meridian or the 'anti-meridian' (the meridian that's opposite the center of the projected coordinate reference system).

I ran 4 points that fall outside South Carolina through the Esri projection engine. For a stress test of 1k or 10k or 1M points, I'll have to code something as my existing similar test just does a 'round-trip'--projected to geographic to projected. 32133 is NAD 1983 State Plane South Carolina (meters). 26917 is 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

So you can see we had two points that came back at 10e-09.

The handling in ArcGIS is complicated by the fact that there's a spatial reference. The spatial reference includes the coordinate system plus some storage and analysis values. By default, coordinate systems that use meters are stored with a precision of a tenth of a millimeter, 0.0001.

Disclosure: I work for Esri.

Answer 3

I think this is a case where you need to test your proposed workflow against some test point features, which are easy to add XY coordinate fields to.

Compare the XY values of your initial points with those that you have projected/transformed (however many times), and you will have quantified the difference.