it looks like the mysterious NAD_1927_To_NAD_1983_6 transform is actually a localized transform for the Quebec area. Furthermore, ESRI has informed us that this is not in any way a default transform, rather the first in the list. In summary, regardless of the software that you are using to do this type of work, always know your data and do your homework on the available methods. Thanks to Matt Wilkie for the detective work.
The Experiment
Quick comparison of the results of a vector projection using both OGR (GDAL wich uses the Proj4 library) and ArcGIS 9.0. Ideally both methods should return just about the same coordinates. The GDAL tools use the NADCON grid corection for datum shifts in the US and Canada, while ArcGIS provides several transformation methods (NADCON being one of them). Figure 1 demonstrates the dialog box where datum transformation parameters are set in ArcGIS. The default (first) method in this list NAD_1927_To_NAD_1983_6 is more than slightly confusing. Perhaps the rationale for this approach is explained somewhere, but I have certainly never come across it. In order to better understand any possible differences in the results of a datum shift operation, the three following operations were performed:
OGR ogr2ogr -s_srs 'prj=latlong datum=NAD83' -t_srs '+proj=utm +zone=11 +datum=NAD27' output.shp input.shp ArcToolbox transform from NAD83 GCS to NAD27 UTM Zone 11 (NADCON) ArcToolbox transform from NAD83 GCS to NAD27 UTM Zone 11 (NAD_1927_To_NAD_1983_6 transformation)
Analysis
Several comparisons where done to assess the three dataum shifting methods. First the output from the cs2cs command (used for datum shift operations among other things) was compared with the same operation on a the NGS NADCON website. The two methods appeard to be the same for at least 10 decimal places. Next, the first coordinate from a polygon projected across datums (with the three methods described above) was compared. The OGR and ArcGIS (NADCON) methods produced identical results, however the ArcGIS (NAD_1927_To_NAD_1983_6) method produced odd results. It appeared as if the coordinates where transformed, but not corectly. Finally, the same coordinate was datum shifted using proj directly: with the same odd results as the NAD_1927_To_NAD_1983_6 method. An analysis of these coordinates as compared to the corectly shifted coordinates is presented in figures 2 through 4.
Results
As discovered in an earlier experiment, vector projection within a single datum by OGR and ArcGIS yielded nearly identical results. This is certainly a good thing if you are using both of these programs in your work flow, or sharing with others critical of open source tools. Based on these simple tests, the output from OGR and ArcGIS (using the NADCON transform) are functionaly identical. This is also a good thing. However, the default transform presented in the project feature dialog within ArcToolbox (Figure 1) is still a bit puzzling. While it appears to perform some soft of coordinate projection, it is certainly not the correct one. Also, it is worth noting that using the proj command directly for a similar operation results in the same odd output! It should be said however, that proj is not designed for datum shifts. Instead, the command cs2cs should be used.
Conclusions
What have we learned? While I wasn't expecting OGR and ArcGIS to produce different results, it is nice to see that when used properly they produce nearly identical output. Perhaps someone can shed some light on the default, and somewhat bizzare NAD_1927_To_NAD_1983_6 datum transform method in ArcGIS 9.0. Details on the tests can be found below:
1. Compare datum shift (NAD83 -> NAD27) Proj 4.9 and NGS NADCON: [exactly the same]
echo "-119.560953035858461 37.06968077687506" | \ cs2cs +proj=latlong +datum=NAD83 +to +proj=latlong +datum=NAD27 -w5
NAD83 NAD27 Lat Lon Lat Lon cs2cs 37.06968077687506 -119.560953035858461 37d4'11.0339"N 119d33'35.93485"W NADCON 37.06968077687506 -119.560953035858461 37d4'11.03390"N 119d33'35.93485"W
2. Datum shift NAD83->NAD27 of first coordinate in sample polygon, into UTM Zone 11:
[note problems with ArcGIS default transform]
NAD83 NAD27 easting northing easting northing 2D shift ogr2ogr 272323.66132 4105670.61614 272404.60847 4105470.64705 215.731 ArcMap NADCON 272323.66132 4105670.61614 272404.60847 4105470.64715 215.731 ArcMap Default Transform 272323.66132 4105670.61614 272318.10172 4105467.33107 203.361
3. Quick test: does proj 4.9 return the same results as ogr2ogr across datums? [NO]
#1st coordinate from polygon GCS NAD83: echo "-119.560953035858461 37.06968077687506" | proj +proj=utm +zone=11 +datum=NAD83 272323.66 4105670.62 echo "-119.560953035858461 37.06968077687506" | proj +proj=utm +zone=11 +datum=NAD27 272318.10 4105467.33 203 m difference in UTM coordinates
4. Comparison of OGR (correct) and ArcGIS defualt (incorrect) datum shifted coordinates in R.
#analyse in R:
#cols 1, 2 are arcmap x,y
#cols 3, 4 are ogr x,y
paste arcmap_9.ascii ogr.ascii | awk ' /\./{print $1, $2, $3, $4}' > raw_coords.dat
#start R
R
a <- read.table('raw_coords.dat')
names(a) <- c('arc_x', 'arc_y', 'ogr_x', 'ogr_y')
attach(a)
#analysis of residuals
par(mfcol=c(2,2)) ; plot(lm(arc_x ~ ogr_x), col='blue', cex=0.3)
par(mfcol=c(2,2)) ; plot(lm(arc_y ~ ogr_y), col='blue', cex=0.3)
#RMSE:
sum(sqrt((ogr_x - arc_x)^2 + (ogr_y - arc_y)^2)) / length(ogr_x)
88.02054 [meters]
Fig 1: projection dialog in ArcMap
Fig 2: observed coordinate shift ArcMap default datum transform output in red.
Fig 3: ArcMap vs. OGR easting values
Fig 4: ArcMap vs. OGR northing values
