Manually transforming rotated lat/lon to regular lat/lon?

Manually reversing the rotation should do the trick; there should be a formula for rotating spherical coordinate systems somewhere, but since I can't find it, here's the derivation ( ' marks the rotated coordinate system; normal geographic coordinates use plain symbols):

First convert the data in the second dataset from spherical (lon', lat') to (x',y',z') using:

x' = cos(lon')*cos(lat')
y' = sin(lon')*cos(lat')
z' = sin(lat')

Then use two rotation matrices to rotate the second coordinate system so that it coincides with the first 'normal' one. We'll be rotating the coordinate axes, so we can use the axis rotation matrices. We need to reverse the sign in the ϑ matrix to match the rotation sense used in the ECMWF definition, which seems to be different from the standard positive direction.

Since we're undoing the rotation described in the definition of the coordinate system, we first rotate by ϑ = -(90 + lat0) = -55 degrees around the y' axis (along the rotated Greenwich meridian) and then by φ = -lon0 = +15 degrees around the z axis):

x   ( cos(φ), sin(φ), 0) (  cos(ϑ), 0, sin(ϑ)) (x')
y = (-sin(φ), cos(φ), 0).(  0     , 1, 0     ).(y')
z   ( 0     , 0     , 1) ( -sin(ϑ), 0, cos(ϑ)) (z')

Expanded, this becomes:

x = cos(ϑ) cos(φ) x' + sin(φ) y' + sin(ϑ) cos(φ) z'
y = -cos(ϑ) sin(φ) x' + cos(φ) y' - sin(ϑ) sin(φ) z'
z = -sin(ϑ) x' + cos(ϑ) z'

Then convert back to 'normal' (lat,lon) using

lat = arcsin(z)
lon = atan2(y, x)

If you don't have atan2, you can implement it yourself by using atan(y/x) and examining the signs of x and y

Make sure that you convert all angles to radians before using the trigonometric functions, or you'll get weird results; convert back to degrees in the end if that's what you prefer...

Example (rotated sphere coordinates ==> standard geographic coordinates):

• southern pole of the rotated CS is (lat0, lon0)

  (-90°, *) ==> (-35°, -15°)

• prime meridian of the rotated CS is the -15° meridian in geographic (rotated 55° towards north)

  (0°, 0°) ==> (55°, -15°)

• symmetry requires that both equators intersect at 90°/-90° in the new CS, or 75°/-105° in geographic coordinates

  (0°, 90°) ==> (0°, 75°)
  (0°, -90°) ==> (0°,-105°)

EDIT: Rewritten the answer thanks to very constructive comment by whuber: the matrices and the expansion are now in sync, using proper signs for the rotation parameters; added reference to the definition of the matrices; removed atan(y/x) from the answer; added examples of conversion.

EDIT 2: It is possible to derive expressions for the same result without explicit tranformation into cartesian space. The x, y, z in the result can be substituted with their corresponding expressions, and the same can be repeated for x', y' and z'. After applying some trigonometric identities, the following single-step expressions emerge:

lat = arcsin(cos(ϑ) sin(lat') - cos(lon') sin(ϑ) cos(lat'))
lon = atan2(sin(lon'), tan(lat') sin(ϑ) + cos(lon') cos(ϑ)) - φ

In case anyone is interested I've shared a MATLAB script on the file exchange transforming regular lat/lon to rotated lat/lon and vice versa: Rotated grid transform

function [grid_out] = rotated_grid_transform(grid_in, option, SP_coor)

lon = grid_in(:,1);
lat = grid_in(:,2);

lon = (lon*pi)/180; % Convert degrees to radians
lat = (lat*pi)/180;

SP_lon = SP_coor(1);
SP_lat = SP_coor(2);

theta = 90+SP_lat; % Rotation around y-axis
phi = SP_lon; % Rotation around z-axis

phi = (phi*pi)/180; % Convert degrees to radians
theta = (theta*pi)/180;

x = cos(lon).*cos(lat); % Convert from spherical to cartesian coordinates
y = sin(lon).*cos(lat);
z = sin(lat);

if option == 1 % Regular -> Rotated

    x_new = cos(theta).*cos(phi).*x + cos(theta).*sin(phi).*y + sin(theta).*z;
    y_new = -sin(phi).*x + cos(phi).*y;
    z_new = -sin(theta).*cos(phi).*x - sin(theta).*sin(phi).*y + cos(theta).*z;

elseif option == 2 % Rotated -> Regular

    phi = -phi;
    theta = -theta;

    x_new = cos(theta).*cos(phi).*x + sin(phi).*y + sin(theta).*cos(phi).*z;
    y_new = -cos(theta).*sin(phi).*x + cos(phi).*y - sin(theta).*sin(phi).*z;
    z_new = -sin(theta).*x + cos(theta).*z;

end

lon_new = atan2(y_new,x_new); % Convert cartesian back to spherical coordinates
lat_new = asin(z_new);

lon_new = (lon_new*180)/pi; % Convert radians back to degrees
lat_new = (lat_new*180)/pi;

grid_out = [lon_new lat_new];

This transformation can also be computed with proj software (either using command-line or programmatically) by employing what proj calls an oblique translation (ob_tran) applied to a latlon transformation. The projection parameters to be set are:

• o_lat_p = north pole latitude => 35° in the example
• lon_0 = south pole longitude => -15° in the example
• o_lon_p = 0

additionally, -m 57.2957795130823 (180/pi) is required in order to consider projected values in degrees.

Replicating the examples proposed by mkadunc gives the same result (notice that here order is lon lat not (lat,lon), coordinates are typed in standard input, output is marked by =>):

invproj -f "=> %.6f" -m 57.2957795130823 +proj=ob_tran +o_proj=latlon +o_lon_p=0 +o_lat_p=35 +lon_0=-15
0 -90
=> -15.000000   => -35.000000
40 -90
=> -15.000000   => -35.000000
0 0
=> -15.000000   => 55.000000
90 0
=> 75.000000    => -0.000000
-90 0
=> -105.000000  => -0.000000

invproj command is used for converting from "projected" (i.e. rotated) coordinates to geographic, while proj for doing the opposite.