Iris does not recognise x and y coordinates of netCDF files with rotated mercator projection

[thomascrocker]

🐛 Bug Report

When loading a model that uses the rotated mercator projection, iris does not recognise the co-ordinate system and the resulting x and y coordinates of the iris cube are anonymous. No warnings or errors are printed.

How To Reproduce

Steps to reproduce the behaviour: Load this file for example: /project/champ/data/cordex/output/SAM-44/ICTP/MOHC-HadGEM2-ES/historical/r1i1p1/RegCM4-3/v4/day/pr/v20151029/pr_SAM-44_MOHC-HadGEM2-ES_historical_r1i1p1_ICTP-RegCM4-3_v4_day_1970010112-1975123012.nc

> c = iris.load_cube('/project/champ/data/cordex/output/SAM-44/ICTP/MOHC-HadGEM2-ES/historical/r1i1p1/R
    ...: egCM4-3/v4/day/pr/v20151029/pr_SAM-44_MOHC-HadGEM2-ES_historical_r1i1p1_ICTP-RegCM4-3_v4_day_197601011
    ...: 2-1980123012.nc')

> print(c)
precipitation_flux / (kg m-2 s-1)                (time: 1800; -- : 199; -- : 189)
    Dimension coordinates:
        time                                          x          -         -
    Auxiliary coordinates:
        latitude                                      -          x         x
        longitude                                     -          x         x

Expected behaviour

The x and y co-ordinates of the cube should be recognised as the projection x and y coordinates and have the associated coordinate system metadata associated with them.

Environment

• OS & Version: RHEL7
• Iris Version: 3.1.0

Additional context

For info:

$ ncdump -h /project/champ/data/cordex/output/SAM-44/ICTP/MOHC-HadGEM2-ES/historical/r1i1p1/RegCM4-3/v4/day/pr/v20151029/pr_SAM-44_MOHC-HadGEM2-ES_historical_r1i1p1_ICTP-RegCM4-3_v4_day_1970010112-1975123012.nc 
netcdf pr_SAM-44_MOHC-HadGEM2-ES_historical_r1i1p1_ICTP-RegCM4-3_v4_day_1970010112-1975123012 {
dimensions:
	time = UNLIMITED ; // (2160 currently)
	time_bnds = 2 ;
	y = 199 ;
	x = 189 ;
variables:
	double time_bnds(time, time_bnds) ;
		time_bnds:units = "days since 1949-12-01 00:00:00 UTC" ;
		time_bnds:calendar = "360_day" ;
	double lat(y, x) ;
		lat:standard_name = "latitude" ;
		lat:long_name = "Latitude at cross points" ;
		lat:units = "degrees_north" ;
	double lon(y, x) ;
		lon:standard_name = "longitude" ;
		lon:long_name = "Longitude at cross points" ;
		lon:units = "degrees_east" ;
	double x(x) ;
		x:standard_name = "projection_x_coordinate" ;
		x:long_name = "x-coordinate in Cartesian system" ;
		x:units = "m" ;
		x:axis = "X" ;
	double y(y) ;
		y:standard_name = "projection_y_coordinate" ;
		y:long_name = "y-coordinate in Cartesian system" ;
		y:units = "m" ;
		y:axis = "Y" ;
	int crs ;
		crs:grid_mapping_name = "rotated_mercator" ;
		crs:latitude_of_projection_origin = -22. ;
		crs:longitude_of_projection_origin = -59. ;
		crs:scale_factor_at_projection_origin = 1. ;
		crs:semi_major_axis = 6371229. ;
		crs:inverse_flattening = 0. ;
		crs:_CoordinateTransformType = "Projection" ;
		crs:_CoordinateAxisTypes = "GeoX GeoY" ;
		crs:false_easting = -25000. ;
		crs:false_northing = -25000. ;
	double time(time) ;
		time:standard_name = "time" ;
		time:long_name = "time" ;
		time:calendar = "360_day" ;
		time:units = "days since 1949-12-01 00:00:00 UTC" ;
		time:bounds = "time_bnds" ;
	float pr(time, y, x) ;
		pr:cell_methods = "time: mean" ;
		pr:standard_name = "precipitation_flux" ;
		pr:long_name = "Precipitation" ;
		pr:units = "kg m-2 s-1" ;
		pr:coordinates = "lat lon" ;
		pr:grid_mapping = "crs" ;

[trexfeathers]

Thanks @thomascrocker. This has come up multiple times in the last few weeks. Unfortunately https://github.com/SciTools/cartopy/issues/1618 blocks us implementing this in our current structure. #4643 offers a solution, but it's a big one - make sure you get your votes in for that!

[thomascrocker]

Thanks @thomascrocker. This has come up multiple times in the last few weeks. Unfortunately SciTools/cartopy#1618 blocks us implementing this in our current structure. #4643 offers a solution, but it's a big one - make sure you get your votes in for that!

Thanks for the heads up, have voted for #4643 !

[vsherratt]

Just to clarify: both available solutions there would still not support rotated_mercator directly, because that is not a recognised grid mapping in the CF conventions.

What they would support is oblique_mercator, which can be used to handle rotated mercator by setting the azimuth to 90.

How important is it to you that iris loads rotated_mercator directly @thomascrocker ? Instead of requiring you to make this modification in the files beforehand.

[thomascrocker]

Just to clarify: both available solutions there would still not support rotated_mercator directly, because that is not a recognised grid mapping in the CF conventions.

What they would support is oblique_mercator, which can be used to handle rotated mercator by setting the azimuth to 90.

How important is it to you that iris loads rotated_mercator directly @thomascrocker ? Instead of requiring you to make this modification in the files beforehand.

Ah I see. Ultimately this is model data that is part of a wider collection of model data available via ESGF. Having to modify the files beforehand is a bit of a pain if using something like ESMValTool for processing, since the model for that software is to work with a centrally downloaded collection of model data which is identical to the "official" files as downloaded from ESGF. My knowledge of projections is not advanced enough to know how easy this would be... but a workaround would be adding a fix to ESMValTool that could modify the rotated_mercator projection to oblique_mercator as part of "on the fly" CMOR fixes that are applied at the point of loading a file. Tagging @nhsavage for info.

[trexfeathers]

Update: Cartopy now includes code supporting Oblique Mercator. Once that feature makes it into a release, this will enable development within Iris.

[trexfeathers]

Big thanks @bjlittle!