Change in projection results relative to 3 years ago

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Change in projection results relative to 3 years ago

MacQueen, Don
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
 
 

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Change in projection results relative to 3 years ago

Roger Bivand
Administrator
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
> [hidden email]
> 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: [hidden email]
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

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: Change in projection results relative to 3 years ago

MacQueen, Don
I get

> list.files("/Users/macqueen1/Library/R/3.4/library/rgdal/proj")
 [1] "CH"           "GL27"         "IGNF"         "epsg"         "esri"        
 [6] "esri.extra"   "nad.lst"      "nad27"        "nad83"        "other.extra"
[11] "proj_def.dat" "world"      

Considerably different indeed, and conus in particular is missing.

I have a somewhat vague memory of some r-sig-geo traffic some months ago that might be relevant.

Thanks
-Don

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

On 10/30/17, 7:55 AM, "Roger Bivand" <[hidden email]> wrote:

    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
    > [hidden email]
    > 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: [hidden email]
    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
   

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Change in projection results relative to 3 years ago

MacQueen, Don
In reply to this post by Roger Bivand
I should have checked this before my last email:

Since I have the KyngChaos frameworks:

[69]% ls /Library/Frameworks/PROJ.framework/unix/share/proj/
./                      WI                      nad.lst                 proj_def.dat
../                     WO                      nad27                   prvi
CH                      alaska                  nad83                   stgeorge
FL                      conus                   ntf_r93.gsb             stlrnc
GL27                    epsg                    ntv1_can.dat            stpaul
IGNF                    esri                    null                    world
MD                      esri.extra              nzgd2kgrid0005.gsb
TN                      hawaii                  other.extra

If my rgdal was installed from the CRAN binary, which I believe it was, then perhaps if I install from source it will pick up those files.

-Don

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

On 10/30/17, 7:55 AM, "Roger Bivand" <[hidden email]> wrote:

    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
    > [hidden email]
    > 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: [hidden email]
    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
   

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Change in projection results relative to 3 years ago

Roger Bivand
Administrator
Fill in from http://download.osgeo.org/proj/proj-datumgrid-1.6.zip, this was temporarily missing in the CRAN OSX binary.

Roger

Roger Bivand
Norwegian School of Economics
Bergen, Norway




On Mon, Oct 30, 2017 at 4:46 PM +0100, "MacQueen, Don" <[hidden email]<mailto:[hidden email]>> wrote:


I should have checked this before my last email:

Since I have the KyngChaos frameworks:

[69]% ls /Library/Frameworks/PROJ.framework/unix/share/proj/
./                      WI                      nad.lst                 proj_def.dat
../                     WO                      nad27                   prvi
CH                      alaska                  nad83                   stgeorge
FL                      conus                   ntf_r93.gsb             stlrnc
GL27                    epsg                    ntv1_can.dat            stpaul
IGNF                    esri                    null                    world
MD                      esri.extra              nzgd2kgrid0005.gsb
TN                      hawaii                  other.extra

If my rgdal was installed from the CRAN binary, which I believe it was, then perhaps if I install from source it will pick up those files.

-Don

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



On 10/30/17, 7:55 AM, "Roger Bivand"  wrote:

    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
    > [hidden email]
    > 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: [hidden email]
    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




        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand
Department of Economics
Norwegian School of Economics
Helleveien 30
N-5045 Bergen, Norway
Reply | Threaded
Open this post in threaded view
|

Re: Change in projection results relative to 3 years ago

MacQueen, Don
That solved the problem, thanks!

-Don

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



From: Roger Bivand <[hidden email]>
Date: Monday, October 30, 2017 at 8:54 AM
To: MacQueen_Don <[hidden email]>
Cc: "[hidden email]" <[hidden email]>
Subject: Re: [R-sig-Geo] Change in projection results relative to 3 years ago

Fill in from http://download.osgeo.org/proj/proj-datumgrid-1.6.zip, this was temporarily missing in the CRAN OSX binary.
Roger
Roger Bivand
Norwegian School of Economics
Bergen, Norway



On Mon, Oct 30, 2017 at 4:46 PM +0100, "MacQueen, Don" <[hidden email]<mailto:[hidden email]>> wrote:

I should have checked this before my last email:



Since I have the KyngChaos frameworks:



[69]% ls /Library/Frameworks/PROJ.framework/unix/share/proj/

./                      WI                      nad.lst                 proj_def.dat

../                     WO                      nad27                   prvi

CH                      alaska                  nad83                   stgeorge

FL                      conus                   ntf_r93.gsb             stlrnc

GL27                    epsg                    ntv1_can.dat            stpaul

IGNF                    esri                    null                    world

MD                      esri.extra              nzgd2kgrid0005.gsb

TN                      hawaii                  other.extra



If my rgdal was installed from the CRAN binary, which I believe it was, then perhaps if I install from source it will pick up those files.



-Don



--

Don MacQueen

Lawrence Livermore National Laboratory

7000 East Ave., L-627

Livermore, CA 94550

925-423-1062

Lab cell 925-724-7509







On 10/30/17, 7:55 AM, "Roger Bivand"  wrote:



    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

    > [hidden email]

    > 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: [hidden email]

    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





        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo