Processing the OLCI `mask` dataset raises an error

[simonrp84]

Describe the bug The Sen-3 / OLCI L2 data files contain a mask variable that provides quality information. When you try to save this dataset an error is raised complaining about bitwise operations on boolean arrays. I suspect this is because we're trying to save a boolean array, which is not supported somehow.

Can we update satpy somehow to work with bool arrays? My guess is that the geotiff writer would need some kind of if statement to catch the bool dtype and then transform it into a zero/one array or something like that. Unfortunately I'm not familiar enough with the writer code to work it out easily though...

To Reproduce

from satpy import Scene
from glob import glob
scn = Scene(glob('./*/*.nc'), reader='olci_l2')
scn.load(['mask']);
scn.save_dataset('mask')

Expected behavior The dataset to be saved.

Actual results This, rather large, error message:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[1], line 5
      3 scn = Scene(glob('./*/*.nc'), reader='olci_l2')
      4 scn.load(['mask']);
----> 5 scn.save_dataset('mask')

File ~\PycharmProjects\satpy\satpy\scene.py:1221, in Scene.save_dataset(self, dataset_id, filename, writer, overlay, decorate, compute, **kwargs)
   1216     writer = self._get_writer_by_ext(os.path.splitext(filename)[1])
   1218 writer, save_kwargs = load_writer(writer,
   1219                                   filename=filename,
   1220                                   **kwargs)
-> 1221 return writer.save_dataset(self[dataset_id],
   1222                            overlay=overlay, decorate=decorate,
   1223                            compute=compute, **save_kwargs)

File ~\PycharmProjects\satpy\satpy\writers\__init__.py:885, in ImageWriter.save_dataset(self, dataset, filename, fill_value, overlay, decorate, compute, units, **kwargs)
    882     dataset = dataset.pint.quantify().pint.to(units).pint.dequantify()
    883 img = get_enhanced_image(dataset.squeeze(), enhance=self.enhancer, overlay=overlay,
    884                          decorate=decorate, fill_value=fill_value)
--> 885 return self.save_image(img, filename=filename, compute=compute, fill_value=fill_value, **kwargs)

File ~\PycharmProjects\satpy\satpy\writers\geotiff.py:278, in GeoTIFFWriter.save_image(self, img, filename, compute, dtype, fill_value, keep_palette, cmap, tags, overviews, overviews_minsize, overviews_resampling, include_scale_offset, scale_offset_tags, colormap_tag, driver, tiled, **kwargs)
    275     tags = {}
    276 tags.update(self.tags)
--> 278 return img.save(filename, fformat='tif', driver=driver,
    279                 fill_value=fill_value,
    280                 dtype=dtype, compute=compute,
    281                 keep_palette=keep_palette, cmap=cmap,
    282                 tags=tags, include_scale_offset_tags=include_scale_offset,
    283                 scale_offset_tags=scale_offset_tags,
    284                 colormap_tag=colormap_tag,
    285                 overviews=overviews,
    286                 overviews_resampling=overviews_resampling,
    287                 overviews_minsize=overviews_minsize,
    288                 tiled=tiled,
    289                 **gdal_options)

File ~\miniconda3\Lib\site-packages\trollimage\xrimage.py:248, in XRImage.save(self, filename, fformat, fill_value, compute, keep_palette, cmap, driver, **format_kwargs)
    246 if fformat in ('tif', 'tiff', 'jp2'):
    247     try:
--> 248         return self.rio_save(filename, fformat=fformat, driver=driver,
    249                              fill_value=fill_value, compute=compute,
    250                              keep_palette=keep_palette, cmap=cmap,
    251                              **format_kwargs)
    252     except ImportError:
    253         logger.warning("Missing 'rasterio' dependency to save GeoTIFF "
    254                        "image. Will try using PIL...")

File ~\miniconda3\Lib\site-packages\trollimage\xrimage.py:432, in XRImage.rio_save(self, filename, fformat, fill_value, dtype, compute, tags, keep_palette, cmap, overviews, overviews_minsize, overviews_resampling, include_scale_offset_tags, scale_offset_tags, colormap_tag, driver, **format_kwargs)
    428     to_store = list(zip(*([to_store] + da_tags)))
    430 if compute:
    431     # write data to the file now
--> 432     res = da.store(*to_store)
    433     to_close = to_store[1]
    434     if not isinstance(to_close, tuple):

File ~\miniconda3\Lib\site-packages\dask\threaded.py:89, in get(dsk, keys, cache, num_workers, pool, **kwargs)
     86     elif isinstance(pool, multiprocessing.pool.Pool):
     87         pool = MultiprocessingPoolExecutor(pool)
---> 89 results = get_async(
     90     pool.submit,
     91     pool._max_workers,
     92     dsk,
     93     keys,
     94     cache=cache,
     95     get_id=_thread_get_id,
     96     pack_exception=pack_exception,
     97     **kwargs,
     98 )
    100 # Cleanup pools associated to dead threads
    101 with pools_lock:

File ~\miniconda3\Lib\site-packages\dask\local.py:511, in get_async(submit, num_workers, dsk, result, cache, get_id, rerun_exceptions_locally, pack_exception, raise_exception, callbacks, dumps, loads, chunksize, **kwargs)
    509         _execute_task(task, data)  # Re-execute locally
    510     else:
--> 511         raise_exception(exc, tb)
    512 res, worker_id = loads(res_info)
    513 state["cache"][key] = res

File ~\miniconda3\Lib\site-packages\dask\local.py:319, in reraise(exc, tb)
    317 if exc.__traceback__ is not tb:
    318     raise exc.with_traceback(tb)
--> 319 raise exc

File ~\miniconda3\Lib\site-packages\dask\local.py:224, in execute_task(key, task_info, dumps, loads, get_id, pack_exception)
    222 try:
    223     task, data = loads(task_info)
--> 224     result = _execute_task(task, data)
    225     id = get_id()
    226     result = dumps((result, id))

File ~\miniconda3\Lib\site-packages\trollimage\xrimage.py:1010, in XRImage._compute_quantile(data, dims, cutoffs)
   1008 data_arr = xr.DataArray(data, dims=dims)
   1009 # delayed will provide us the fully computed xarray with ndarray
-> 1010 left, right = data_arr.quantile([cutoffs[0], 1. - cutoffs[1]], dim=['x', 'y'])
   1011 logger.debug("Interval: left=%s, right=%s", str(left), str(right))
   1012 return left.data, right.data

File ~\miniconda3\Lib\site-packages\xarray\core\dataarray.py:5093, in DataArray.quantile(self, q, dim, method, keep_attrs, skipna, interpolation)
   4983 def quantile(
   4984     self: T_DataArray,
   4985     q: ArrayLike,
   (...)
   4990     interpolation: QuantileMethods | None = None,
   4991 ) -> T_DataArray:
   4992     """Compute the qth quantile of the data along the specified dimension.
   4993 
   4994     Returns the qth quantiles(s) of the array elements.
   (...)
   5090        The American Statistician, 50(4), pp. 361-365, 1996
   5091     """
-> 5093     ds = self._to_temp_dataset().quantile(
   5094         q,
   5095         dim=dim,
   5096         keep_attrs=keep_attrs,
   5097         method=method,
   5098         skipna=skipna,
   5099         interpolation=interpolation,
   5100     )
   5101     return self._from_temp_dataset(ds)

File ~\miniconda3\Lib\site-packages\xarray\core\dataset.py:7656, in Dataset.quantile(self, q, dim, method, numeric_only, keep_attrs, skipna, interpolation)
   7650     if name not in self.coords:
   7651         if (
   7652             not numeric_only
   7653             or np.issubdtype(var.dtype, np.number)
   7654             or var.dtype == np.bool_
   7655         ):
-> 7656             variables[name] = var.quantile(
   7657                 q,
   7658                 dim=reduce_dims,
   7659                 method=method,
   7660                 keep_attrs=keep_attrs,
   7661                 skipna=skipna,
   7662             )
   7664 else:
   7665     variables[name] = var

File ~\miniconda3\Lib\site-packages\xarray\core\variable.py:2311, in Variable.quantile(self, q, dim, method, keep_attrs, skipna, interpolation)
   2306         raise ValueError(
   2307             f"Interpolation method '{method}' requires numpy >= 1.22 or is not supported."
   2308         )
   2309     kwargs = {"q": q, "axis": axis, "interpolation": method}
-> 2311 result = apply_ufunc(
   2312     _wrapper,
   2313     self,
   2314     input_core_dims=[dim],
   2315     exclude_dims=set(dim),
   2316     output_core_dims=[["quantile"]],
   2317     output_dtypes=[np.float64],
   2318     dask_gufunc_kwargs=dict(output_sizes={"quantile": len(q)}),
   2319     dask="parallelized",
   2320     kwargs=kwargs,
   2321 )
   2323 # for backward compatibility
   2324 result = result.transpose("quantile", ...)

File ~\miniconda3\Lib\site-packages\xarray\core\computation.py:1207, in apply_ufunc(func, input_core_dims, output_core_dims, exclude_dims, vectorize, join, dataset_join, dataset_fill_value, keep_attrs, kwargs, dask, output_dtypes, output_sizes, meta, dask_gufunc_kwargs, *args)
   1205 # feed Variables directly through apply_variable_ufunc
   1206 elif any(isinstance(a, Variable) for a in args):
-> 1207     return variables_vfunc(*args)
   1208 else:
   1209     # feed anything else through apply_array_ufunc
   1210     return apply_array_ufunc(func, *args, dask=dask)

File ~\miniconda3\Lib\site-packages\xarray\core\computation.py:761, in apply_variable_ufunc(func, signature, exclude_dims, dask, output_dtypes, vectorize, keep_attrs, dask_gufunc_kwargs, *args)
    756     if vectorize:
    757         func = _vectorize(
    758             func, signature, output_dtypes=output_dtypes, exclude_dims=exclude_dims
    759         )
--> 761 result_data = func(*input_data)
    763 if signature.num_outputs == 1:
    764     result_data = (result_data,)

File ~\miniconda3\Lib\site-packages\xarray\core\variable.py:2298, in Variable.quantile.<locals>._wrapper(npa, **kwargs)
   2296 def _wrapper(npa, **kwargs):
   2297     # move quantile axis to end. required for apply_ufunc
-> 2298     return np.moveaxis(_quantile_func(npa, **kwargs), 0, -1)

File <__array_function__ internals>:200, in quantile(*args, **kwargs)

File ~\miniconda3\Lib\site-packages\numpy\lib\function_base.py:4461, in quantile(a, q, axis, out, overwrite_input, method, keepdims, interpolation)
   4459 if not _quantile_is_valid(q):
   4460     raise ValueError("Quantiles must be in the range [0, 1]")
-> 4461 return _quantile_unchecked(
   4462     a, q, axis, out, overwrite_input, method, keepdims)

File ~\miniconda3\Lib\site-packages\numpy\lib\function_base.py:4473, in _quantile_unchecked(a, q, axis, out, overwrite_input, method, keepdims)
   4465 def _quantile_unchecked(a,
   4466                         q,
   4467                         axis=None,
   (...)
   4470                         method="linear",
   4471                         keepdims=False):
   4472     """Assumes that q is in [0, 1], and is an ndarray"""
-> 4473     return _ureduce(a,
   4474                     func=_quantile_ureduce_func,
   4475                     q=q,
   4476                     keepdims=keepdims,
   4477                     axis=axis,
   4478                     out=out,
   4479                     overwrite_input=overwrite_input,
   4480                     method=method)

File ~\miniconda3\Lib\site-packages\numpy\lib\function_base.py:3752, in _ureduce(a, func, keepdims, **kwargs)
   3749             index_out = (0, ) * nd
   3750             kwargs['out'] = out[(Ellipsis, ) + index_out]
-> 3752 r = func(a, **kwargs)
   3754 if out is not None:
   3755     return out

File ~\miniconda3\Lib\site-packages\numpy\lib\function_base.py:4639, in _quantile_ureduce_func(a, q, axis, out, overwrite_input, method)
   4637     else:
   4638         arr = a.copy()
-> 4639 result = _quantile(arr,
   4640                    quantiles=q,
   4641                    axis=axis,
   4642                    method=method,
   4643                    out=out)
   4644 return result

File ~\miniconda3\Lib\site-packages\numpy\lib\function_base.py:4756, in _quantile(arr, quantiles, axis, method, out)
   4754     result_shape = virtual_indexes.shape + (1,) * (arr.ndim - 1)
   4755     gamma = gamma.reshape(result_shape)
-> 4756     result = _lerp(previous,
   4757                    next,
   4758                    gamma,
   4759                    out=out)
   4760 if np.any(slices_having_nans):
   4761     if result.ndim == 0 and out is None:
   4762         # can't write to a scalar

File ~\miniconda3\Lib\site-packages\numpy\lib\function_base.py:4573, in _lerp(a, b, t, out)
   4559 def _lerp(a, b, t, out=None):
   4560     """
   4561     Compute the linear interpolation weighted by gamma on each point of
   4562     two same shape array.
   (...)
   4571         Output array.
   4572     """
-> 4573     diff_b_a = subtract(b, a)
   4574     # asanyarray is a stop-gap until gh-13105
   4575     lerp_interpolation = asanyarray(add(a, diff_b_a * t, out=out))

TypeError: numpy boolean subtract, the `-` operator, is not supported, use the bitwise_xor, the `^` operator, or the logical_xor function instead.

Environment Info:

• OS: Windows 11
• Satpy Version: Latest github main

[djhoese]

You need to define an enhancement that is a no-op. So operations: [] and then it might work. At least there wouldn't be any scaling. There might be conditions in trollimage that check for np.issubdtype(X.dtype, np.integer) that, if that condition is False for booleans, should maybe be not np.issubdtype(X.dtype, np.floating) instead.