[R-sig-Geo] Using spTransform() to reproduce another software package's transformation
Roger Bivand
Roger.Bivand at nhh.no
Fri Aug 29 09:26:25 CEST 2014
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.
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#//003r00000010000000:
>
> <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