[R-sig-Geo] spTransform: unable to convert coordinates

[Roger Bivand]

On Wed, 15 Aug 2012, Roger Bivand wrote:

> On Wed, 15 Aug 2012, Martin Ivanov wrote:
>
>> Dear R users,
>> I have a rotated latitude-longitude grid with coordinates of the north pole 
>> POLE_LAT=39.25 and POLE_LON=-162 on a sphere with radius 6371 km. My data 
>> are on a conventional longlat grid on the same sphere. I found that I can 
>> transform my data to the rotated grid by means of proj4's ptransform. For 
>> example:
>> crds <- matrix(data=c(9.05, 48.52), ncol=2); # a pair of geographical 
>> coordinates
>> 
>> rot_crds <- 180/pi*ptransform(data=crds*pi/180, dst.proj="+proj=longlat 
>> +ellps=sphere +no_defs",
>>  src.proj="+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 
>> +lon_0=180 +ellps=sphere +no_defs")[1:2]
>> > rot_crds
>> [1] -5.917698 -1.871950 # these are the correct coordinates in the rotated 
>> grid
>> You may notice that to get the correct result I have to swap the source and 
>> destination CRS.
>> 
>> I would be happy if I could reproduce that result with the sp function 
>> spTransform, but I have no luck:
>> spPoint <- SpatialPoints(coords=crds, proj4string=CRS("+proj=longlat 
>> +ellps=sphere +no_defs"));
>> try1 <- spTransform(spPoint, CRS("+proj=ob_tran +o_proj=longlat 
>> +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"));
>> > try1
>> SpatialPoints:
>>      coords.x1 coords.x2
>> [1,]  2.896087   1.37327# this result is not correct
>>
>>  I also tried with swappped source and destination:
>> 
>> > spPoint <- SpatialPoints(coords=crds, proj4string=CRS("+proj=ob_tran 
>> +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere 
>> +no_defs"));
>> + try2 <- spTransform(spPoint, CRS("+proj=longlat +ellps=sphere 
>> +no_defs"));
>> + try2
>> SpatialPoints:
>>      coords.x1 coords.x2
>> [1,] -170.7407 -46.63283 # this is also incorrect
>> 
>> My question is, is it possible to use spTransform for the kind of 
>> coordinate conversion I need? Or do I have to stick to
>> ptransform of proj4?
>
> The ob_tran is indeed an awkward beast, as witnessed by a lot of traffic on 
> mailing lists, for example http://trac.osgeo.org/gdal/ticket/4285. It is as 
> you say counter-intuitive. It is typically used by modellers rather than 
> cartographers.
>
> I see that (using paste to try to overcome mailer linebreaks):
>
> crds <- matrix(data=c(9.05, 48.52), ncol=2)
> library(rgdal)
> #Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
> project(crds*(pi/180), paste("+proj=ob_tran +o_proj=longlat",
> "+o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"),
> TRUE)
>
> gives the result you stipulate (-5.917698  -1.87195), as does:
>
> spPoint <- SpatialPoints(coords=crds*(pi/180),
> proj4string=CRS(paste("+proj=ob_tran +o_proj=longlat +o_lon_p=-162",
> "+o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs")))
> spTransform(spPoint, CRS("+proj=longlat +ellps=sphere +no_defs"))
>
> by casting the input coordinates to radians first, but I don't understand 
> why.

It is the reversal of source and destination projection definitions in 
rgdal::project() and spTransform() methods that is defeating the logic - 
since the source definition is not +proj=longlat but +proj=ob_trans, the 
input coordinates are not multiplied by DEG_TO_RAD.

I'll look at adding an argument to rgdal::project() and spTransform() 
methods to permit the user to say that +proj=ob_tran should be respected, 
which will then give a warning and be ignored if no ob_tran is found. Does 
that seem like a resolution? Can you install rgdal from source from 
R-Forge (once I'm done), or would you prefer a Windows binary from 
Win-builder?

Roger

> I've seen reports of the use of -m=180/pi for command line proj, but 
> cannot see how to apply the scaling factor within a proj4 string.
>
> If this needs to be done a lot, we'll need a way of multiplying the 
> coordinates by (pi/180) for Spatial* objects, possibly by extending elide() 
> methods in maptools.
>
> Does this make sense?
>
> Roger
>
>> 
>> Best regards,
>> 
>> 
>
>

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