[R-sig-Geo] Using spTransform() to reproduce another software package's transformation
MacQueen, Don
macqueen1 at llnl.gov
Wed Sep 3 19:02:05 CEST 2014
Roger, Hermann, Frede,
My thanks to all who took a look at this. Your responses, including
Roger¹s re-expression using cs2cs on the proj list and the discussion
there (subject: NAD_1983_To_WGS_1984_5), are very helpful.
The proj.4 parameters that Roger has provided do better than I had been
able to achieve, such that I am completely comfortable telling my
co-workers that I can reproduce their coordinate transformations. And at a
programmatic level, that is what I really needed, so thanks!
-Don
--
Don MacQueen
Lawrence Livermore National Laboratory
7000 East Ave., L-627
Livermore, CA 94550
925-423-1062
On 8/30/14, 12:00 PM, "Roger Bivand" <Roger.Bivand at nhh.no> wrote:
>On Fri, 29 Aug 2014, Roger Bivand wrote:
>
>> On Fri, 29 Aug 2014, Frede Aakmann Tøgersen wrote:
>>
>>> Hi
>>>
>>> It seems to me that you think that ESRI performs the "correct"
>>> transformation.
>>
>> Thanks, Frede. My guess is that this is a question for the proj-devel
>>list,
>> probably using cs2cs as the workhorse. Please re-cast your example to
>>use
>> cs2cs from the console prompt - I think that will be the most efficient
>> resolution. I can post on the proj list if you like.
>
>Having posted on the proj list, and after being helped by Hermann Peifer,
>this now does better:
>
>crs.step1.cf <- CRS(paste("+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 +units=us-ft +no_defs",
>"+towgs84=-0.991,1.9072,0.5129,0.025789908,0.0096501,0.0116599,0.0"))
>crs.step1.cf
>locs.step1.cf <- spTransform(locs.ll, crs.step1.cf)
>coordinates(locs.step1.esri)-coordinates(locs.step1.cf)
># coords.x1 coords.x2
>#[1,] -3.469177e-06 -5.122274e-08
>#[2,] -3.418885e-06 -7.846393e-08
>
>
>locs.tmp <- locs.step1.cf
>suppressWarnings(proj4string(locs.tmp) <- CRS(paste("+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 +units=us-ft",
>"+no_defs +nadgrids=@null")))
>locs.step2.cfb <- spTransform(locs.tmp, CRS("+init=epsg:26743"))
>coordinates(locs.ref)-coordinates(locs.step2.cfb)
># coords.x1 coords.x2
>#[1,] -1.028087e-05 -5.030306e-07
>#[2,] -1.081964e-05 -5.360343e-07
>
>This with (forthcoming) PROJ.4 4.9.0 and (forthcoming) rgdal 0.9-1, not
>checked with earlier versions, but should work. The sources are not very
>forthcoming on why the +nadgrids=@null trick works, but it appears to
>force the omission a step in datum transformation, which itself is a
>multi-step process internally.
>
>Hope this helps,
>
>Roger
>
>>
>> I suspect that the +/- 1m is what is available, but if it was possible
>>to
>> help in the PROJ.4 framework, all other downstream software components
>>would
>> benefit.
>>
>> I do see that Frank Warmerdam posted:
>>
>> http://lists.maptools.org/pipermail/proj/2008-September/003833.html
>>
>> 6 years ago in answer to a question about NAD_1983_To_WGS_1984_5.
>>
>> Hope this helps,
>>
>> Roger
>>
>>>
>>>> From
>>>>
>>>>http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//003r0000
>>>>0010000000:
>>>
>>> <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!
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>
>>> _______________________________________________
>>> 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; fax +47 55 95 91 00
>e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list