Preserving dtype of element properties (continuation of #1496)

[vincentcasseau]

Hello,

In opendrift/elements/elements.py, there are a few places where the variable dtype is not preserved. This happens whenever np.ones is called. I am wondering if this is by design and if so why that is.

For instance:

# Status must always be array
if not type(self.status) == np.ndarray:
    self.status = self.status*np.ones(self.lon.shape)

could become

# Status must always be array
if not isinstance(self.status, np.ndarray):
    dtype = self.variables['status']['dtype']
    self.status = self.status*np.ones(self.lon.shape, dtype=dtype)

This prevents self.status from being an np.float64.

I am attaching a file showing all instances (please diff against current version of this file).

elements.txt

Thanks

[knutfrode]

Hi,
Internally in OpenDrift, float is used for all variables for simplicity.
However, the data type is converted to the specified type in the output Xarray Dataset (self.result) and in the output netCDF file.

Do you experience any issues with the internal, temporary change of datatype?

[vincentcasseau]

Thank you, Knut.

One downside of your approach is that we cannot use any of the integer variables defined in Elements to perform integer array indexing in the models without first casting them to int.

---

However, the data type is converted to the specified type in the output Xarray Dataset (self.result) and in the output netCDF file.

This is what I am looking at right now. What I uploaded does not entirely fix my problem as status is still exported as a float64.

I ran your example examples/example_leeway.py. I exported the solution as

# Running model
ds = o.run(duration=timedelta(hours=48), time_step=1800, time_step_output=3600, outfile='leeway.nc')
exit()

Then I opened it as

import opendrift
oa = opendrift.open("leeway.nc")

I am getting

16:06:50 INFO    opendrift:82: Returning <class 'opendrift.models.leeway.Leeway'> object
>>> oa.result['status']
<xarray.DataArray 'status' (trajectory: 1000, time: 49)> Size: 392kB
array([[ 0.,  0.,  0., ..., nan, nan, nan],
       [ 0.,  0.,  0., ..., nan, nan, nan],
       [ 0.,  0.,  0., ..., nan, nan, nan],
       ...,
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ..., nan, nan, nan],
       [ 0.,  0.,  0., ..., nan, nan, nan]], shape=(1000, 49))
Coordinates:
  * trajectory  (trajectory) int32 4kB 0 1 2 3 4 5 6 ... 994 995 996 997 998 999
  * time        (time) datetime64[ns] 392B 2015-11-16 ... 2015-11-18
Attributes:
    valid_range:    [0 1]
    flag_values:    [0 1]
    flag_meanings:  active stranded
>>> oa.result['status'].dtype
dtype('float64')

In addition, in opendrift/models/basemodel/__init__.py, status is excluded from the loop that computes minval and maxval. minval and maxval used to be computed for status in earlier versions of OD, which was convenient -- no need for valid_range and status could be dealt with like any other variable.
Also here, you are using the variable dtype from self.result which I understand is not the real dtype. I don't think you are then re-casting minval and maxval later.

if final is True or buffer_is_full:
    logger.debug('Updating minval and maxval')
    # Update min and max values
    for varname, var in self.result.data_vars.items():
        if varname in ['status']: 
            continue
        minval = var.min(skipna=True).item()
        maxval = var.max(skipna=True).item()
        if 'minval' in var.attrs:
            minval = np.minimum(minval, var.minval)
            maxval = np.maximum(maxval, var.maxval)
        # Attributes shall have same datatype as in output file
        if 'dtype' in var.attrs:
            att_type = var.attrs['dtype']
        else:
            att_type = np.float32
        minval = att_type(minval)
        maxval = att_type(maxval)
        self.result[varname] = self.result[varname].assign_attrs({'minval': minval, 'maxval': maxval})

Thanks again!

PS: After writing my comment, I thought this may have to do with NaNs but even this 2-time-step simulation w/o NaNs is returning float64

ds = o.run(duration=timedelta(hours=1), time_step=1800, time_step_output=1800, outfile='leeway.nc')

[knutfrode]

I recall now that the reason for using float internally also for int-variables (e.g. status) is that integer dtypes do not support NaN, which is needed for e.g. unseeded elements.

But int-variables are still written as int to the output netCDF-file (i.e. respecting dtypes from module specification), but again converted to float when imported with opendrift.open (i.e. just like during/after a simulation).

I do not remember why status was excluded from the loop for minval and maxval, perhaps I thought it was unnecessary as we anyway include the CF standard attribute valid_range (which is "obligatory").
However, I see that all tests are passing also when again including minval/maxval for status, so if that helps you, I can put those back (which should then be identical to valid_range).

[vincentcasseau]

I recall now that the reason for using float internally also for int-variables (e.g. status) is that integer dtypes do not support NaN, which is needed for e.g. unseeded elements.

That makes sense, so it's perfect this way.

But int-variables are still written as int to the output netCDF-file (i.e. respecting dtypes from module specification), but again converted to float when imported with opendrift.open (i.e. just like during/after a simulation).

I can confirm that opening the dataset as follows works as intended.

import netCDF4
ds = netCDF4.Dataset("leeway.nc")
print(ds['status'].dtype)  # int32
print(ds['origin_marker'].dtype)  # int32

To get the dtype with opendrift.open:

oa.result[varname].minval.dtype
oa.result['status'].attrs['valid_range'].dtype

This answers all my questions, thank you very much!

when again including minval/maxval for status, so if that helps you, I can put those back (which should then be identical to valid_range).

If you don't see any benefits of having them, then please don't do it. This is just one extra conditional statement. I didn't know valid_range was mandatory.

[vincentcasseau]

Hi Knut,

It turns out there is still something that could be improved when filling the array that will then be saved to file (in export/io_netCDF.py). You have on line 53:

for varname in self.result.data_vars:
    if 'time' in self.result[varname].dims:
        var = self.outfile.variables[varname]
        var[:, numtimes:numtimes + self.result.sizes['time']] = self.result[varname]

self.result is always float and may contain NaNs as you said, while var can be of type int.
You are properly dealing with dtype and NaNs on lines 32-47 during file export but that occurs later. In this assignment, no fill value are given for NaNs during float to int casting. Thus, the resulting value in your int array may be different from encoding[varname]['_FillValue'].

You can catch that as follows

First, turn xarray warnings into errors

import warnings
warnings.filterwarnings("error", category=RuntimeWarning, module="xarray")

Then let the code crash (you need a simulation long enough for the outfile to be present on your disk, say > 200 Dt)

10:12:43 INFO     Traceback (most recent call last):                                                                                                    
                    File "/home/vince/opendrift/opendrift/models/basemodel/__init__.py", line 2099, in run                                   
                      self.state_to_buffer()  # Append status to history array                                                                          
                    File "/home/vince/opendrift/opendrift/models/basemodel/__init__.py", line 2318, in state_to_buffer                       
                      self.io_write_buffer()                                                                                                            
                    File "/home/vince/opendrift/opendrift/export/io_netcdf.py", line 60, in write_buffer                                     
                      var[:, numtimes:numtimes + self.result.sizes['time']] = self.result[varname]                                                      
                    File "src/netCDF4/_netCDF4.pyx", line 5571, in netCDF4._netCDF4.Variable.__setitem__                                                
                    File "/home/vince/miniconda3/envs/opendrift/lib/python3.10/site-packages/xarray/core/common.py", line 181, in __array__             
                      return np.array(self.values, dtype=dtype, copy=copy)                                                                              
                  RuntimeWarning: invalid value encountered in cast                                                                                     
                                                                                                                                                        
10:12:43 INFO     'The simulation stopped before requested end time was reached.'

Something along these lines will help get rid of the warning and guarantee the correct fill value for integer fields:

for varname in self.result.data_vars:
        if 'time' in self.result[varname].dims:
            var = self.outfile.variables[varname]
            if np.issubdtype(var.dtype, np.integer):
                fill_value = getattr(var, '_FillValue', np.iinfo(var.dtype).max)
                var[:, numtimes:numtimes + self.result.sizes['time']] = self.result[varname].fillna(fill_value)   # implicit casting to int
            else:
                var[:, numtimes:numtimes + self.result.sizes['time']] = self.result[varname]

Thanks