Handle grids where axis is reversed

[karlmsmith]

Reported by @AnsleyManke on 20 Sep 2013 18:23 UTC Our GFDL friends wrote us with the question quoted below, about Ferret's handling of axes with decreasing coordinates. We will be discussing how to resolve this in general.

An element of this which we can immediately address is to have Ferret WRITE data and axes in the same reversed order as in the dataset. Any variable defined on a grid which includes an axis that Ferret marked as reversed should be written with the coordinates and data reversed. In fact there is a NetCDF call to write data with a mapping vector for ordering the data, just as there is a call to read data in a permuted order.

I think we should do this by default, perhaps issuing a NOTE. This goes along with our efforts in recent Ferret releases to write netCDF files that conform to the original dataset as much as possible.

(from Matt, @AndrewWittenberg, Balaji, Keith, Alistair Adcroft 9/17/2013)

Hi Ansley,

We have identified an issue with Ferret's handling of vertical 
coordinates containing negative values.  The code below provides an 
example.

This issue came up in the development of  MOM6 diagnostics where we 
are now using height as opposed to depth as the vertical coordinate. 
According to the current CF documentation, this is allowed by the 
standard. Unfortunately, we hit on an oddity in ferret which, if we 
don't address, will cause all sorts of confusion for ferret users 
using our output.  Although plotting works, the data seems to be 
read in upside down into memory so that k=1 is at the bottom instead 
of the surface where z=0. Users get the expected result when 
referring to z=0 but not when referring to k=1. Is this a known 
issue and is there a workaround? or do we need to patch Ferret? 

Thanks,

Matt (and Alistair)

The following steps illustrate the problem:
cat <<EOFA >! test.cdl
netcdf test {
dimensions:
zl = 5 ;
variables:
double zl(zl) ;
zl:long_name = "Height Z" ;
zl:units = "m" ;
zl:positive = "up" ;
zl:cartesian_axis = "Z" ;
double temp(zl) ;
temp:long_name = "Temperature" ;
temp:units = "degrees C" ;
// global attributes:
:filename = "test.nc" ;
data:
zl = -5, -15, -25, -35, -45 ;
temp = 20, 16, 13, 9, 8 ;
}
EOFA
ncgen -b test.cdl
ferret
use test.nc
plot temp
list temp[k=1]
list temp[k=5]

You get a plot that IS correct but the data evidently has been inverted:
yes? list temp[k=1]
             VARIABLE : Temperature (degrees C)
             FILENAME : test.nc
             Z (m)    : -45
          8.000
yes? list temp[k=5]
             VARIABLE : Temperature (degrees C)
             FILENAME : test.nc
             Z (m)    : -5
          20.00

Migrated-From: http://dunkel.pmel.noaa.gov/trac/ferret/ticket/2093

[karlmsmith]

Modified by @AnsleyManke on 26 Sep 2013 17:57 UTC

[karlmsmith]

Comment by @AnsleyManke on 26 Sep 2013 18:03 UTC Here is a proposed set of changes to complete the work, changing the interpretation of indices on reverse-ordered axes. Continuing to support the reversed axes as in-memory, forward-ordered variables, but making this transparent to users:

• get the user's "/K=lo:hi" to be interpreted in a manner consistent with the reverse ordered axis *FLESH_OUT_AXIS seems the best place to do this. It is at this point that there is a guarantee that both the grid and the user's K= (or other index) request are known

  *add new variable to xcontext:   cx_reversal_has_been_applied(maxcx, maxferrdims)   (the long name is for clarity here)
  
  * IF    ( axis_is_reversed
  AND .NOT.cx_reversal_has_been_applied(cx, idim)
  AND   cx_by_ss(idim,cx)  ) THEN
        cx_lo_ss(cx,idim) = axlen-cx_lo_ss(cx,idim)+1
        cx_hi_ss(cx,idim) = axlen-cx_hi_ss(cx,idim)+1
        cx_reversal_has_been_applied(cx, idim) = TRUE
  ENDIF

• correct outputs back to the native axis ordering LIST (netCDF files and ASCII tables) SHOW AXIS, SHOW GRID
• New qualifier DEFINE AXIS/REVERSED (outputs from such axes will be reverse-ordered) needed to enable the normal range of control for reverse-ordered axes
• transformations such as RSUM and IIN that are direction-dependent ==> a high value solution to this: Ferret has always lacked versions of RSUM and IIN that worked in reverse order. We could implement new transformations @RRSUM and @RIIN that work in the reverse order. Then the support of reversed axes becomes mostly a matter of swapping between the "R" transforms and their forward-ordered counterparts (keeping the title comments straight as this is done)
• pseudo-variables : the index pseudo-variables, I, J, K, L, ... need to have the axlen-i+1 transformation applied. Ditto for some of the RETURN= immediate mode results.