FMRC aggregation broken since v7.41

[AndrewWittenberg]

FMRC aggregation and regridding seems to have broken between Ferret v7.4 (4/18/18) and v7.41 (7/9/18). Perhaps due to fixes for #1868, #1649, or #1657? Or something else?

Marking this as urgent, since it breaks something we use -- and we don't have any good fallback versions that don't contain other bugs.

Here's what we had in v7.4:

    NOAA/PMEL TMAP
    FERRET v7.4 (optimized)
    Linux 2.6.32-696.23.1.el6.x86_64 64-bit - 04/18/18
     7-Sep-18 16:02     

yes? let agg_lister = "seq 1997 1999 | xargs -I@YEAR@ sh -c 'ls -1 /archive/rgg/SPEAR/SPEAR_c96_o1_pred_@YEAR@{01,07}/pp_ensemble/atmos/ts/monthly/1yr/atmos.@YEAR@??-??????.t_surf.nc'"
yes? fmrc ens_mean.agg = spawn(agg_lister)
yes? let/unit="degC"/title="NINO3 SST" nino3 = t_surf[x=150w:90w@ave,y=5s:5n@ave] - 273.15
yes? can view; ppl cross 1
yes? set view ul; shade/nokey/trans nino3
yes? set view ur; shade/nokey/trans/lev nino3[gt(tf_times)=tf_cal_t@fmrc]
yes? set view lr; plot/step=con/along=t/nokey nino3[gt(tf_times)=tf_cal_t@fmrc]

The resulting plot:

But in v7.41 through 7.421 (8/27/18) we get:

yes? set view ur; shade/nokey/trans/lev nino3[gt(tf_times)=tf_cal_t@fmrc]
 **ERROR: invalid command: auxiliary variable TF_TIMES has a F axis not found on variable NINO3
shade/nokey/trans/lev nino3[gt(tf_times)=tf_cal_t@fmrc]
Command file, command group, or REPEAT execution aborted

An error!

Let's look at the axes. In v7.4 we had:

yes? sho grid tf_times
    GRID (G014)
 name       axis              # pts   start                end                 subset
 normal    X
 normal    Y
 normal    Z
 TF_LAG_T  MODEL ELAPSED TIME (12 r   15.182               349.18              full
 normal    E
 TF_CAL_F  FORECAST             6 i   01-JAN-1997 00:00    01-JUL-1999 00:00   full

yes? sho grid t_surf
    GRID (G013)
 name       axis              # pts   start                end                 subset
 LON       LONGITUDE          288mi   0.625E               0.625W              full
 LAT       LATITUDE           180 i   89.5S                89.5N               full
 normal    Z
 TF_LAG_T  MODEL ELAPSED TIME (12 r   15.182               349.18              full
 normal    E
 TF_CAL_F  FORECAST             6 i   01-JAN-1997 00:00    01-JUL-1999 00:00   full

But in v7.41 we get:

yes? sho grid tf_times
    GRID (G015)  Forecast Aggregation Grid
 name       axis              # pts   start                end                 subset
 normal    X
 normal    Y
 normal    Z
 TF_LAG_T  MODEL ELAPSED TIME (12 r   15.083               346.92              full
 normal    E
 TF_CAL_F  FORECAST            36 r   01-JAN-1997 10:00    23-NOV-1999 06:00   full

yes? sho grid t_surf
    GRID (G013)  Forecast Aggregation Grid
 name       axis              # pts   start                end                 subset
 LON       LONGITUDE          288mi   0.625E               0.625W              full
 LAT       LATITUDE           180 i   89.5S                89.5N               full
 normal    Z
 AGG_LAG_T MODEL ELAPSED TIME (12 r   15.182               349.18              full
 normal    E
 RUN       FORECAST             6 r   01-JAN-1997 10:00    25-JUN-1999 10:00   full

It used to be that TF_TIMES and the aggregated variable shared their T and F axes. But now they have different axis names, and TF_CAL_F has a different size (36 rather than 6).

A related issue: we need test datasets for aggregations, so that testing, teaching, and bug-reporting are easier & more reproducible. A good place to do this would be in the documentation, e.g.

https://ferret.pmel.noaa.gov/Ferret/documentation/users-guide/commands-reference/DEFINE#_define_data_agg

The docs should either provide the example files "tmp/fcst_1.nc" etc., or provide commands to create them.

[karlmsmith]

I am checking if this may be an issue from changes to add the "spacing" as described at the end of #1663 - I suspect so.

[karlmsmith]

This is due to the changes made in #1663 - obviously needs to be documented better.

A minimal superset axis of both the forecasted times and the times each forecast was made is now used, and the user decides how to fill in the gaps.

So adding in

let nino3fine = nino3[GT=TF_LAG_T@BIN,GF=TF_CAL_F@BIN]

and then plotting nino3fine gives the pyferret plot:

[karlmsmith]

By the way, the @BI in the label gets deleted because I was using a non-Hershey font in PyFerret, and it thinks that is a Hershey font-change command that should be ignored.

[karlmsmith]

This change is in the documentation under https://ferret.pmel.noaa.gov/Ferret/documentation/users-guide/Grids-Regions/GRIDS#_4.2.7 but obviously also needs to be documented under https://ferret.pmel.noaa.gov/Ferret/documentation/users-guide/commands-reference/DEFINE#_define_data_agg along with a script to generate simple data files to play with.

[karlmsmith]

The PyFerret plot using Hershey font and telling it to just plot N=1:37:6 in the last plot.

[AndrewWittenberg]

Thanks very much Karl, especially for taking the time to work with the actual GFDL files! Glad it was a new "feature" rather than a bug -- and your plot is getting close to what we want.

The F-coordinates of the aggregation aren't quite right though. Our forecasts start at 00Z on 1-Jan and 1-July of each year, as indicated by the starting TBOXLO values in the constituent datasets (d=7.1, 7.2, ..., 7.6). You can also see these starting values with, for example, SHOW GRID/T t_surf[d=7.1]. The F-axis coordinates (the "cell centers") should thus be an uneven number of days apart, at exactly 1-jan-1997, 1-jul-1997, 1-jan-1998, 1-jul-1998, etc. But instead they're evenly spaced (exactly 181 days apart), and a bit too close together, so that case #6 is at F=25-jun-1999 instead of F=1-jul-1999.

The F "cell centers" ought to default to the minimum TBOXLO of each forecast case. The F "cell bounds" are undefined, and I'd suggest just assigning these bounds to sit midway between the cell centers, using whatever algorithm Ferret usually relies on when it sees an axis with no "bounds" attribute. A user could then always regrid the F-axis to something else (e.g. the 1-month-wide strips you see above) using @bin.

In short, there are two well-defined and completely independent calendar axes: (1) the axis of start times F, defined using cell centers; and (2) the axis of verification times, defined as an inter-forecast union of T bounds (for time-averaged diagnostics like monthly-means) or T centers (for instantaneous snapshots). The "diagonal plot" at top-right (if the strips were infinitesimally thin in the vertical) would be the fundamental representation of the forecast data -- always well-defined and accurate.

Then for convenience, we also have a constructed, fudged time axis of "approximate elapsed times". This assumes that all forecast cases have "roughly similar" elapsed times (e.g. 1 month spacing), so that the "inter-forecast average" elapsed times will reasonably approximate any individual forecast's true elapsed times.

[karlmsmith]

Agreeing that the result of the forecast aggregation should be the "diagonal" form - see: https://github.com/NOAA-PMEL/Ferret/issues/1665#issuecomment-428775469 Seems like this change would make things easier and more accurate; not sure why Steve chose to go with the compact form.