Supprising type of nc.variables[foo][:]. (ndarray or MaskedArray depending on its data)

[MeraX]

Hi all,

for me it is surprising, that the return of nc.variables[foo][:] is either a ndarray or a MaskedArray depending on its data content: MaskedArrays are used as long as there is any "missing" entry in the data, ndarray otherwise. Even more surprising this data depending rules also apply to slices:

import netCDF4
import numpy as np

print(netCDF4.__version__)

with netCDF4.Dataset('test.nc', 'w') as nc:
    dim = 'test'
    nc.createDimension(dim, size=10)

    var = nc.createVariable('test-var', np.int,
        dim, fill_value=5)
    var[:] = np.arange(10)
    print('var[:]', var[:], type(var[:]))
    print('var[:5]', var[:5], type(var[:5]))
    print('var[:6]', var[:6], type(var[:6]))

output:

1.3.1
var[:] [0 1 2 3 4 -- -- -- -- --] <class 'numpy.ma.core.MaskedArray'>
var[:5] [0 1 2 3 4] <class 'numpy.ndarray'>
var[:6] [0 1 2 3 4 --] <class 'numpy.ma.core.MaskedArray'>

On the first glance it might seem convenient to get a MaskedArray, only if there is any missing data, but it confuses me very much, that one slice has another type (and behavior) as the other. And it is the same for the whole [:]: Imagine one day-file with missing entries, but another day without any.

If one likes to work with MaskedArrays, this means that one always have to check for the type and add a full False mask in case of no missing entries.

Is there any chance to make the return of var[:] more consistent and less surprising?

[MeraX]

A) One could suggest always returning MaskedArrays, as long as the auto-masking is turned on.

But this would introduce MaskedArray into code, that was never aware of having any missing values in its input data. And this would lead to very undesired results.
Before: read data does not contain any missing data, so it was returned as an array called old.
Now: Read data would be returned as MaskedArray called new with full False mask.
Imagine this old code printing [2]:

old = np.arange(3)
old_log = np.log(old-1)
print(old[np.isfinite(old_log)])

But with a MaskedArray, the result would be [0, 1, 2]:

new = np.ma.MaskedArray(np.arange(3))
new_log = np.log(new-1)
print(new[np.isfinite(new_log)])

B) Or one could always return MaskedArray, as long as the auto-masking is turned on and turn auto-masking off by default.

C) Or as an extension of B): the auto-masking is turned off for variable which have no specific FillValue (or valid_range etc.) set in the netcdf file. I assume, that a file author expecting no missing values, does normally not set a FillValue. Therefore, programs reading these files do not expect any missing data. But authors aware of missing data in their files, probably set an explicit FillValue.

I think, I would prefer B or C. Personally I would go for B, because I don't like working with MaskedArrays, but I can imagine B breaking code, that depends of reading MaskedArrays (so having at least one missing value).

I'm interested in how others think about this issue.

[jswhit]

In retrospect, I perhaps should not have made this the default behavior. I don't think it's wise to change it now though, since too much code relies on it. However, you can change it so that numpy arrays are always returned by using the set_auto_mask Dataset method. For example

nc = Dataset('test.nc')
nc.set_auto_mask(False)

will cause all data to be returned as numpy arrays. Additionally,

nc.set_auto_scale(False)

will suppress the auto-rescaling of data packed as short integers.

nc.set_auto_maskandscale(False)

will do both.

[jswhit]

Why is it so bad to get a masked array back when the slice contains missing data, and a numpy array when it does not? numpy.ma is pretty much a drop in replacement for numpy.ndarray - is there a particular situation where getting a masked array instead of a numpy array (or vice versa) causes a problem?

[akrherz]

Why is it so bad to get a masked array back when the slice contains missing data, and a numpy array when it does not?

I'm not the sharpest tool in the toolbox, but for me, this means that everywhere I read data from a netcdf variable, I have to do mask checking to see if I need to deal with masked data or not in my code logic. So I have to add code like

if isinstance(my_array, np.ma.MaskedArray):

[MeraX]

As I pointed out in my example, MaskedArrays are unfortunately not a drop-in replacement for arrays, at least in these edge cases, that rely on functions returning nan for invalid input. I could imagine, that np.log(-1) -> np.nan is undefined, but I checked the docs and it is actually noted

For real-valued input data types, log always returns real output. For
each value that cannot be expressed as a real number or infinity, it
yields nan and sets the invalid floating point error flag.

Off course I can use @akrherz if-case or always set_auto_mask(False), but I would prefer changing the default behavior of this pretty useful library to make it less surprising. (Actually I was surprised first, when I realized that netCDF4 started to apply the scale and offset after some update some yeas ago)

I don't think it's wise to change it now though, since too much code relies on it.

What about my option C? Turning auto_mask off if there is no explicit _FillValue set in the data.

Honestly, I have no idea, if any or how many netcdf files mark missing entries with the default FillValue. From my experience I don't know any, but this could off course be a misleading sample. So instead of using the default FillValue, if there is no explicit FillValue use no FillValue, which results in no detected missing entries, which results in returning an array in the old logic.

Second point, is there really code, that relays on getting an array for a complete slice and a MaskedArray for an incomplete slice?

I totally agree, that it is always risky to change default behavior because old code might relay exactly on this behavior. But it is good to have this discussion to elucidate all pros and cons of probable side effects and if they are real side effects.

[jswhit]

There is always a _FillValue - if it's not set explicitly, then the C library uses the type-dependent default value. If you create a variable and don't write any data to it, then read it back it will always be masked.

[jswhit]

The only potential option I see for changing the default behavior is to always return a masked array.

[MeraX]

There is always a _FillValue - if it's not set explicitly, then the C library uses the type-dependent default value. If you create a variable and don't write any data to it, then read it back it will always be masked.

Yes and no. In the current implementation, non filled entries are populated with the default fill value and are interpreted as such but the _FillValue attribute is not defined:

import netCDF4
import numpy as np

print(netCDF4.__version__)

with netCDF4.Dataset('test_without_fillvalue.nc', 'w') as nc:
    nc.createDimension('dim', size=10)
    var = nc.createVariable('test', 'f8', 'dim')
    print("without fill value `hasattr(var, '_FillValue')`:", hasattr(var, '_FillValue'))
    print("without fill value `nc.variables['test'][:]`:", nc.variables['test'][:])

with netCDF4.Dataset('test_without_fillvalue.nc', 'r') as nc:
    print("without fill value `hasattr(nc.variables['test'], '_FillValue')`:", hasattr(nc.variables['test'], '_FillValue'))
    print("without fill value `nc.variables['test'][:]`:", nc.variables['test'][:])


with netCDF4.Dataset('test_with_fillvalue.nc', 'w') as nc:
    nc.createDimension('dim', size=10)
    var = nc.createVariable('test', 'f8', 'dim', fill_value=1.)
    print("with fill value `hasattr(var, '_FillValue')`:", hasattr(var, '_FillValue'))
    print("with fill value `nc.variables['test'][:]`:", nc.variables['test'][:])

with netCDF4.Dataset('test_with_fillvalue.nc', 'r') as nc:
    print("with fill value `hasattr(nc.variables['test'], '_FillValue')`:", hasattr(nc.variables['test'], '_FillValue'))
    print("with fill value `nc.variables['test'][:]`:", nc.variables['test'][:])


with netCDF4.Dataset('test_with_default_fillvalue.nc', 'w') as nc:
    nc.createDimension('dim', size=10)
    var = nc.createVariable('test', 'f8', 'dim', fill_value=netCDF4.default_fillvals['f8'])
    print("with default fill value `hasattr(var, '_FillValue')`:", hasattr(var, '_FillValue'))
    print("with default fill value `nc.variables['test'][:]`:", nc.variables['test'][:])

Output:

1.3.1
without fill value `hasattr(var, '_FillValue')`: False
without fill value `nc.variables['test'][:]`: [-- -- -- -- -- -- -- -- -- --]
without fill value `hasattr(nc.variables['test'], '_FillValue')`: False
without fill value `nc.variables['test'][:]`: [-- -- -- -- -- -- -- -- -- --]
with fill value `hasattr(var, '_FillValue')`: True
with fill value `nc.variables['test'][:]`: [-- -- -- -- -- -- -- -- -- --]
with fill value `hasattr(nc.variables['test'], '_FillValue')`: True
with fill value `nc.variables['test'][:]`: [-- -- -- -- -- -- -- -- -- --]
with default fill value `hasattr(var, '_FillValue')`: True
with default fill value `nc.variables['test'][:]`: [-- -- -- -- -- -- -- -- -- --]

ncdump:

netcdf test_without_fillvalue {
dimensions:
        dim = 10 ;
variables:
        double test(dim) ;

// global attributes:
                :_NCProperties = "version=1|netcdflibversion=4.4.1.1|hdf5libversion=1.8.18" ;
data:

 test = _, _, _, _, _, _, _, _, _, _ ;
}

netcdf test_with_fillvalue {
dimensions:
        dim = 10 ;
variables:
        double test(dim) ;
                test:_FillValue = 1. ;

// global attributes:
                :_NCProperties = "version=1|netcdflibversion=4.4.1.1|hdf5libversion=1.8.18" ;
data:

 test = _, _, _, _, _, _, _, _, _, _ ;
}

netcdf test_with_default_fillvalue {
dimensions:
        dim = 10 ;
variables:
        double test(dim) ;
                test:_FillValue = 9.96920996838687e+36 ;

// global attributes:
                :_NCProperties = "version=1|netcdflibversion=4.4.1.1|hdf5libversion=1.8.18" ;
data:

 test = _, _, _, _, _, _, _, _, _, _ ;
}

The test_without_fillvalue does not contain an explicit fill value. This is tested with hasattr().
The test_with_fillvalue contains an explicit fill value as one would expect.
And test_with_default_fillvalue uses the default fill value explicitly, but in contrast to the first case, this value is explicitly saved to the file.

At least as far as I can test it, there is a difference between an implicit use of the default fill value and an explicit use.

[jswhit]

The point is that whether or not it is set explicitly, there is an expectation that values equal to the _FillValue are missing.

[MeraX]

But do people expect that values equal to the default fill value are actually missing? Do people at all expect to have missing values, if _FillValue is not set explicitly?

That's something that I don't know.

From my experience: If there are missing values, some _FillValue is set. Or people fill missing values with NaN but don't set an explicit FillValue. (In the latter case, the auto mask doesn't mask out these NaNs anyway, since they are different from the default fill value)

[jswhit]

yes and yes. _FillValue is meant for pre-filling as yet unwritten data. missing_value is meant to indicate invalid or missing data. See https://www.unidata.ucar.edu/software/netcdf/docs/BestPractices.html.

[jswhit]

If the need is to know exactly what kind of array the module will return, the I propose we always return a masked array by default. If you always want a numpy array, you can use set_auto_mask(False). The only change to the current behavior would then be that a masked array would be returned even if the slice does not contain any masked values.

[shoyer]

I am strongly in support of type stable interfaces. It should not be data dependent whether indexing returns a NumPy arrays or a masked array.

If you always want a numpy array, you can use set_auto_mask(False).

I like synchronizing this behavior with set_auto_mask().

However, it's also worth considering whether to change the default behavior to not do masking and always return NumPy arrays intead.