[R-sig-Geo] Using spTransform() to reproduce another software package's transformation
Roger Bivand
Roger.Bivand at nhh.no
Sat Aug 30 21:00:31 CEST 2014
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#//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