[R-sig-Geo] Change in projection results relative to 3 years ago
Roger Bivand
Roger.Bivand at nhh.no
Mon Oct 30 15:54:44 CET 2017
On Mon, 30 Oct 2017, MacQueen, Don wrote:
> 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
>
Briefly, please check the current contents of:
list.files("/Users/macqueen1/Library/R/3.4/library/rgdal/proj")
I get:
> locs.26743 <- spTransform(locs.ll, CRS("+init=epsg:26743"))
> coordinates(locs.ref)-coordinates(locs.26743)
coords.x1 coords.x2
[1,] 3.746539 -1.876668
[2,] 3.746607 -1.876466
with
> library(rgdal)
Loading required package: sp
rgdal: version: 1.2-15, (SVN revision 691)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
Path to GDAL shared files: /usr/local/share/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
Path to PROJ.4 shared files: (autodetected)
Linking to sp version: 1.2-5
(1.2-15 is the development version, but nothing changed in this respect
vis. 1.2-13)
and
> list.files("/usr/local/share/proj")
[1] "alaska" "CH" "conus"
[4] "epsg" "esri" "esri.extra"
[7] "FL" "GL27" "hawaii"
[10] "IGNF" "MD" "nad.lst"
[13] "nad27" "nad83" "ntf_r93.gsb"
[16] "ntv1_can.dat" "null" "nzgd2kgrid0005.gsb"
[19] "other.extra" "proj_def.dat" "prvi"
[22] "stgeorge" "stlrnc" "stpaul"
[25] "TN" "WI" "WO"
[28] "world"
Maybe you are being thrown back just to the ellipsoid if the NAD grids are
absent?
Hope this helps,
Roger
>
>
> 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
>
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
Editor-in-Chief of The R Journal, https://journal.r-project.org/index.html
http://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
More information about the R-sig-Geo
mailing list