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

Question

First I should clarify I don't have previous experience with the field, so I don't know the technical terminology. My question is as follows:

I have two weather datasets:

• The first one has the regular coordinate system (I don't know if it has an specific name), ranging from -90 to 90 and -180 to 180, and the poles are at latitudes -90 and 90.

• In the second one, although it should correspond to the same region, I noticed something different: latitude and longitude were not the same, as they have another reference point (in the description is called a rotated grid). Together with the lat/lon pairs, comes the following information: southern pole lat: -35.00, southern pole lon: -15.00, angle: 0.0.

I need to transform the second pair of lon/lat to the first one. It could be as simple as add 35 to the latitudes and 15 to the longitudes, since the angle is 0 and it seems a simple shifting, but I'm not sure.

Edit: The information I have about the coordinates is the following

http://rda.ucar.edu/docs/formats/grib/gribdoc/llgrid.html

Apparently, the second coordinate system is defined by a general rotation of the sphere

"One choice for these parameters is:

• The geographic latitude in degrees of the southern pole of the coordinate system, thetap for example;

• The geographic longitude in degrees of the southern pole of the coordinate system, lambdap for example;

• The angle of rotation in degrees about the new polar axis (measured clockwise when looking from the southern to the northern pole) of the coordinate system, assuming the new axis to have been obtained by first rotating the sphere through lambdap degrees about the geographic polar axis, and then rotating through (90 + thetap) degrees so that the southern pole moved along the (previously rotated) Greenwich meridian."

but still I don't know how to convert this to the first one.

• Kirk Kuykendall: So is this GRIB data? If so, maybe we need a grib tag.
• user36959: @skd the ECMWF links do not seem to be valid. Can you edit ?
• skd: @gansub I've edited the links. I don't know if the information is exactly the same since it has been a long time, but I believe the new link can provide some context for future reference.
• FaCoffee: @skd when you say angle=0.0, do you mean the bearing? I have a netcdf file with the rotated pole coordinates, but there's no mention of any angle.
• skd: @CF84 I'm actually not sure. I guess if there is no mention of the angle then it is the same as angle = 0

Answer 1

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(θ)) - φ

• whuber: The idea's good, but some of the details need fixing. lon0 = -15, not +15. All three lines in the expansion of the matrix product are incorrect. ATan2 (or its equivalent) must be used, modified to return any reasonable longitude when x=y=0. Note that because x^2+y^2+z^2 = 1, at the end you get simply lat = Arcsin(z).
• mkadunc: Thanks. I fixed the answer to at least make the math correct. The rotations should now match the description in the CS definition, but it's hard to be certain about their sign without an example (other than the position of the south pole).
• whuber: Well done! I'm surprised this reply isn't getting more votes, because it provides useful and hard-to-find material.
• skd: This is indeed very hard to find material, thank you very much for the answer. I ended up using this software code.zmaw.de/projects/cdo to convert from a rotated grid to a regular grid. My guess is that it first transforms the coordinates as in this answer and then interpolates them in order to give the results at the points of a rectangular grid. Although a little late, I leave this her for future reference.
• mkadunc: @alfe I am not an expert on Bloch spheres, but the principle looks very similar to what I've done, but instead of converting to cartesian space with 3 real coordinates, the hint suggests converting to a space with 2 imaginary coordinates (which means 4 real components) and executing the rotation there. Triggered by your comment, I put all the expressions together and added a result in which the intermediate cartesian step is not apparent any more.

Answer 2

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

Answer 3

https://www.mathworks.com/matlabcentral/fileexchange/43435-rotated-grid-transform

PYTHON:

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;
rotated_grid_transform((0,0), 1, (0,30))

• Berke Şentürk: It works with small differences for example I'm using a cordex data with domain EUR-11 from model DMI-HIRHAM5, and inputs were ``` example_coords_from_cordex={"lat": 40.977142333984375 "lon": 28.675996780395508 "rlat": -9.1850004196167 "rlon": 8.145000457763672}``` and runned the function like with inputs SP_coor=[-162.0,39.25] grid_in=[28.675996780395508,40.977142333984375] and >>>rotated_grid_transform(grid_in, 1,SP_coor) results; lon:-8.145000071911468, lat:9.185000943924845
• Berke Şentürk: The results surprisingly inverse idk why however precision looks acceptable right?
• Li Yu: I think the SP_coor here actually refers to the north pole coordinates, so the correct one should be [-162-180, -39.25]

Answer 4

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.

Answer 5

I recommend using gdal. One of the main geodata libraries, which is used under the hood of every other geoprocessing tool.

Given you have a model run like *evspsbl_EUR-11_ICHEC-EC-EARTH_rcp85_r12i1p1_CLMcom-CCLM4-8-17_v1_day_20160101-20201231.nc. This is a climate Projection from the Cordex Initiative. (See meta data here: https://cera-www.dkrz.de/WDCC/ui/cerasearch/header?acronym=CXEU11CLECr8121001CLD1evsbl) and you find something like this in the metadata:

double rlat(rlat) ; 
 rlat:units = "degrees" ; 
 rlat:axis = "Y" ; 
 rlat:long_name = "latitude in rotated pole grid" ; 
 rlat:standard_name = "grid_latitude" ; 
double rlon(rlon) ; 
 rlon:units = "degrees" ; 
 rlon:axis = "X" ; 
 rlon:long_name = "longitude in rotated pole grid" ; 
 rlon:standard_name = "grid_longitude" ; 
int rotated_latitude_longitude ; 
 rotated_latitude_longitude:grid_mapping_name = "rotated_latitude_longitude" ; 
 rotated_latitude_longitude:grid_north_pole_latitude = 39.25 ; 
 rotated_latitude_longitude:grid_north_pole_longitude = -162. ; 
 rotated_latitude_longitude:north_pole_grid_longitude = 0. ; 

You could use gdal straight from the commandline as:

gdalwarp -geoloc NETCDF:evspsbl_EUR-11_ICHEC-EC-EARTH_rcp85_r12i1p1_CLMcom-CCLM4-8-17_v1_day_20160101-20201231.nc:evspsbl outputfilename.tif

You can check out the options for gdalwarp in the documentation here: https://gdal.org/programs/gdalwarp.html

However if you are keen on batch processing/calling everything from within python, this should work. Note that the filename usually gives away a lot of information. In this case evspsbl refers to the datatype within the data. Here evapotranspiration. EUR-11 is the Cordex Domain (see: https://cordex.org/domains/)

from osgeo import gdal
netcdf_file = 'evspsbl_EUR-11_ICHEC-EC-EARTH_rcp85_r12i1p1_CLMcom-CCLM4-8-17_v1_day_20160101-20201231.nc'
# optional: complete dataset & inspect the metadata
ds = gdal.Open(netcdf_file)
for key, value in ds.GetMetadata().items():
 print("{:35}: {}".format(key, value))
#open the specified subset
evspsbl = gdal.Open('NETCDF:'+netcdf_file+':evspsbl')
# warp into EPSG:4326, 
gdal.Warp('outputfile.tif', evspsbl, dstSRS='EPSG:4326', geoloc=True)

I'm not entirely sure what the -geoloc flag does, but it works... . dstSRS defines the destination CRS. Also check out the gdal Python API https://gdal.org/python/. You can now inspect the .tif in QGIS or your prefered method of inspecting and will see it has different bands. In my case those refer to the timestamps of the original .nc file.

Answer 6

Which software are you using? Every GIS software will have the facility to show you the current cordinate system/projection information. , which can help you in getting the name of your current coordinate system.

Additionally, if you're using ArcGIS, you can use the Project tool to re-project the second dataset, importing settings from the first one.

• skd: Unfortunatly I'm not using any software. Those are just grid datasets and they come with the following information: - For the first one: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/… - For the second one: ecmwf.int/publications/manuals/d/gribapi/fm92/grib1/detail/…
• ujjwalesri: Since the angle of rotation is 0, I think a simple translation should align the second dataset to the first one, like you said adding 15 to X and 35 to Y

Answer 7

I have developed an asp.net page for converting coordinate from rotated to non-rotated based on CORDEX domains.

It based on above methods. You can use it freely in this link:

Manually transforming rotated lat/lon to regular lat/lon

• Sohrab kolsoomi ayask: Cordex Data Extractor is a windows desktop software for extracting data from CORDEX NetCDF file. Cordex Data Extractor doesn't need help file because of all processes have been done in behind codes and the user just enters dates, coordinate, and variable name. Please watch this video: youtu.be/RmpZblZPXjI agrimetsoft.com/cordexDataExtractor.aspx

Answer 8

The cartopy module for python has a rotated pole facility. It understands the transformations it seems.

https://scitools.org.uk/cartopy/docs/latest/gallery/rotated_pole.html#sphx-glr-gallery-rotated-pole-py

Might be helpful I guess.

Answer 9

Since I really liked @Magnus idea of using GDAL, but also was missing some answer using R in general, I'd like to provide another option hereby.

Let's get some rotated grids first, e.g. CanRCM4 data covering Canada (see docs).

# get some file with a rotated grid
f <- "snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc"
download.file(url = paste0("https://crd-data-donnees-rdc.ec.gc.ca/CCCMA/products/CanSISE/output/CCCma/CanRCM4/NAM-44_CCCma-CanESM2_historical-r1/mon/atmos/snw/r1i1p1/", f),
 destfile = f,
 mode = "wb")
# gdalinfo, basically
terra::describe(f)[46:55]
#> [1] " rlat#long_name=latitude in rotated pole grid" 
#> [2] " rlat#standard_name=grid_latitude" 
#> [3] " rlat#units=degrees" 
#> [4] " rlon#axis=X" 
#> [5] " rlon#long_name=longitude in rotated pole grid" 
#> [6] " rlon#standard_name=grid_longitude" 
#> [7] " rlon#units=degrees" 
#> [8] " rotated_pole#grid_mapping_name=rotated_latitude_longitude"
#> [9] " rotated_pole#grid_north_pole_latitude=42.5" 
#> [10] " rotated_pole#grid_north_pole_longitude=83"
# inspect raw data
r1 <- terra::rast(f)
r1
#> class : SpatRaster 
#> dimensions : 130, 155, 12 (nrow, ncol, nlyr)
#> resolution : 0.4400001, 0.44 (x, y)
#> extent : -34.1, 34.1, -28.82, 28.38 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> source : snw_NAM-44_CCCma-CanESM2_historical-r1_r1i1p1_CCCma-CanRCM4_r2_mon_195001-195012.nc:snw 
#> varname : snw (Surface Snow Amount) 
#> names : snw_1, snw_2, snw_3, snw_4, snw_5, snw_6, ... 
#> unit : kg m-2, kg m-2, kg m-2, kg m-2, kg m-2, kg m-2, ... 
#> time (days) : 1954年08月08日 to 1955年07月08日
terra::plot(r1[[1]])

This will be the starting point when working with this kind of data, noticing a significant misplacement in comparison to EPSG:4326. Making use of GDAL we can now warp this raster data to the desired CRS quite conveniently:

# get some reference data
adm <- geodata::gadm("Canada", path = tempdir())
# gdalwarp
sf::gdal_utils(util = "warp", 
 source = f, 
 destination = "r_wgs84.nc", 
 options = c("-t_srs", "EPSG:4326"))
# inspect warped data and verify georeference
r2 <- terra::rast("r_wgs84.nc")
r2
#> class : SpatRaster 
#> dimensions : 113, 262, 12 (nrow, ncol, nlyr)
#> resolution : 0.5647745, 0.5642747 (x, y)
#> extent : -170.9646, -22.99372, 12.11698, 75.88002 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
#> sources : r_wgs84.nc:Band1 
#> r_wgs84.nc:Band2 
#> r_wgs84.nc:Band3 
#> ... and 9 more source(s)
#> varnames : Band1 (Surface Snow Amount) 
#> Band2 (Surface Snow Amount) 
#> Band3 (Surface Snow Amount) 
#> ...
#> names : Band1, Band2, Band3, Band4, Band5, Band6, ... 
#> unit : kg m-2, kg m-2, kg m-2, kg m-2, kg m-2, kg m-2, ...
terra::plot(r2[[1]])
terra::plot(adm, add = TRUE)

^{Created on 2024年03月27日 with reprex v2.1.0}