Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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? If so, how much -- centimeters, meters, kilometers?

Would the answer change if the second projection was in a different datum (e.g. NAD27)?

share|improve this question
    
Because of: transformation resolution, coordinate system resolution differences, and geometry storage resolution and tolerance. Each one of these "variables" are different. So you need to go read the documentation for each. – GISI Feb 10 at 1:55
    
...and if you are using ArcGIS there's potentially hundreds of transformation equations listed in order of highest resolution transformations for your data's spatial domain. – GISI Feb 10 at 1:59
1  
The usual result from A -> B -> A' is A ~= A', but adding datum transformation to the mix could really foul things up if you do it wrong. A great deal depends on how the coordinate references are defined (and therefore the truncation in the map units of each coordinate system). – Vince Feb 10 at 2:15
up vote 17 down vote accepted

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.

share|improve this answer
    
A great answer! – Deer Hunter Feb 10 at 3:27
    
This is even more thorough than I was hoping for, and very encouraging. Thank you for the test, and for the example script to replicate it with my own data! – Erica Feb 10 at 11:21
    
Can you include what you mean by the NAD83 UTM corners? If they're at the extremes of the zone (high latitude, for instance), using points within the US will likely give even better results. – mkennedy Feb 10 at 17:42
    
I assume, that the WGS84 bounds shipped with EPSG 26917 at spatialreference.org/ref/epsg/nad83-utm-zone-17n .. WGS84 Bounds: -84.0000, 24.0000, -78.0000, 83.0000 is the right region of interest. Do I made a mistake? – huckfinn Feb 10 at 17:56
    
@huckfinn Duh, I should have seen the values in the code! Sorry for the stupid question. Great values for a generalized answer about UTM. – mkennedy Feb 10 at 18:10

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.

share|improve this answer

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.

share|improve this answer
1  
Agree. Also, be aware ArcGIS displays 6 decimal places of a Double datatype by default in the Table View. You can change the properties of the field to display 12 decimal places, in the table view. Geographic xy values are typically 9 decimal places or so, double precision. – klewis Feb 10 at 18:34

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.