value scaling wrong in special cases

[forman]

For the same netCDF file used in #821, the value scaling seems to be wrongly applied to compute float64 surface temperature values from a (signed) short variable analysed_sst:

short analysed_sst(time=1, lat=3600, lon=7200);
  :_FillValue = -32768S; // short
  :units = "kelvin";
  :scale_factor = 0.01f; // float
  :add_offset = 273.15f; // float
  :long_name = "analysed sea surface temperature";
  :valid_min = -300S; // short
  :valid_max = 4500S; // short
  :standard_name = "sea_water_temperature";
  :depth = "20 cm";
  :source = "ATSR<1,2>-ESACCI-L3U-v1.0, AATSR-ESACCI-L3U-v1.0, AVHRR<12,14,15,16,17,18>_G-ESACCI-L2P-v1.0, AVHRRMTA-ESACCI-L2P-v1.0";
  :comment = "SST analysis produced for ESA SST CCI project using the OSTIA system in reanalysis mode.";
  :_ChunkSizes = 1, 1196, 2393; // int

Values are roughly -50 to 600 Kelvin instead of 270 to 310 Kelvin. It seems like the problem arises from misinterpreting the signed short raw values in the netCDF file.

Here is a notebook that better explains the issue: https://github.com/CCI-Tools/sandbox/blob/4c7a98a4efd1ba55152d2799b499cb27027c2b45/notebooks/norman/xarray-sst-issues.ipynb

[fmaussion]

It seems that the problem is happening at the NetCDF4 level already:

import netCDF4
import numpy as np

f = '/home/mowglie/Downloads/20100101120000-ESACCI-L4_GHRSST-SSTdepth-OSTIA-GLOB_LT-v02.0-fv01.1.nc'
d = netCDF4.Dataset(f)
da = d['analysed_sst'][:]
print(np.max(da), np.min(da))
print(netCDF4.num2date(d['time'][:], d['time'].units))

prints

600.82 -54.53
[datetime.datetime(1947, 5, 12, 9, 58, 14)]

[rabernat]

@fmaussion then it's a netCDF bug
https://github.com/Unidata/netcdf4-python/issues

[fmaussion]

For what its worth, I've tested your file with my IDL library, and I get:

IDL> print, min(data), max(data)
     -54.5300      306.440
IDL> print, TIME_to_STR(time)
01.01.2010 12:00:00

you should fill a report in netcdf4 as @rabernat suggests.

[shoyer]

So we actually use our own value scaling logic, independent of netCDF4. It's interesting that we have the same issue, though!

Possibly it's because we don't use the valid_min/valid_max arguments in determining what should be masked? I can download your dataset later and take a more in depth look.

[shoyer]

I just opened this up with decode_cf=False. Here's what the raw data looks like:

ds.analysed_sst[0].plot.imshow()

with attributes:

Attributes:
    _FillValue: -32768
    units: kelvin
    scale_factor: 0.01
    add_offset: 273.15
    long_name: analysed sea surface temperature
    valid_min: -300
    valid_max: 4500
    standard_name: sea_water_temperature
    depth: 20 cm
    source: ATSR<1,2>-ESACCI-L3U-v1.0, AATSR-ESACCI-L3U-v1.0, AVHRR<12,14,15,16,17,18>_G-ESACCI-L2P-v1.0, AVHRRMTA-ESACCI-L2P-v1.0
    comment: SST analysis produced for ESA SST CCI project using the OSTIA system in reanalysis mode.

To me, It looks like somebody just mis-prepared this dataset, given the smooth transitions in the ocean from dark red to dark blue, which would correspond to numeric overflow. If not, there are lots of places in the ocean where the temperature is in negative degrees Kelvin.

I don't have a strong opinion on whether or not we should automatically masking values outside valid_min/valid_max.

[forman]

Just found that the valid_min/valid_max attributes are not directly part of the CF convention but the NUG convention and applying to byte only according to CF Section 2.2 Data Types

All integer types are treated by the netCDF interface as signed. It is possible to treat the byte type as unsigned by using the NUG convention of indicating the unsigned range using the valid_min, valid_max, or valid_range attributes.

As for for #821, Panoply shows the correct values for the same file:

[fmaussion]

Note that ncview also displays the output correctly.

It seems that the problem occurs at the necdf4 level already:

import netCDF4
import matplotlib.pyplot as plt
f = '/home/mowglie/Downloads/20100101120000-ESACCI-L4_GHRSST-SSTdepth-OSTIA-GLOB_LT-v02.0-fv01.1.nc'
d = netCDF4.Dataset(f)
v = d['analysed_sst']
v.set_auto_maskandscale(False)
plt.imshow(v[0, ...], origin='lower')
plt.show()

To be compared to the equivalent output from IDL at the netcdf backend level:

[forman]

Ok, I'll submit a netCDF issue then.

[forman]

With h5py the data is read correctly too.

[shoyer]

@forman Just a note -- if h5py can read the data correctly, you can read the data into xarray using h5netcdf (with engine='h5netcdf' in open_dataset).

[forman]

Fantastic, thanks!

[stale]

In order to maintain a list of currently relevant issues, we mark issues as stale after a period of inactivity
If this issue remains relevant, please comment here; otherwise it will be marked as closed automatically

[dcherian]

Fixed upstream.