cfplot can't cope with southern hemisphere rotated pole data, or else nzlam orog is broken

[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.