PP Proleptic Gregorian

[rcomer]

🚀 Pull Request

Description

A proposal to resolve #3561, hopefully without creating breaking changes for the user. Quoting myself from that issue:

Noting that

• 1st January 1970 is the same date in the Proleptic Gregorian and Standard calendars
• Any time point expressed as hours since 1st January 1970 is also therefore the same in the Proleptic Gregorian and Standard calendars

So, once we have converted our pp time headers into a coordinate with unit "hours since 1970-01-01 00:00:00", it is somewhat arbitrary whether we label it "proleptic_gregorian" or "standard". Would it therefore be reasonable to use the Proleptic Gregorian calendar to do the conversion from pp-header to coordinate, but still label the coordinate with the "standard" calendar? It would be a bit of a fudge, but would mean any change in behaviour would be limited to dates before 15th October 1582, and therefore a bugfix.

Similarly, when saving we can convert a standard calendar unit to Proleptic Gregorian with the cf_units.Unit.change_calendar method, and then use that new unit to do the num2date conversion.

I've done a fair bit of test wrangling. I think the last few failures reflect that the years 0 and 1 are genuinely different in the two calendars, so the expected values could just be updated. Leaving that for now and opening as draft to see what people think.

---

Consult Iris pull request check list

[trexfeathers]

📢 @SciTools/iris-devs 📢

Need to make sure we get this one right. If you have any experience with calendars / PP-format it would be great to have your informed opinion ❤

[rcomer]

Full disclosure: this is my first time doing anything in fileformats so, um 😨

[trexfeathers]

@nhsavage are you able to test this branch against any of your scripts, to see if it addresses your concerns?

[nhsavage]

I am struggling a bit to do this - the problem is that the bug mainfested when using ANTS not directly in Iris. trying to mimic it just in Iris doesn't seem to work and at the moment, ANTS is at Iris 2.3.0. Would it be useful to see if I can replicate the problem with ANTS 1.0 and go from there?

[nhsavage]

Would it be useful to see if I can replicate the problem with ANTS 1.0 and go from there? OK. I have been able to come up with a simple piece of ANTs code which demonstrates the problem. However, now reading the ANTS code it is clear that fixing this issue is simply a way to unblock the changes needed in ANTS. The ANTS code which Generates pp fields from a cube. has:

if time_coords[0].units.calendar in ["proleptic_gregorian", "standard"]:
       time_coord = ants.utils.cube.find_time_coordinates(cube)[0]
       new_units = cf_units.Unit(time_coord.units.name, calendar="gregorian")
       time_coord.units = new_units

This willl be wrong for proleptic gregorian calendars before 1582 (more specifically, it changes the calendar of the base time "days since 0001-01-01" without changing the time points to ensure that they match. If this PR passes, then some careful refactoring of this part of the ANTS code would then be needed - specifically to raise an error if a standard or gregorian calendar has data before 1582. I probably need to now make a pure Iris version of this to see if I can replicate the bug without ANTS and understand impacts of this change

Our current workaround is to update the basis time of the NetCDF file prior to running the ants.save

new=cf_units.Unit('days since 1850-1-1', calendar='proleptic_gregorian') sst.coord('time').convert_units(new)

which avoids the bug in ANTS as then gregrian and proleptic gregorian are the same

[nhsavage]

so I've boiled the Iris part of the problem down to something simple. Take any netCDF file with a proleptic gregorian calendar and run:

cube=iris.load_cube(infile)
# save to a pp file
iris.save(cube,'sst_out_fix.pp')

if you look at the output of this, then the calendar is set to 0 which is not a valid. Header from pumf:

(13) lbtim : 20

it should be 21 as the calendar is proleptic gregorian.

[trexfeathers]

if you look at the output of this, then the calendar is set to 0 which is not a valid. Header from pumf:

(13) lbtim : 20

it should be 21 as the calendar is proleptic gregorian.

I can't replicate this, I get the below, both on this PR and with v3.4.0:

  (13) lbtim   :                       0

Here's my code

from cf_units import CALENDAR_PROLEPTIC_GREGORIAN, Unit
import numpy as np
from numpy import random

import iris
from iris.coords import AuxCoord, DimCoord
from iris.cube import Cube

dim_names_lens = [("longitude", 2), ("latitude", 3)]
dim_coords = [DimCoord(np.arange(l), standard_name=n) for n, l in dim_names_lens]

time_unit = Unit("days since 2001-01-01", calendar=CALENDAR_PROLEPTIC_GREGORIAN)
time_coord = AuxCoord(1, standard_name="time", units=time_unit)

prep_cube = Cube(
    random.rand(*[c.shape[0] for c in dim_coords]),
    standard_name="sea_surface_temperature",
    units="K",
    dim_coords_and_dims=[(c, ix) for ix, c in enumerate(dim_coords)],
    # Make this a scalar coord.
    aux_coords_and_dims=[(time_coord, None)],
)

infile = "sst_in.nc"
iris.save(prep_cube, infile)

cube=iris.load_cube(infile)
# save to a pp file
iris.save(cube,'sst_out_fix.pp')

[rcomer]

Looks like I missed a bit 👀

[nhsavage]

First of all LBTIM is defined thus (in my well thumbed copy of UMDP F03, excellent reading if you can't sleep) BTIM is coded as (100IA + 10IB + IC) where: for time aggregated fields IA is the time interval in hours between the individual fields from which the mean was computed [zero in both our cases. yours isn't a time mean, the netCDF I read doesn't include this info

IB = 2 if the field is a time mean between T1 and T2 [my data is time mean yours isn't]

IC (the only relevant part for this PR)

= 0 if ’model time’ is used for T1 and T2 (i.e. only day number, hour and minute are set). – = 1 if the Proleptic Gregorian calendar is used for T1 and T2. This calendar is equivalent to the Gregorian calendar for dates after 1582. Applications producing data using the Proleptic Grego- rian calendar should extend back in time to BC years in a way that make the date continuous: year 1 is 1 AD, year 0 is 1 BC, year -1 is 2 BC etc. (Note some ancillary files use year 0 in this calendar for ancillary fields that do not vary with time.) – = 2 if the ’360-day’ calendar (i.e. 12 30-day months) is used for T1 and T2. – = 3 if ’model time’ is used for T1 and T2 (i.e. only day number, hour and minute are valid; year, month and day in month are to be ignored if set). – = 4 if the ’365-day’ (no leap year) calendar is used for T1 and T2

so ignore whether IB is set or not and only consider the IC - this should be 1 in your case and my case.

To do this all in memory, it should be possibel to use the Iris routine iris.fileformats.pp.as_fields adn check the the LBTIM value in that, which should make it a more useful unit test

[nhsavage]

this code can replace the save and pumf steps

fields = iris.fileformats.pp.as_fields(prep_cube)

for field in fields:
    print(field.lbtim)

currently gives 0 (no calendar defined) which is wrong

[rcomer]

My latest commit gives lbtim of 1 with @trexfeathers' code.

[trexfeathers]

My latest commit gives lbtim of 1 with @trexfeathers' code.

I can replicate this. Thanks @nhsavage for the help debugging.

I'm concerned about whether I would have spotted this even if I were to carefully read UMDP F03 (and any other applicable resources); it's just not where my talents lie. I can do my best (resource allowing), but it would be great to have more users trying this out too (it would be mean to lean exclusively on Nick).

I will see if I can apply the changes here to an older version of Iris to allow more testing by people close to the ANTS tool.

[github-actions[bot]]

In order to maintain a backlog of relevant PRs, we automatically label them as stale after 500 days of inactivity.

If this PR is still important to you, then please comment on this PR and the stale label will be removed.

Otherwise this PR will be automatically closed in 28 days time.

[trexfeathers]

Testing this is still on my personal to-do list

[nhsavage]

I note that #3561 is about to be closed as stale. Is there an update on this?

[trexfeathers]

I note that #3561 is about to be closed as stale. Is there an update on this?

Sorry, no.

Nothing unique to this topic, but it becoming increasingly difficult to deliver improvements to Iris:

• A gradual increase in the complexity of the codebase.
• A marked increase in frequency of upstream breaking changes - unplanned work ™. We have discovered 3 in the last two weeks, and I have personally needed to 'down tools' many times in the last year.

As Iris 'primary lead' I am thinking carefully about changes we might make to our design and our working practices to free up more resource/focus for Iris improvements.

[nhsavage]

thanks for clarifying that. I will let #3561 close then. It's unfortunate, but I understand. Is it worth closing this PR?

[trexfeathers]

Is it worth closing this PR?

Good point. Since @rcomer has put in the work I would say it ought to get higher priority than other things (which might still be abstract ideas). I will try to treat it with more importance going forward.

[rcomer]

I think the reason this stalled was that we didn't have a reproducing example to demonstrate the problem described in the issue. So I cannot actually prove that this fixes it. If someone could provide such an example, we could pick this up again.

[arjclark]

From the ants end, I'm not worried about you testing it with a backported version. We'll catch it up (and many other fixes!) next time we update our dependencies. Original reporter in #3561 is long moved on, so @nhsavage is probably best person to know what's needed for testing. Can confirm from that original ticket that climate calendars under the "proper" calendar are starting to appear (the standard vs gregorian naming change since Iris 3.3 caught out the afterburner climate model monitor package recently).

[nhsavage]

I'll try and manufacture some data. I am struggling to remember exactly the circumstances it failed in though. I think it was when the base time of the calendar was 1/1/1 AD in some Aussie files. I need to think exactly what we want to capture here.

[nhsavage]

making some test data with mule doesn't seem to replicate the issue. Next I am going to go back to the data that gave use the problems in the first place SST and SeaIce cmip5 ACCESS 1.3.

[nhsavage]

I have tried and failed to reproduce the problem. It is so long ago now that probably most people with this issue have implemented work arounds, so I am wondering if we should just close this now as won't fix

[rcomer]

I agree there is no sense keeping this open if we can no longer reproduce the problem.

@arjclark do you have an opinion on this?

[arjclark]

@hdyson - I think I'm happy to let this lapse and get closed if you are?

[hdyson]

We have a longstanding workaround in ANTS for this which we've not needed to update in a while (it does follow our "workaround" working practice so a grep for "3561" is sufficient to find it - I don't want to link to the trac location since we're mid migration to GitHub). So I'm happy to let this be closed - it's not something that's actively causing us issues.

I will add that we do have unittests for our workaround. If this does become a priority due to other users, I think we should be able to adapt those unittests to generate an example if that would be helpful? It looks like the trick is to have a time coordinate with " since 0000-01-01".

[rcomer]

OK, I will close this but keep the branch in case we ever want to pick this up again.