Cartopy rotated pole definition is inconsistent with CORDEX definitions

[bnlawrence]

Description

I have followed the examples for generated regional rotated model maps. They work for those defined by CORDEX, for the northern hemisphere, but break with a southern hemisphere example (the Australasia domain).

A simple fix is to manually refine the pole by moving it, but that suggests that either Cartopy or CORDEX have got something wrong.

Code to reproduce

import cartopy.crs as ccrs
from matplotlib import pyplot as plt

# official CORDEX definition
defn = {'RotPole': (321.38, -60.31), 'TLC': (142.16, 33.44), 'Nx': 200, 'Ny':129, 'delta':0.44}

x1, y1 = defn['TLC']
xn = x1 + (defn['Nx']-1)*defn['delta']
yn = y1 - (defn['Ny']-1)*defn['delta']

nplon, nplat = defn['RotPole']

#### this gives the wrong answer

rotated_pole = ccrs.RotatedPole(pole_latitude = nplat, pole_longitude=nplon, central_rotated_longitude=180.)
ax = plt.subplot(211, projection=rotated_pole)
ax.stock_img()
ax.coastlines()
ax.gridlines()
ax.set_title('wrong')

xl1 = x1-180
xl2 = xn-180
yl1, yl2 = y1, yn

ax.set_extent([xl1,xl2,yl1,yl2], crs=rotated_pole)

### this gives the right answer

nplon = nplon-180
nplat = - nplat

rotated_pole = ccrs.RotatedPole(pole_latitude = nplat, pole_longitude=nplon, central_rotated_longitude=180.)
ax = plt.subplot(212, projection=rotated_pole)
ax.stock_img()
ax.coastlines()
ax.gridlines()
ax.set_title('correct ')

xl1 = x1-180
xl2 = xn-180
yl1, yl2 = y1, yn

ax.set_extent([xl1,xl2,yl1,yl2], crs=rotated_pole)

plt.show()

Full environment definition

python 3.9.13

Cartopy version

0.20.3

[dopplershift]

Well this line 'lon_0', 180 + pole_longitude looks suspicious given how you're able to correct things:

https://github.com/SciTools/cartopy/blob/11987150a8350b0d66f9f18f7bee1dccc47e4ea0/lib/cartopy/crs.py#L1928-L1933

@pelson This dates back to your initial commit for Cartopy, any ideas why you included the shift? Otherwise, I'm inclined to remove it and hand things directly to PyPROJ.

[dopplershift]

Ugh https://github.com/SciTools/cartopy/issues/525#issuecomment-68840970 seems to indicate that this was an intentional design decision, so I'm really not inclined to change it. I'm also getting a bit out of my depth here.

Essentially it seems like meteorology (WMO) decided to define the rotated pole projection differently than cartography (PROJ). Based on this post it also sounds like the netCDF CF conventions are different than the WMO. The different "modes" are also noted here by the OGC.

Maybe the right answer is to add a flag to our RotatedPole class (e.g. use_wmo_definition=True) to make the "escape hatch" more clear. Documenting this discrepency in the docstring would also be a really good idea.

[pelson]

Maybe the right answer is to add a flag to our RotatedPole class (e.g. use_wmo_definition=True) to make the "escape hatch" more clear. Documenting this discrepency in the docstring would also be a really good idea.

I definitely think it can be more clearly documented - and the references I cited back in 2015 be integrated into the docs. In terms of changing the default - I could imagine that causing a lot of pain.

Based on this post it also sounds like the netCDF CF conventions are different than the WMO

I don't know how helpful it is, but the rules implemented in Iris are going to be fairly well battle tested on this:

https://github.com/SciTools/iris/blob/d8f8dac7c69f765c75690ab540d8c961adbb921c/lib/iris/fileformats/_nc_load_rules/helpers.py#L306

and

https://github.com/SciTools/iris-grib/blob/a4bed1a0a671f4ca5770e9df9d3c6c161db4f417/iris_grib/_load_convert.py#L333

[nhsavage]

I also wonder if these conflicting definitions need to be more explitictly recognised in CORDEX community. Perhaps by updating https://cordex.org/wp-content/uploads/2012/11/CORDEX-domain-description_231015.pdf ?

[nhsavage]

Having read the documents linked here, we conclude that Cartopy is functioning correctly but the documentaiton does need some further improvement.

Starting with the CF conventions, Table F.1. Grid Mapping Attributes we have the following 3 relevant parameters

Attribute	Type	Description
grid_north_pole_latitude	N	True latitude (degrees_north) of the north pole of the rotated grid.
grid_north_pole_longitude	N	True longitude (degrees_east) of the north pole of the rotated grid.
north_pole_grid_longitude	N	Longitude (degrees) of the true north pole in the rotated grid.

For more context on this see this OGC document

This has the mappings from CF names to proj:

Identifier proposal	Name proposal	CF-convention name	PROJ name	Cartopy argument
urn:ogc:def:parameter:OGC::111	Latitude of rotated pole	grid_north_pole_latitude	+o_lat_p	pole_latitude
urn:ogc:def:parameter:OGC::112	Longitude of rotated pole	grid_north_pole_longitude	+o_lon_p	pole_longitude
urn:ogc:def:parameter:OGC::113	Axis rotation	north_pole_grid_longitude	+lon_0 antipodal	central_rotated_longitude

Note on PROJ parameters: The value given to PROJ +lon_0 parameter needs to be the axis rotation value with a shift of ±180°. If there is no axis rotation, then +lon_0=180 needs to be specified. Alternatively, the -I option can be avoided by swapping the +lon_0 and +o_lon_p values (the +180 offset stay on +lon_0).

Comparing this to the code in cartopy

it has the following arguments: pole_longitude, pole_latitude, central_rotated_longitude=0.0

these are unfortunately, not the same names as used in the CF conventions. However the mapping through to proj follows the table above, if we assume that these names translate to the CF conventions as:

• pole_latitude = grid_north_pole_latitude
• pole_longitude = grid_north_pole_longitude
• central_rotated_longitude = north_pole_grid_longitude

as the code is:

        proj4_params = [('proj', 'ob_tran'), ('o_proj', 'latlon'),
                        ('o_lon_p', central_rotated_longitude),
                        ('o_lat_p', pole_latitude),
                        ('lon_0', 180 + pole_longitude),
                        ('to_meter', math.radians(1) * (
                            globe.semimajor_axis or WGS84_SEMIMAJOR_AXIS))]

It would however, be good if the docstring for this function could clarify these issues by

• noting the correct mapping from CF names to the arguments of the function
• linking to the OGC description and
• clarify that the WMO definition used in Grib is based on the rotated south pole (as also described in the OGC document) but CF and hence Cartopy use the roatated north pole convention.

[larsbuntemeyer]

Thanks @bnlawrence for pointing this issue out here with the pole definition in Australasia. I had some similar issue as well when creating cartopy plots, that's when i recognized this inconsistency in the pole definition. I think for the other southern hemisphere plots the pole definitions seem to be correct.

[bnlawrence]

Does that mean there is an inconsistency between the way the way some of the SH domains are defined? If some work, and some don't? (I can't remember what I did when I was doing this ... but am keen to track the outcome because at some point I will come back to this space :-) )

[nhsavage]

I have been looking into this and for Australasia it seems to be a typo in the rotated pole. It plots correctly if I use 321.38; 60.31 and not if use the version on the CORDEX page 321.38; -60.31

I double checked this with the Met Office tool lampos

I also had some issues with S. America that I am working on too.

There are some complicated issues as well with the best way to plot the domains on a cartopy map as well with the simplest approach of setting the extent of the map giving misleading results. I hope to be able to fully document this soon.

[larsbuntemeyer]

@nhsavage i found this helpful. And here is how i did it. Of course, it depends a little bit on your projection and if you want to plot gridded data or just the domain as a kind of map.

[nhsavage]

@nhsavage i found this helpful.

thanks, that's basically what I have done, but trying to add some refinements to double check things - in particular the true lat and lon of corners to check against the CORDEX docs.

Polar domains also look a bit odd in current plotting and I am trying to understand those as well. Trying to get a really well researched picture here to give good evidence for any queries about the CORDEX docs