Issues with contour plot of rotated pole data on native grid

[bnlawrence]

Grenville provided me with a datafile, which I can't upload here, but it's small. I'll send it to you via slack.

The surface altitude should be easily plotted using the standard methods, but it creates garbage.

Either the data is not reporting the coordinates correctly, or cf-plot can't quite cope with this situation.

Platform: macOS-14.7.3-arm64-arm-64bit 
HDF5 library: 1.14.4 
netcdf library: 4.9.2 
udunits2 library: /Users/bnl28/mambaforge/envs/madpy/bin/../lib/libudunits2.dylib 
esmpy/ESMF: not available 
Python: 3.11.11 /Users/bnl28/mambaforge/envs/madpy/bin/python3.11
dask: 2024.8.0 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/dask/__init__.py
netCDF4: 1.7.2 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/netCDF4/__init__.py
psutil: 6.1.1 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/psutil/__init__.py
packaging: 24.2 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/packaging/__init__.py
numpy: 1.26.2 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/numpy/__init__.py
scipy: 1.15.1 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/scipy/__init__.py
matplotlib: 3.10.0 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/matplotlib/__init__.py
cftime: 1.6.4 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/cftime/__init__.py
cfunits: 3.3.7 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/cfunits/__init__.py
cfplot: 3.3.0 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/cfplot/__init__.py
cfdm: 1.11.1.0 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/cfdm/__init__.py
cf: 3.16.2 /Users/bnl28/mambaforge/envs/madpy/lib/python3.11/site-packages/cf/__init__.py

MWE:

import cf
import cfplot as cfp
f = cf.read('nzlam.orog')
ff = f[5]
cfp.mapset(proj='rotated')
cfp.con(ff)

with or without the mapset it's garbage.

[bnlawrence]

So, I think the problem is that in this case the auxiliary coordinates are being ignored.

The problems arise, I think, in method cf_data_assign where in the loop after (my line 3741 why it's not 1543 like this one is a puzzle for another day),

cf-plot/cfplot/cfplot.py

Line 1543 in 0f57800

if f.ref("grid_mapping_name:rotated_latitude_longitude", default=False):

we go into the rotated pole loop where ptype is set to 6. The grid_latitude and grid_longitude are read.

I guess there are theoretically two ways forward possible, if you know your rotated coordinates, then you just ignore these and plot nice rectangular grids, and then draw the coordinates on afterwards, or, if you have auxillary coordinates, you can use those. These files have auxillary coordinates, but the next step in the code is:

cf-plot/cfplot/cfplot.py

Lines 1568 to 1573 in 0f57800

	# Extract auxiliary lons and lats if they exist
	if ptype == 1 or ptype is None:
	if plotvars.proj != "rotated" and not rotated_vect:
	aux_lons = False
	aux_lats = False
	for mydim in list(f.auxiliary_coordinates()):

and it's clear we are ignoring the auxillary coordinates (because ptype=6 doesn't go there). I guess the code assumes we can always do the first method. Assuming that it works for NH files set up the same (to be tested when we get data from Grenville), then the bug must arise after this point.

[bnlawrence]

OK, continuing. At this point we have xpole = 172.0 and ypole=50.0, obtained via:

cf-plot/cfplot/cfplot.py

Lines 1547 to 1548 in 0f57800

	xpole = rotated_pole["grid_north_pole_longitude"]
	ypole = rotated_pole["grid_north_pole_latitude"]

That's kind of interesting, given the NZ region is

Auxiliary coords: 
     latitude(grid_latitude(220), grid_longitude(210)) =
        [[-50.92822893616319, ..., -27.202321288454904]] degrees_north
     longitude(grid_latitude(220), grid_longitude(210)) = 
       [[153.9040741025498, ..., -175.4001007803285]] degrees_east

(those being the real coordinates of interest) and so we're on the anti-polar side of that ... but the CORDEX domain is similar with pole at (321.38, -60.31). I am slightly suspicious that the cordex pole is the anti-pole of the one defined here, and so that might lead to the problems here. I vaguely remember a definitions problem which was previously solved by

nplon = nplon-180
nplat = - nplat

that ended up as a Cartopy issue: SciTools/cartopy#2091

[bnlawrence]

Sadly, I don't think that's the problem here.

[bnlawrence]

So, continuing the detective work (which I am recording in bits because I don't get long on this in any given chunk of time):

The next interesting bit of code which is executed for ptype=6:

   # Extract x and y grid points
        if plotvars.proj == 'cyl':
            xpts = x
            ypts = y
        else:
            xpts = np.arange(np.size(x))
            ypts = np.arange(np.size(y))

which feels wrong, insofar as what this will do is expect to plot something which is rectangular innx,ny, but I bet is then converted somehow using xpole and ypole. It does look like this code hasn't quite gone where it needs to.

I can also confirm that the same problem exists with another orography file from a northern hemisphere file so it seems like the capability itself is broken.

[sadielbartholomew]

Thanks Bryan for the report. I've possibly got this fixed via some bug fixes I have put in or are going to be merged today. I'll check if you can give me the dataset - will ask over on Slack.

[sadielbartholomew]

(Can confirm it isn't the dataset in question - the basic example from the docs doesn't work, displaying a figure but a blank plot with dodgy axes labelling, at v3.3.0 until my bug fix goes in hopefully for v3.4.0, showing this with cfp.mapset(proj="rotated"):

though note if you try without the mapset call it may handle it naturally - see the tested-as-working example https://ncas-cms.github.io/cf-plot/build/rotated_pole.html#example-21-rotated-pole-data-plot which is the same as the broken one above but without the proj setting. But if you want to plot it on its native rotated grid, that is indeed broken in v3.3.0 and I will put in a bug fix shortly.)

[sadielbartholomew]

I have worked out what the issue is, at least with respect to the example from the docs and that data - when I get access to your dataset I can ensure it works for that too. PR to go up shortly - I'll link it to the Issue.

[sadielbartholomew]

Assuming it's OK, I am changing the title of the Issue so that I can reference the specific scope in the changelog for the new release to avoid writing a separate Issue out.