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_Transformations_10.1.zip
(2) http://resources.arcgis.com/en/help/main/10.2/index.html#/Equation_based_methods/003r00000012000000/
(3) http://resources.arcgis.com/en/help/main/10.2/index.html#/Grid_based_methods/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