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

[Roger Bivand]

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