Skip to content

spTransform: unable to convert coordinates

3 messages · Martin Ivanov, Roger Bivand

Original
Martin Ivanov
# 
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?

Best regards,
Roger Bivand
#
On Wed, 15 Aug 2012, Martin Ivanov wrote:

            

      
      
      
    

        
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

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Roger Bivand
    

    

      
      #
    

  


  

      
On Wed, 15 Aug 2012, Roger Bivand wrote:

        

            

      
      
      
    

        
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