[R-sig-Geo] Different result from sptransform() and ptransform()

[Roger Bivand]

On Sat, 28 Jan 2012, Agustin Lobo wrote:

> Right, grid1km is the polygon object from whence cents comes from and
> I copied the wrong line.
> Sorry for the confusion.

The confusion is in proj4 in the ptransform() function, as:

> centslonlat<- spTransform(cents,CRS("+init=epsg:4326"))
> coordinates(centslonlat)[1:3,]
             X        Y
[1,] 1.315311 42.23956
[2,] 1.327326 42.24060
[3,] 1.339342 42.24163
> ptransform(data = coordinates(cents),
  src.proj = CRSargs(CRS("+init=epsg:3035")),
  dst.proj = "+proj=longlat +ellps=WGS84 +datum=WGS84"
  )[1:3,]*57.29577951308232
          [,1]     [,2] [,3]
[1,] 1.315311 42.23956    0
[2,] 1.327326 42.24060    0
[3,] 1.339342 42.24163    0

ptransform() omits to test for geographical coordinates in the dst.proj= 
argument, and consequently does not multiply by RAD_TO_DEG. The correct 
code is in rgdal/src/projectit.cpp, from line 255. ptransform() is more 
recent than project() in proj4, so the fact that it isn't the same (either 
source or destination as longlat) got overlooked.

Hope this clarifies,

Roger

> The part with ptransform() should be:
>
>> a <- ptransform(data = coordinates(cents), src.proj = CRSargs(CRS("+init=epsg:3035")), dst.proj = CRSargs(CRS("+init=epsg:4326")))
>> coordinates(cents)[1:3,]
>           X       Y
> [1,] 3603500 2167500
> [2,] 3604500 2167500
> [3,] 3605500 2167500
>> a[1:3,]
>           [,1]      [,2] [,3]
> [1,] 0.02295650 0.7372194    0
> [2,] 0.02316621 0.7372375    0
> [3,] 0.02337593 0.7372556    0
>
> Agus
>
> 2012/1/28 Rainer Hurling <rhurlin at gwdg.de>:
>> On 28.01.2012 19:30 (UTC+1), Agustin Lobo wrote:
>>>
>>> Hi!
>>> I'm reprojecting and get different results from sptransform() than
>>> from ptransform():
>>> (objec cents in https://sites.google.com/site/filestemp2/home/cents.rda)
>>>
>>>> str(cents,max.level=2)
>>>
>>> Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
>>>   ..@ data       :'data.frame':        234 obs. of  4 variables:
>>>   ..@ coords.nrs : num(0)
>>>   ..@ coords     : num [1:234, 1:2] 3603500 3604500 3605500 3606500
>>> 3607500 ...
>>>   .. ..- attr(*, "dimnames")=List of 2
>>>   ..@ bbox       : num [1:2, 1:2] 3603500 2155500 3620500 2167500
>>>   .. ..- attr(*, "dimnames")=List of 2
>>>   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
>>>
>>>> proj4string(cents)
>>>
>>> [1] " +init=epsg:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000
>>> +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
>>>>
>>>> centslonlat<- spTransform(cents,CRS("+init=epsg:4326"))
>>>> coordinates(centslonlat)[1:3,]
>>>
>>>             X        Y
>>> [1,] 1.315311 42.23956
>>> [2,] 1.327326 42.24060
>>> [3,] 1.339342 42.24163
>>>
>>>> a<- ptransform(data = as.data.frame(coordinates(grid1km)), src.proj =
>>>> CRSargs(CRS("+init=epsg:3035")), dst.proj = CRSargs(CRS("+init=epsg:4326")))
>>
>>
>> Hmm, grid1km is not known here until now?
>>
>> Greetings,
>> Rainer
>>
>>
>>>> coordinates(a)[1:3,]
>>>
>>>               x         y z
>>> [1,] 0.02295650 0.7372194 0
>>> [2,] 0.02316621 0.7372375 0
>>> [3,] 0.02337593 0.7372556 0
>>>
>>> The correct results are those from sptransform(), but I would like to
>>> know what I'm doing wrong with ptansform()
>>>
>>> Thanks
>>>
>>> Agus
>>
>>
>
> _______________________________________________
> 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, NHH Norwegian School of Economics,
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