[R-sig-Geo] spTransform: wrong results??

[Roger Bivand]

On Fri, 10 Jul 2009, Agustin Lobo wrote:

> But if I use:
> meuse.lonlatintlWGS84  <- spTransform(meuse.utm, CRS("+proj=lonlat
> +datum=WGS84"))
>
> get the same result. Also, don't know where I should put the towgs84=
> field, have searched the R doc for
> towgs and do not find anything. I've looked in the doc you mention,
> and towgs84= must include a set of
> parameters that cannot find. The doc states
> "For instance, the following demonstrates converting from the Greek GGRS87
> datum to WGS84.<p>
>
> <pre>
> % cs2cs +proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62 \
>    +to +proj=latlong +datum=WGS84
> "
> but where can I find those 3 parameters for my case? Also, why is not
> spTransform() taking care
> of that? The spTransform() doc does not mention that towgs84 is
> required, actually there are several
> examples using the meuse dataset and trasforming to different datums:
>
> proj4string(meuse) <- CRS("+proj=stere +lat_0=52.15616055555555
> +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000
> +ellps=bessel +units=m")
> meuse.utm <- spTransform(meuse, CRS("+proj=utm +zone=32 +ellps=WGS84"))
>
> and the towgs84 is not there. Does this mean that coordinates in
> meuse.utm remain in ellpsoid Bessel instead of WGS84
> as specified in spTransform() ?
>
> On the other hand, I've made:
> EPSG <- make_EPSG()
> and found
> SG[27,3]
> [1] "+proj=longlat +ellps=WGS84 +no_defs"
>> EPSG[2446,]
>      code                  note
> 2446 23031 # ED50 / UTM zone 31N
>                                                 prj4
> 2446 +proj=utm +zone=31 +ellps=intl +units=m +no_defs
> and
>> EPSG[27,3]
> [1] "+proj=longlat +ellps=WGS84 +no_defs"
>
> no +datum field.
>
> This is so confusing...

Nobody can sort it out for anyone. EPSG don't give +towgs84 unless the 
authority is unquestioned. Until we had satellites, most surveying was 
done to local or national datum values, and ED50 pre-dates satellites. 
WGS84 is the first authoritative world geodetic system, hence its name, 
and is the same one that GPS satellites use. Because cartography is much 
older, and often military (secret), the conversion parameters are not well 
established or tabulated.

The way of doing things at the moment is to try to find an unequivocal 
authority (Grids & Datums is as good as they come), and use that to 
supplement EPSG - oil companies need to know where their oilfields and 
wells are. Once you think you are there, match with GPS or GE in 
geographical coordinates.

Roger

>
> Agus
>
> 2009/7/10 Roger Bivand <Roger.Bivand at nhh.no>:
>> On Fri, 10 Jul 2009, Agustin Lobo wrote:
>>
>>> I'm finding conflicting results at comparing
>>> results of converting UTM ED50 coordinates
>>> to lon,lat WGS84 using spTransform() (thus, proj4)
>>> and other programs such as GlobalMapper and Compegps.
>>>
>>> My main surprise is that spTransform() yields the same result
>>> whatever the ellps field for the inverse transform.
>>
>> Wrong tree to bark up, I'm afraid. The tags you need and are missing are
>> +datum= and +towgs84= (from memory, you'll see from the EPSG tags). The
>> +ellps= tag sets the ellipsoid, but doesn't change the default +datum= from
>> WGS84. Note that ED50 is a collection of datum values, and in general EPSG
>> does not give 3 or 7 parameter datum transformation values unless they are
>> known with certainty. On the Proj.4 side, you'll find some documentation on:
>>
>> http://svn.osgeo.org/metacrs/proj/trunk/proj/html/gen_parms.html
>>
>> otherwise use makeEPSG() in rgdal to search the bundled EPSG database, and
>> do check the "Grids & Datums" essays on:
>>
>> http://www.asprs.org/resources/grids/
>>
>> for further help in tracking down possible +towgs84= parameter values.
>>
>> Repeating, unless you give a +datum=, it will be assumed to be WGS84, so you
>> shouldn't be surprised when you get similar results.
>>
>> Sorry to disappoint, but surveying and cartography have their own
>> traditions, some of which are localised as well.
>>
>> Hope this helps,
>>
>> Roger
>>
>>> This is
>>> what I'm doing with the meuse dataset as an example, you can
>>> see that the results of the inverse transform to lon,lat WGS84 are
>>> identical to
>>> those to lon,lat ED50, both starting from UTM50 or from UTMWGS84
>>>
>>> At the end
>>> I include results with GlobalMapper for the first point. Results
>>> with GlobalMapper do differ for lon,lat WGS84 and lon,lat ED50
>>>
>>> Any help would be appreciated, as I'm doing all my work
>>> in UTM ED50 but have to pass the coordinates as lon,lat WGS84
>>> for the flight.
>>>
>>> data(meuse)
>>> coordinates(meuse) <- c("x", "y")
>>> proj4string(meuse) <- CRS("+proj=stere +lat_0=52.15616055555555
>>> +lon_0=5.38763888888889 +k=0.999908 +x_0=155000 +y_0=463000 +ellps=bessel
>>> +units=m")
>>> meuse.utm <- spTransform(meuse, CRS("+proj=utm +zone=32 +ellps=WGS84"))
>>> coordinates(meuse.utm)[1:3,]
>>>           x       y
>>> [1,] 272571.9 5653978
>>> [2,] 272522.4 5653927
>>> [3,] 272661.2 5653899
>>>
>>> meuse.utmintl <- spTransform(meuse, CRS("+proj=utm +zone=32 +ellps=intl"))
>>> coordinates(meuse.utmintl)[1:3,]
>>>           x       y
>>> [1,] 272561.0 5654094
>>> [2,] 272511.5 5654043
>>> [3,] 272650.3 5654015
>>>
>>> meuse.lonlatWGS84WGS84 <- spTransform(meuse.utm, CRS("+proj=lonlat
>>> +ellps=WGS84"))
>>> meuse.lonlatWGS84intl  <- spTransform(meuse.utm, CRS("+proj=lonlat
>>> +ellps=intl"))
>>> coordinates(meuse.lonlatWGS84WGS84)[1:3,]
>>>           x        y
>>> [1,] 5.759053 50.99238
>>> [2,] 5.758380 50.99190
>>> [3,] 5.760373 50.99171
>>>
>>> coordinates(meuse.lonlatWGS84intl)[1:3,]
>>>           x        y
>>> [1,] 5.759053 50.99238
>>> [2,] 5.758380 50.99190
>>> [3,] 5.760373 50.99171
>>>
>>> meuse.lonlatintlWGS84 <- spTransform(meuse.utmintl, CRS("+proj=lonlat
>>> +ellps=WGS84"))
>>> meuse.lonlatintlintl  <- spTransform(meuse.utmintl, CRS("+proj=lonlat
>>> +ellps=intl"))
>>> coordinates(meuse.lonlatintlWGS84)[1:3,]
>>>           x        y
>>> [1,] 5.759053 50.99238
>>> [2,] 5.758380 50.99190
>>> [3,] 5.760373 50.99171
>>>
>>> coordinates(meuse.lonlatintlintl)[1:3,]
>>>           x        y
>>> [1,] 5.759053 50.99238
>>> [2,] 5.758380 50.99190
>>> [3,] 5.760373 50.99171
>>>
>>> Test with Global Mapper
>>>
>>> UTM ED50 272561.0 5654094 to lon,lat:
>>> WGS84   5° 45' 27.5632" E       50° 59' 29.5996" N
>>> ED50    5° 45' 32.5898" E       50° 59' 32.5665" N
>>>
>>> UTM WGS84 272571.9 5653978 to lon,lat:
>>> WGS84   5° 45' 32.5900" E       50° 59' 32.5621" N
>>> ED50    5° 45' 37.6165" E       50° 59' 35.5289" N
>>>
>>> Also, utmWGS84 to utmED50 is different in GlobalMapper
>>>
>>> UTM WGS84  272571.9     5653978
>>> UTM ED50   272662.989   5654181.171
>>>
>>>
>>> Thanks,
>>>
>>> Agus
>>>
>>
>> --
>> Roger Bivand
>> Economic Geography Section, Department of Economics, Norwegian School of
>> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
>> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
>> e-mail: Roger.Bivand at nhh.no
>>
>
>
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no