CDAT / cdms

8 stars 10 forks source link

Change longitude range from -180--180 to 0--360 #371

Open Xunius opened 4 years ago

Xunius commented 4 years ago

Hi,

I noticed a change in behaviour from the old cdat (I was using cdat-lite) to the newer cdat8 in this line of code:

var=var(longitude=(0, 360))

The old version changes my data in -180 -- 180 longitude range to 0 -- 360, and also does the zonal shift, so the size remains the same.

The newer version clips the 0 -- 180 hemisphere out for me.

I wonder what's the new recommended way of zonal shifting and range change? Thanks.

jypeter commented 4 years ago

I have not noticed this suspicious behavior. Maybe there is something wrong or missing in the metadata of your longitude axis. You would probably get what you describe if your longitude was not recognized as a longitude axis (therefore cdms2 would not know about the 360 degrees periodicity)

Could you give the output of var.getLongitude(), and maybe the definition of the variable and the associated axes in the original netCDF file (part of the ncdump -h output), if var was loaded from an nc file

Xunius commented 4 years ago

Thanks for replying.

Here is the output from var.getLongitude():

>>> var.getLongitude()
   id: x
   Designated a longitude axis.
   units:  degree_east
   Length: 480
   First:  -180.0
   Last:   179.25
   Other axis attributes:
      modulo: 360.0
      realtopology: linear
      topology: circular
      axis: X
   Python id:  0x7efbe43b2090

Output from ncdump -h:

netcdf test {
dimensions:
    axis_23 = 1 ;
    y = 241 ;
    bound = 2 ;
    x = 480 ;
variables:
    int axis_23(axis_23) ;
        axis_23:units = "" ;
    double y(y) ;
        y:bounds = "bounds_y" ;
        y:axis = "Y" ;
        y:units = "degree_north" ;
        y:realtopology = "linear" ;
    double bounds_y(y, bound) ;
    float x(x) ;
        x:bounds = "bounds_x" ;
        x:axis = "X" ;
        x:modulo = 360. ;
        x:topology = "circular" ;
        x:units = "degree_east" ;
        x:realtopology = "linear" ;
    double bounds_x(x, bound) ;
    float tcw(axis_23, y, x) ;
        tcw:missing_value = 1.e+20f ;
        tcw:_FillValue = 1.e+20f ;

// global attributes:
        :Conventions = "CF-1.0" ;
}

This is an intermediate output file I wrote using cdat-lite so not so many metadata.

The code I used to reproduce my issue:

import cdms2 as cdms
fin=cdms.open('test.nc', 'r')
var=fin('tcw')
var2=var(longitude=(0, 360, 'co'))
print(var2.info())

In cdat-lite, var2 has 0--359.25 longitude range. In cdat8 has 0--179.25 longitude range.

jypeter commented 4 years ago

Your longitude axis has indeed a circular topology with modulo: 360.0, so it should work, unless maybe there is a problem with your actual longitude values?

I have just tried rewrapping on a test file in CDAT 8.1, and it works fine. You can try to do the same and see if it works on the test file

jypeter@obelix2 - ...jypeter - 43 >wget https://cdat.llnl.gov/cdat/sample_data/c                     lt.nc
--2019-10-17 10:55:09--  https://cdat.llnl.gov/cdat/sample_data/clt.nc
Resolving cdat.llnl.gov (cdat.llnl.gov)... 198.128.245.146
Connecting to cdat.llnl.gov (cdat.llnl.gov)|198.128.245.146|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1718908 (1.6M) [application/x-netcdf]
Saving to: 'clt.nc'

100%[======================================>] 1,718,908    661KB/s   in 2.5s

2019-10-17 10:55:12 (661 KB/s) - 'clt.nc' saved [1718908/1718908]

jypeter@obelix2 - ...jypeter - 45 >conda activate cdat-8.1_py3

(cdat-8.1_py3) jypeter@obelix2 - ...jypeter - 46 >ncdump -v longitude clt.nc
[...]
        float longitude(longitude) ;
                longitude:units = "degrees_east" ;
                longitude:long_name = "Longitude" ;

[...]
 longitude = -180, -175, -170, -165, -160, -155, -150, -145, -140, -135,
    -130, -125, -120, -115, -110, -105, -100, -95, -90, -85, -80, -75, -70,
    -65, -60, -55, -50, -45, -40, -35, -30, -25, -20, -15, -10, -5, 0, 5, 10,
    15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100,
    105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165, 170, 175 ;
}

(cdat-8.1_py3) jypeter@obelix2 - ...jypeter - 47 >python
Python 3.6.7 | packaged by conda-forge | (default, Feb 28 2019, 09:07:38)
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import cdms2
>>> f = cdms2.open('clt.nc')
>>> clt_orig = f('clt')
>>> clt_orig.getLongitude()[:]
array([-180., -175., -170., -165., -160., -155., -150., -145., -140.,
       -135., -130., -125., -120., -115., -110., -105., -100.,  -95.,
        -90.,  -85.,  -80.,  -75.,  -70.,  -65.,  -60.,  -55.,  -50.,
        -45.,  -40.,  -35.,  -30.,  -25.,  -20.,  -15.,  -10.,   -5.,
          0.,    5.,   10.,   15.,   20.,   25.,   30.,   35.,   40.,
         45.,   50.,   55.,   60.,   65.,   70.,   75.,   80.,   85.,
         90.,   95.,  100.,  105.,  110.,  115.,  120.,  125.,  130.,
        135.,  140.,  145.,  150.,  155.,  160.,  165.,  170.,  175.],
      dtype=float32)

>>> clt_rewrap = clt_orig(longitude=(0, 360, 'co'))

>>> clt_rewrap.getLongitude()[:]
array([  0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,
        55.,  60.,  65.,  70.,  75.,  80.,  85.,  90.,  95., 100., 105.,
       110., 115., 120., 125., 130., 135., 140., 145., 150., 155., 160.,
       165., 170., 175., 180., 185., 190., 195., 200., 205., 210., 215.,
       220., 225., 230., 235., 240., 245., 250., 255., 260., 265., 270.,
       275., 280., 285., 290., 295., 300., 305., 310., 315., 320., 325.,
       330., 335., 340., 345., 350., 355.], dtype=float32)

>>> clt_wrap_when_reading = f('clt', longitude=(0, 360, 'co'))

>>> clt_wrap_when_reading.getLongitude()[:]
array([  0.,   5.,  10.,  15.,  20.,  25.,  30.,  35.,  40.,  45.,  50.,
        55.,  60.,  65.,  70.,  75.,  80.,  85.,  90.,  95., 100., 105.,
       110., 115., 120., 125., 130., 135., 140., 145., 150., 155., 160.,
       165., 170., 175., 180., 185., 190., 195., 200., 205., 210., 215.,
       220., 225., 230., 235., 240., 245., 250., 255., 260., 265., 270.,
       275., 280., 285., 290., 295., 300., 305., 310., 315., 320., 325.,
       330., 335., 340., 345., 350., 355.], dtype=float32)
>>>
Xunius commented 4 years ago

@jypeter Thanks for the info. I tried the test file and it worked.

I also tried re-creating the longitude axis in cdat8, using the original values, and re-assigning the axis, then it worked. Really weird.

I've uploaded my nc file to a repo (https://github.com/Xunius/dummy_repo), would you like to have a look?