Using spTransform() to reproduce another software package's transformation

Hi

It seems to me that you think that ESRI performs the "correct" 
transformation.

From http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//003r00000010000000:

<quote start>
Converting between NAD 1983 and WGS 1984

Originally, NAD 1983 and WGS 1984 were considered coincident. To minimize coordinate changes, NAD 1983 is tied to the North American and Pacific (for Hawaii, and so on) plates. WGS 1984 is tied to the International Terrestrial Reference System (ITRF), which is independent of the tectonic plates. Over time, the two coordinate systems have become increasingly different.

NAD_1983_To_WGS_1984_1: Published accuracy from EPSG is 2 meters. This transformation applies to the entire North American continent. This transformation uses the geocentric translation method, with the transformation's parameters (dx, dy, and dz) all equal to zeroes. This transformation treats the NAD 1983 and WGS 1984 datums as though they are equivalent.
NAD_1983_To_WGS_1984_2: Calculated by the U. S. Defence Mapping Agency (DMA), now known as the National Geospatial Intelligence Agency (NGA), for the Aleutian islands. Accuracy is listed by EPSG at +/-8 m.
NAD_1983_To_WGS_1984_3: Calculated by the NGA for Hawaii. Accuracy is listed by EPSG at +/-4 m.
NAD_1983_To_WGS_1984_4: Formerly applied within the 48 contiguous states, but is superseded by _5. This transformation method should no longer be used.
NAD_1983_To_WGS_1984_5: Transformation parameters calculated by the U.S. National Geodetic Survey (NGS) using CORS stations, and ties WGS 1984 to ITRF96. Accuracy according to EPSG is +/- 1 meter.
NAD_1983_To_WGS_1984_6, _7, and _8: Canadian NTv2 transformations, for the Quebec, Saskatchewan, and Alberta provinces, respectively.
<quote end>

So the precision of NAD_1983_To_WGS_1984_5 according to EPSG is +- 1 meter so if your results are in feet then they are within that limit.

In R the parameters for PROJ.4 are:

CRS("+init=epsg:26743")

CRS arguments:
+init=epsg:26743 +proj=lcc +lat_1=38.43333333333333
+lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5
+x_0=609601.2192024384 +y_0=0 +datum=NAD27 +units=us-ft +no_defs
+ellps=clrk66 +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat +towgs=-3.746315,1.876856

This corresponds to what you can find on www.epsg-registry.org

Is there a corresponding EPSG code for the ESRI NAD_1983_To_WGS_1984_5 transformation? Or can you by any other means find something similar to the CRS arguments? If so we can compare the CRS arguments from R and ESRI.




Yours sincerely / Med venlig hilsen


Frede Aakmann T?gersen
Specialist, M.Sc., Ph.D.
Plant Performance & Modeling

Technology & Service Solutions
T +45 9730 5135
M +45 2547 6050
frtog at vestas.com
http://www.vestas.com

Company reg. name: Vestas Wind Systems A/S
This e-mail is subject to our e-mail disclaimer statement.
Please refer to www.vestas.com/legal/notice
If you have received this e-mail in error please contact the sender.

-----Original Message-----
From: r-sig-geo-bounces at r-project.org [mailto:r-sig-geo-bounces at r-
project.org] On Behalf Of MacQueen, Don
Sent: 29. august 2014 01:58
To: r-sig-geo at r-project.org
Subject: [R-sig-Geo] Using spTransform() to reproduce another software
package's transformation

The program I work for has specified the use of a local coordinate reference
system and a method for transforming and projecting from WGS84 long/lat
to the local system. They use ESRI products to convert from long/lat to the
local system.

Since I do everything in R, naturally I wish to use spTransform() to replicate
their conversion. I've been using spTransform() for a number of years now,
and thought I understood what I've been doing.

But I have run into trouble. I would appreciate any advice.

I believe I have a reproducible example. Toward the end of this email are R
expressions (based on dput) that will create two SpatialPoints objects that
are used in the example. They need to be created first, before running the
example.

####
## before adding further detail and the example, here are some references
####

(1)
http://downloads2.esri.com/support/TechArticles/Geographic_Transformati
ons_10.1.zip
(2)
http://resources.arcgis.com/en/help/main/10.2/index.html#/Equation_base
d_methods/003r00000012000000/
(3)
http://resources.arcgis.com/en/help/main/10.2/index.html#/Grid_based_m
ethods/003r00000013000000/



The programs's specified CRS is epsg 26743 = California State Plane Zone 3
NAD27 US feet (out of my control!)

The specified method for transforming and projecting from WGS84 long/lat
to the local CRS consists of two steps:
 1) transform and project to epsg 2227 = California State Plane Zone 3 NAD83
US feet
 2) transform to epsg 26743 = California State Plane Zone 3 NAD27 US feet

When doing the steps in the ESRI software's projection tool:
 step 1) use what ESRI calls "NAD_1983_To_WGS_1984_5"  (wkid 1515 in
reference 1)
 step 2) use what ESRI calls "NAD_1927_To_NAD_1983_NADCON"  (wkid 1241
in reference 1)

According to reference 1 "NAD_1983_To_WGS_1984_5" is a "coordinate
frame" transformation.
Based on reference 2, this means it is a 7 parameter Bursa-Wolf method
Also based on reference 1, "NAD_1927_To_NAD_1983_NADCON" is a grid-
based method


As I hope you will see, a naive use of spTransform() produces coordinates
that differ from the ESRI two-step process by approximately 3.7 ft (x) and -
1.9 ft (y). This is too large for our use case. I also believe as a matter of
principle that it should be possible to do better (I'd like to believe that any
transformation possible in ESRI is also possible using PROJ.4).


##
## reproducible example begins
##

## define two example points in WGS84 long/lat
locs.xy <- cbind(
             c(-121.524764291826,-121.523480804667),
             c(37.6600366036405,37.6543604613483)
             )

locs.ll <- SpatialPoints(locs.xy, proj4string=CRS("+proj=longlat
+datum=WGS84") )

## source the expressions near the bottom of this email to create
##    locs.ref
##    locs.step1.esri

## use spTransform to go from WGS84 to the local system in one step:
locs.26743 <- spTransform(locs.ll, CRS("+init=epsg:26743"))

## not close enough:
coordinates(locs.ref)-coordinates(locs.26743)
##      coords.x1 coords.x2
## [1,]  3.746539 -1.876668
## [2,]  3.746607 -1.876466


## spTransform equivalent of ESRI step 1
locs.step1.proj4 <- spTransform(locs.ll, CRS("+init=epsg:2227"))

## not close enough, essentially the same difference as above
coordinates(locs.step1.esri)-coordinates(locs.step1.proj4)
##      coords.x1 coords.x2
## [1,]  3.746244 -1.877057
## [2,]  3.746315 -1.876856


## next, try the spTransform equivalent of ESRI step 1, but specifying the
seven parameters
## note, had to reverse the sign of the rotation args from wkid 1515 in
reference 1;
## evidently the PROJ.4 default is the "position vector" method (reference 2)

crs.step1.cf <- CRS('+proj=lcc +lat_1=38.43333333333333
+lat_2=37.06666666666667 +lat_0=36.5 +lon_0=-120.5\
 +x_0=2000000.0 +y_0=500000.0 +ellps=GRS80 +datum=NAD83\
 +units=us-ft +no_defs\
 +towgs84=-0.991,1.9072,0.5129,0.025789908,0.0096501,0.0116599,0.0')

## by the way, this alternative to specifying the CRS gives the same result
##  crs.step1.cf <- CRS("+init=epsg:2227 +towgs84=-
0.991,1.9072,0.5129,0.025789908,0.0096501,0.0116599,0.0")

locs.step1.cf <- spTransform(locs.ll, crs.step1.cf)

## good enough (hooray!)
coordinates(locs.step1.esri)-coordinates(locs.step1.cf)
##          coords.x1     coords.x2
## [1,] -3.469177e-06 -5.122274e-08
## [2,] -3.418885e-06 -7.380731e-08



## now try for step 2 using spTranform()
locs.step2.cf <- spTransform(locs.step1.cf, CRS("+init=epsg:26743"))

## the original difference is back!
coordinates(locs.ref)-coordinates(locs.step2.cf)
##      coords.x1 coords.x2
## [1,]  3.746539 -1.876668
## [2,]  3.746608 -1.876466

## the implication is that in doing the transformation to epsg 26743, it
reversed the effect of the 7-parameter method

## attempt to prevent that:

locs.tmp <- locs.step1.cf
proj4string(locs.tmp) <- CRS("+init=epsg:2227")
## Warning message:
## In `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) :
##   A new CRS was assigned to an object with an existing CRS:
## +proj=lcc +lat_1=38.43333333333333 +lat_2=37.06666666666667
+lat_0=36.5 +lon_0=-120.5 +x_0=2000000.0 +y_0=500000.0 +ellps=GRS80
+datum=NAD83 +units=us-ft +no_defs +towgs84=-
0.991,1.9072,0.5129,0.025789908,0.0096501,0.0116599,0.0
## without reprojecting.
## For reprojection, use function spTransform in package rgdal

locs.step2.cfb <- spTransform(locs.tmp, CRS("+init=epsg:26743"))

## This actually works; the difference is now acceptably small
coordinates(locs.ref)-coordinates(locs.step2.cfb)
##         coords.x1    coords.x2
## [1,] 0.0003266879 0.0006651825
## [2,] 0.0003261503 0.0006651356


## Another way, perhaps more appropriate (and with no warning message)
is recreate the SpatialPoints
## object instead of modifying it:
locs.tmp <- SpatialPoints(coordinates(locs.step1.cf),
proj4string=CRS("+init=epsg:2227"))
locs.step2.cfb <- spTransform(locs.tmp, CRS("+init=epsg:26743"))

coordinates(locs.ref)-coordinates(locs.step2.cfb)
##         coords.x1    coords.x2
## [1,] 0.0003266879 0.0006651825
## [2,] 0.0003261503 0.0006651356



This suggests to me that in PROJ.4 the specifications for the CRS are in some
way mixed in with the specifications for how the geographic transformation
is performed. Apparently, one cannot specify the transformation method
independent of the CRS specification. To put it another way, I have two ways
of converting from the original CRS to the first step's target CRS. The target
CRS is the same either way. But a second conversion to another CRS is
affected by the method of the first one.

Have I interpreted correctly? If so, I guess it doesn't seem appropriate-the
conversion from one CRS to another should depend only on what the CRS is,
not on how it got there.

Is there something I don't understand so that this kind of dependency is
appropriate?

In the end, I guess I do have a solution, but I kind of don't like it. I have to
insert a "correction" to the CRS. Is there a better way?




#####
### source the following to create objects used by the reproducible example
above
#####

## the two points converted from long/lat using the complete ESRI "two-
step process"
## saved as a shapefile, loaded into R using readOGR(). Then just the
coordinates
## were "dput"
locs.ref <- new(
              "SpatialPoints",
              coords = structure(c(1703671.30566227, 1704020.20113366,
                424014.398045834, 421943.708664294), .Dim = c(2L, 2L),
                .Dimnames = list(NULL, c("coords.x1", "coords.x2")))
              , bbox = structure(
                  c(1703671.30566227, 421943.708664294,
                    1704020.20113366, 424014.398045834),
                  .Dim = c(2L, 2L),
                  .Dimnames = list(c("coords.x1",  "coords.x2"), c("min", "max")))
              , proj4string =
              new("CRS",
                  projargs = "+proj=lcc +lat_1=37.06666666666667
+lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5
+x_0=609601.2192024384 +y_0=0 +datum=NAD27 +units=us-ft +no_defs
+ellps=clrk66 +nadgrids=@conus, at alaska, at ntv2_0.gsb, at ntv1_can.dat"
                  )
              )


## the points converted to epsg 2227 using ESRI's step 1
## saved as a shapefile, loaded into R using readOGR
## Then just the coordinates were "dput"
locs.step1.esri <- new(
                     "SpatialPoints",
                     coords = structure(c(6265039.1378244, 6265388.04257557,
                       2064418.92932968, 2062348.22239488), .Dim = c(2L, 2L),
                       .Dimnames = list(NULL, c("coords.x1", "coords.x2")))
                     , bbox = structure(
                         c(6265039.1378244, 2062348.22239488,
                           6265388.04257557, 2064418.92932968),
                         .Dim = c(2L, 2L),
                         .Dimnames = list(c("coords.x1", "coords.x2"), c("min", "max")))
                     , proj4string = new("CRS",
                         projargs = "+proj=lcc +lat_1=37.06666666666667
+lat_2=38.43333333333333 +lat_0=36.5 +lon_0=-120.5 +x_0=2000000
+y_0=500000.0000000001 +datum=NAD83 +units=us-ft +no_defs
+ellps=GRS80 +towgs84=0,0,0"
                         )
                     )


#####
###package and session information
#####


Loading required package: sp
Checking rgeos availability: TRUE

Loading required package: rgdal
rgdal: version: 0.8-16, (SVN revision 498)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.11.0, released 2014/04/16
Path to GDAL shared files:
/Users/macqueen1/Library/R/3.1/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files:
/Users/macqueen1/Library/R/3.1/library/rgdal/proj

sessionInfo()

R version 3.1.1 (2014-07-10)
Platform: x86_64-apple-darwin13.1.0 (64-bit)

locale:
[1] C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] sp_1.0-15       rgdal_0.8-16    maptools_0.8-30 xlsx_0.5.5
[5] xlsxjars_0.6.0  rJava_0.9-6     rmacq_1.3-1

loaded via a namespace (and not attached):
[1] foreign_0.8-61  grid_3.1.1      lattice_0.20-29 tools_3.1.1

----
two final comments:

I understand that this is really a PROJ.4 question, so I hope that R-sig-geo
folks don't mind too much; apologies in advance if so.

I hope my email software truly sends a plain text email like it claims. But I'm
not sure I trust it!