cdms2 generating wrong longitude axis topology when reading grads dat/ctl data ?

jypeter commented 3 years ago

@jasonb5 I don't have time to upload a test file (before next week), but maybe you can test this anyway

I have converted a grads dat/ctl file using cdms2

The ctl file has the following content

dset   ^topography-read.dat
options big_endian
undef  -1.0000E+20
title ECBILT orography
xdef  64 linear    0.00   5.62
ydef  32 levels
 -85.7606 -80.2688 -74.7445 -69.2130 -63.6786
 -58.1430 -52.6065 -47.0696 -41.5325 -35.9951
 -30.4576 -24.9199 -19.3822 -13.8445 -8.30670
 -2.76890 2.76890 8.30670 13.8445 19.3822
  24.9199 30.4576 35.9951 41.5325 47.0696
  52.6065 58.1430 63.6786 69.2130 74.7445
  80.2688 85.7606
zdef   1 linear    0.00   1.00
tdef    1 linear 1jan0001  1YR
vars 1
orog 1  99 orography ECBILT
endvars

Things worked fine, except that the generated lon axis in the generated file has conflicting topology/realtopology information that ends up with my lon axis losing its periodicity. That is, the generated file has longitudes from 0 to 360, but when I request longitudes from -180 to 180, I end up with a longitude axis that only goes from 0 to 180

>>> v = f('orog')
>>> v.shape
(32, 64)
>>> v2 = f('orog', longitude=(-180,180,'co'))
>>> v2.shape
(32, 33)
>>> v.getLongitude()
   id: longitude
   Designated a longitude axis.
   units:  degrees_east
   Length: 64
   First:  0.0
   Last:   354.06
   Other axis attributes:
      axis: X
      modulo: 360.0
      topology: circular
      realtopology: linear
   Python id:  0x2b10141ff2e0

>>> v2.getLongitude()
   id: longitude
   Designated a longitude axis.
   units:  degrees_east
   Length: 33
   First:  0.0
   Last:   179.84
   Other axis attributes:
      axis: X
      modulo: 360.0
      topology: circular
      realtopology: linear
   Python id:  0x2b1014251f70

Any idea?

jypeter commented 3 years ago

I have tried the following, in order to fix the faulty files:

ncatted -a realtopology,longitude,d,, ICE-6G-C/orog_linear.nc ICE-6G-C/orog.nc: completely remove the realtopology attribute. Did not work. I was still getting only the [0, 180] range of the data
ncatted -a realtopology,longitude,m,c,circular ICE-6G-C/orog_linear.nc ICE-6G-C/orog.nc: replace the value with circular. Success! I was able to get data in the [-180, 180] range

So I'm also wondering how this realtopology attribute is handled, even when it is missing

jypeter commented 3 years ago

I have uploaded a grads ctl and dat file to our LSCE server and you will find below a session reproducing the issue discussed in this thread. Hopefully someone can use this to check why cdms2 is not setting the longitude realtopology correctly when importing a ctl/dat file

 > conda list | grep cdms
cdms2                     3.1.5                    pypi_0    pypi
libcdms                   3.1.2              h981a4fd_113    conda-forge

> python
Python 3.8.8 | packaged by conda-forge | (default, Feb 20 2021, 16:22:27)
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import cdms2

>>> f_grads = cdms2.open('topography-read.ctl')
>>> orog = f_grads('orog')
>>> f_grads.close()

>>> orog.info()
*** Description of Slab orog ***
id: orog
shape: (1, 1, 32, 64)
filename:
missing_value: -1e+20
comments:
grid_name: <None>
grid_type: generic
time_statistic:
long_name:
units:
tileIndex: None
title: orography ECBILT
Grid has Python id 0x2ac757a40940.
Gridtype: generic
Grid shape: (32, 64)
Order: yx
** Dimension 1 **
   id: time
   Designated a time axis.
   units:
   Length: 1
   First:  0.0
   Last:   0.0
   Other axis attributes:
      axis: T
      calendar: gregorian
   Python id:  0x2ac75efceca0
** Dimension 2 **
   id: level
   Designated a level axis.
   units:  lev
   Length: 1
   First:  0.0
   Last:   0.0
   Other axis attributes:
      axis: Z
   Python id:  0x2ac7dfdc44c0
** Dimension 3 **
   id: latitude
   Designated a latitude axis.
   units:  degrees_north
   Length: 32
   First:  -85.7606
   Last:   85.7606
   Other axis attributes:
      axis: Y
      realtopology: linear
   Python id:  0x2ac7dfdc44f0
** Dimension 4 **
   id: longitude
   Designated a longitude axis.
   units:  degrees_east
   Length: 64
   First:  0.0
   Last:   354.06
   Other axis attributes:
      axis: X
      modulo: 360.0
      topology: circular
      realtopology: linear
   Python id:  0x2ac7dfdc4520
*** End of description for orog ***

>>> lon = orog.getLongitude()
>>> lon
   id: longitude
   Designated a longitude axis.
   units:  degrees_east
   Length: 64
   First:  0.0
   Last:   354.06
   Other axis attributes:
      axis: X
      modulo: 360.0
      topology: circular
      realtopology: linear
   Python id:  0x2ac7dfdc4520

>>> lon.topology
'circular'
>>> lon.realtopology
'linear'

>>> orog.shape
(1, 1, 32, 64)

>>> orog_shift = orog(longitude=(-180,180, 'co'))

>>> orog_shift.shape
(1, 1, 32, 33)

>>> orog_shift.getLongitude()
   id: longitude
   Designated a longitude axis.
   units:  degrees_east
   Length: 33
   First:  0.0
   Last:   179.84
   Other axis attributes:
      axis: X
      modulo: 360.0
      topology: circular
      realtopology: linear
   Python id:  0x2ac7def7dd00

>>> lon.realtopology = 'circular'
>>> orog_shift_bis = orog(longitude=(-180,180, 'co')) 

>>> orog_shift_bis.shape
(1, 1, 32, 64)
>>> orog_shift_bis.getLongitude()
   id: longitude
   Designated a longitude axis.
   units:  degrees_east
   Length: 64
   First:  -174.54
   Last:   179.83997
   Other axis attributes:
      axis: X
      modulo: 360.0
      topology: circular
      realtopology: circular
   Python id:  0x2ac7e0152580

jypeter commented 3 years ago

Oh, I think this realtopology problem is probably why @Xunius was getting #371

CDAT / cdms

cdms2 generating wrong longitude axis topology when reading grads dat/ctl data ? #433