[R-sig-Geo] Change in projection results relative to 3 years ago

Mon Oct 30 15:38:32 CET 2017

Three years ago Roger and Hermann Peifer graciously helped me solve a projection issue in R. Something has since changed; results are different now.

I'd once again appreciate assistance... and am more or less hoping that someone will have some idea about what has changed, hence what I can do to return to my previously blissful state.

Thanks in advance,
-Don

My original request is here
   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021611.html
and Roger's solution is here
   https://stat.ethz.ch/pipermail/r-sig-geo/2014-August/021623.html
(actually, they did more than help, the pretty much did all the work!)

I recently discovered the coordinates produced by that solution are now
roughly 100m different from what they were then.
(I've done quite a bit of checking to make sure I'm not making some
dumb mistake, and feel confident I haven't. Hopefully I'm right!)

A small set of example locations is transformed from long/lat to a local coordinate system using
  (A) a reference method deemed correct
  (B) an R method using sp::spTransform()
The goal was to have B reproduce A. Three years ago, it did. Now it does not.

(The accuracy of method A is not the issue. It uses standard methods implemented in ESRI software.)

The reproducible example below (same one as 3 years ago) does not use the full version of method B, but it does illustrate how something has changed.

##
## reproducible example begins
##

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

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

## Outside of R, reference method A converts from long/lat
## to a local system, EPSG 26743, which is:
##   California State Plane Zone III, NAD27, feet;
##   http://spatialreference.org/ref/epsg/26743/
## using an ESRI "two-step process" (some detail at the end),
## saved as a shapefile, loaded into R using readOGR().

## Here is the "dput" of those coordinates from three years ago (the example location transformed by method A)
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"
        )
)

## locs.ref and locs.ll are created as above, the same as they were 3 years ago
## (same reproducible example data as then)

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

## distances relative to the reference transformation in August 2014
coordinates(locs.ref)-coordinates(locs.26743)
##      coords.x1 coords.x2
## [1,]  3.746539 -1.876668
## [2,]  3.746607 -1.876466

## distances relative to the reference transformation now
coordinates(locs.ref)-coordinates(locs.26743)
##      coords.x1 coords.x2
## [1,]  309.9325  20.05890
## [2,]  309.9136  20.03086

##
## a transformation that formerly was within a few feet of the example location is now some 300 feet away
##

## the ESRI "two step" starts with the lon/lat,
## Step 1 converts transforms them using what ESRI calls
##     NAD_1983_To_WGS_1984_5       (wkid 1515)
## Step 2 transforms the resulting coordinates using
##     NAD_1927_To_NAD1983_NADCON   (wkid 1241)

--------  Session info now:

library(rgdal)
Loading required package: rgdal
Loading required package: sp
rgdal: version: 1.2-13, (SVN revision 686)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
Path to GDAL shared files: /Users/macqueen1/Library/R/3.4/library/rgdal/gdal
Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
Path to PROJ.4 shared files: /Users/macqueen1/Library/R/3.4/library/rgdal/proj
Linking to sp version: 1.2-5

sessionInfo()
R version 3.4.2 (2017-09-28)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] C

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

other attached packages:
[1] openxlsx_4.0.17 Rkml_1.6-5      mymaps_1.4-2    rmacq_1.3-5     rgdal_1.2-13    sp_1.2-5      

loaded via a namespace (and not attached):
[1] compiler_3.4.2  tools_3.4.2     Rcpp_0.12.13    grid_3.4.2      lattice_0.20-35

------ Session info then:
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    

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509

--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
Lab cell 925-724-7509