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

[Roger Bivand]

On Mon, 30 Jan 2012, Simon Urbanek wrote:

>
> On Jan 30, 2012, at 3:48 AM, Agustin Lobo wrote:
>
>> Roger,
>> While the explanation of the error is clear, I am actually not sure
>> about the practical consequences.
>> What I understand is that we must not use ptransform() by now (and use
>> spTransform() instead), and
>> that ptransform() will be fixed in the future. Am I right?
>>
>
> Roger's explanation is wrong - unlike project() there is no "degrees" 
> argument to so ptransform() expects radians, not degrees - see the 
> supplied examples. This is simply a user error.

Yes, given that the user reads the example on ?ptransform carefully, and 
sees that you are converting degrees to radians on the input for input 
geographical coordinates transforming to projected coordinates. From that 
the astute user might infer that you need to invert the conversion if the 
input coordinates were projected, and the output geographical, but the 
help page doesn't say this.

I'm not aware of any data source providing geographical coordinates with a 
radian metric. I would argue that since there is no established use case 
for the radian metric in geographical data, this is an implementation 
error, because neither input nor output correspond to general practice, 
and it is more like a trap for users than a user error.

I guess we agree to differ?

Best wishes,

Roger

>
> Cheers,
> Simon
>
>
>> Thanks
>>
>> Agus
>>
>>
>> 2012/1/29 Roger Bivand <Roger.Bivand at nhh.no>:
>>> 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
>>>
>>
>>
>
>

-- 
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