[R-sig-Geo] spTransform: wrong results??
Agustin Lobo
alobolistas at gmail.com
Fri Jul 10 12:18:26 CEST 2009
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...
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
>
--
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
Lluis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
e-mail Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster
More information about the R-sig-Geo
mailing list