I am interested in transforming regular coordinates to rotated coordinates similar as done in this post: Manually transforming rotated lat/lon to regular lat/lon?
These functions base the rotation on the coordinates of the south pole. However, I have a dataset where I am provided north pole coordinates instead (CORDEX data). I am stuck on what I should change in the code. Anyone would know what to do?
South pole rotation code:
from math import *
def rotated_grid_transform(grid_in, option, SP_coor):
lon = grid_in[0]
lat = grid_in[1];
lon = (lon*pi)/180 # Convert degrees to radians
lat = (lat*pi)/180
SP_lon = SP_coor[0]
SP_lat = SP_coor[1]
theta = 90+SP_lat # Rotation around y-axis
phi = SP_lon # Rotation around z-axis
theta = (theta*pi)/180
phi = (phi*pi)/180 # Convert degrees to radians
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;
else: # 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
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
print lon_new,lat_new
