Update And process to cope with array and a number

[danielFlemstrom]

Either I have misunderstood the openeo cookbook, or we need to do some work on and_.
I would expect mask =( band & 2) = 2 to work.
In the generated code, the operator ´=´becomes EQ which can handle array:int, but the operator ´&´ becomes logic and_ which assumes that both are arrays as I can see. Although & and and in python seems to give the same result, i suspect things get worse when handling arrays.

WATER_FLAG_MASK=2
res=connection.load_collection(s3.s3_olci_l2wfr,
                                 spatial_extent=s3.bbox.karlstad_large,
                                 bands=['oa04','oa06','oa09','wqsf_lsb'],
                                 temporal_extent=["2022-05-01", "2022-05-02"])
                                
wqsf_band = res.band('wqsf_lsb')  
mask = (wqsf_band & WATER_FLAG_MASK ) == WATER_FLAG_MASK  
res = res.mask(mask)
display(res.flat_graph())

-->

{'loadcollection1': {...},
 'reducedimension1': {'process_id': 'reduce_dimension',
  'arguments': {'data': {'from_node': 'loadcollection1'},
   'dimension': 'bands',
   'reducer': {'process_graph': {'arrayelement1': {'process_id': 'array_element',
      'arguments': {'data': {'from_parameter': 'data'}, 'index': 3}},
     'and1': {'process_id': 'and',
      'arguments': {'x': {'from_node': 'arrayelement1'}, 'y': 2}},
     'eq1': {'process_id': 'eq',
      'arguments': {'x': {'from_node': 'and1'}, 'y': 2},
      'result': True}}}}},
 'mask1': {'process_id': 'mask',
  'arguments': {'data': {'from_node': 'loadcollection1'},
   'mask': {'from_node': 'reducedimension1'}},
  'result': True}}

Which generates the following script:

...
_arrayelement1_2 = oeop.array_element(
    **{"data": _loadcollection1_0, "index": 3, "dimension": "bands"}
)
_and1_3 = oeop.and_(**{"x": _arrayelement1_2, "y": 2})
_eq1_4 = oeop.eq(**{"x": _and1_3, "y": 2})
_reducedimension1_1 = oeop.reduce_dimension(
    **{"data": _eq1_4, "dimension": "bands", "reducer": {}}
)
_mask1_5 = oeop.mask(**{"data": _loadcollection1_0, "mask": _reducedimension1_1})
...

... which results in some error text:

Traceback (most recent call last):
  File "code.py", line 69, in <module>
    _and1_3 = oeop.and_(**{"x": _arrayelement1_2, "y": 2})
  File "/user/ubuntu/.venv/lib/python3.8/site-packages/openeo_processes/utils.py", line 116, in fun_wrapper
    return cls_fun(*args, **kwargs)
  File "/user/ubuntu/.venv/lib/python3.8/site-packages/openeo_processes/logic.py", line 96, in exec_xar
    y_nan = y.where(y == True, False)
AttributeError: 'int' object has no attribute 'where'

/All the best
The DigitalEarthSweden team!

[danielFlemstrom]

Thinking this over, logical and between 2 and an array is a bit weird. Maybe we should add a bit-wise version, or have I only missed this in the implementation? The client should emit something else than logical and for the & operator perhaps...

[clausmichele]

Hi @danielFlemstrom. Looking at what you're doing using an and operator with integer values doesn't indeed make much sense. If I understood correctly, couldn't you do something like this?

mask = (wqsf_band == WATER_FLAG_MASK ) * 2

You firstly get a boolean mask, where your data is equal to 2 and then multiply it with 2 so that the result is either 0 or 2.

[danielFlemstrom]

The wqsf band contains more info, so each bit means something like this

and more than one bit may be true

Take two was to go through each pixel in the mask and mask out the first bit like this:

(tried to follow the cook-book to work with individual pixels).
Obviously this does not work either, but maybe it makes it clearer what I want to do
I thought x would be a pixel, but it is some array? There is something here that i have missed....

[clausmichele]

and more than one bit may be true

I would not know how to handle that case actually! The problem is that you from opendatacube/xarray side you don't see the individual bits, but integer values that are represented by the combination of those bit flags.

Therefore, if we consider the case where only bit 2 is set to True/1, you can create a boolean mask comparing the values with 2 using the eq operator.

WATER_FLAG_MASK=2
res=connection.load_collection(s3.s3_olci_l2wfr,
                                 spatial_extent=s3.bbox.karlstad_large,
                                 bands=['oa04','oa06','oa09','wqsf_lsb'],
                                 temporal_extent=["2022-05-01", "2022-05-02"])
                                
wqsf_band = res.band('wqsf_lsb')  
mask = (wqsf_band == WATER_FLAG_MASK  )
res = res.mask(mask)

[danielFlemstrom]

Each pixel in the wqsf band corresponds to a pixel int the "image". And if the mask is 3 pixels like this:
100010 == 2 False
000010 == 2 True
010010 == 2 False
So the water flag is set in all three, but only in the second case it would be added to the final mask.
Something like this

It seems like both apply and &/and gets an integer and an array to process. I suggest I setup my environment and try to add and_ between an integer and an array, just as the np.bitwise does. I assume xr array can do the same. This way we could do either of the solutions above. At least sentinel3 images have such flag fields, so I think this would be useful in the long run.

[clausmichele]

Yes I see your point! @m-mohr this is currently not possible right? We would need a new process to perform a bitwise and operation like numpy does:

Computes the bit-wise AND of the underlying binary representation of the integers in the input arrays. This ufunc implements the C/Python operator &.

https://numpy.org/doc/stable/reference/generated/numpy.bitwise_and.html

[danielFlemstrom]

Well its all a mess.. I have also seen this:
https://numpy.org/doc/stable/reference/generated/numpy.logical_and.html so the logical and bitwise is "the same" but not the same in python.... crystal clear right :)

BUT, I am working with an implementation, to keep things simple, as i understand it:

CASE 1: (both ways obviously)
x = int, y = (xr or np)  => nan if e is nan else x & e for all e in y.
CASE 2:
x = (xr or np), y = (xr or np)  => nan if ex or ey is nan else ex & ey for all ex,xy in x,y. 

If you have True or False, you get True or False back, if you have 1 and 0, you get 1 or 0 back.....
I'll do some experimenting and get back to you..

[m-mohr]

Yes, indeed. bitwise operations are not yet defined.
What you could probably do is to convert them to normal integers during the load process. The STAC metadata can help to describe this.

[clausmichele]

@m-mohr those values are indeed visible as normal integers at the back-end side (when loaded from opendatacube for instance).
However, I would not know how to easily get the same result as the bitwise and operation using integers and other mathematical operations (and expressing it with openEO processes)

[m-mohr]

What are you trying to achieve? I couldn't really follow the discussion...

[danielFlemstrom]

The mask-band from sentinel3 consists of several bits. Each bit means different things, e.g bit 1 means WATER and another bit indicates SNOW. So to create a mask for water, we need to mask out the other values (only keeping bit 1) and then make this a boolean mask. Only using ==we would miss out snowy water.

[m-mohr]

My thinking would be that creating these masks automatically and offering them to users might be beneficial and user-friendly. But it might also be expensive. Expose each of them as (virtual) bands and if requested, compute them on the fly.

The path with "less resistance" would probably be to define bitwise operations for openEO instead.

[danielFlemstrom]

Gr8! On it!
...But the current implementation differs for np and xr arrays: For numpy arrays, the result is already OK (the & operator is used), but for xArrays there is a something else going on with NaNs:

The "set NaN to false" discards anything that is not exactly 1. I got confused because I found no way of inserting np.nan into an array that can be used with the and operation. This NaN masking is the only thing hindering a bitwise-and as well...

I need some elaboration on the background here before braking the code:

1. what do you refer to when you say NaN ? Is it any negative number (you cannot and two floats) ?
   2 xArray IS an numpy array right? What is the difference to xArray (xArray & xArray is now the recommended way to go)
2. Adding np.nan would change the type of the array to float right ?

If we could consolidate numpy and xArray implementation, the current implementation would work "out of the box", however I suspect the xArray difference was not a random choice :D
Current Numpy impl:

[clausmichele]

it seems that xArray already wraps the numpy way of doing it:
https://docs.xarray.dev/en/stable/generated/xarray.ufuncs.logical_and.html

The only issue would be if the input is a float. So, either we automatically try to cast the floats to integers if needed, or throw an exception if the input is not boolean/integer.

import xarray as xr
import numpy as np

np.random.seed(0)
temperature = 15 + 8 * np.random.randn(2, 2)
lon = [[-99.83, -99.32], [-99.79, -99.23]]
lat = [[42.25, 42.21], [42.63, 42.59]]

da_float = xr.DataArray(
    data=temperature,
    dims=["x", "y"],
    coords=dict(
        lon=(["x", "y"], lon),
        lat=(["x", "y"], lat)
    )

)

print(da_float)

<xarray.DataArray (x: 2, y: 2)>
array([[29.11241877, 18.20125767],
       [22.82990387, 32.92714559]])
Coordinates:
    lon      (x, y) float64 -99.83 -99.32 -99.79 -99.23
    lat      (x, y) float64 42.25 42.21 42.63 42.59
Dimensions without coordinates: x, y

# bitwise/logical and with floats doesn't work:

da_float & da_float

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Input In [144], in <cell line: 2>()
      1 # bitwise/logical and with floats doesn't work:
----> 2 da_float & da_float

File ~/anaconda3/envs/eurac_env/lib/python3.8/site-packages/xarray/core/_typed_ops.py:227, in DataArrayOpsMixin.__and__(self, other)
    226 def __and__(self, other):
--> 227     return self._binary_op(other, operator.and_)

File ~/anaconda3/envs/eurac_env/lib/python3.8/site-packages/xarray/core/dataarray.py:3098, in DataArray._binary_op(self, other, f, reflexive)
   3094 other_variable = getattr(other, "variable", other)
   3095 other_coords = getattr(other, "coords", None)
   3097 variable = (
-> 3098     f(self.variable, other_variable)
   3099     if not reflexive
   3100     else f(other_variable, self.variable)
   3101 )
   3102 coords, indexes = self.coords._merge_raw(other_coords, reflexive)
   3103 name = self._result_name(other)

File ~/anaconda3/envs/eurac_env/lib/python3.8/site-packages/xarray/core/_typed_ops.py:417, in VariableOpsMixin.__and__(self, other)
    416 def __and__(self, other):
--> 417     return self._binary_op(other, operator.and_)

File ~/anaconda3/envs/eurac_env/lib/python3.8/site-packages/xarray/core/variable.py:2467, in Variable._binary_op(self, other, f, reflexive)
   2464 attrs = self._attrs if keep_attrs else None
   2465 with np.errstate(all="ignore"):
   2466     new_data = (
-> 2467         f(self_data, other_data) if not reflexive else f(other_data, self_data)
   2468     )
   2469 result = Variable(dims, new_data, attrs=attrs)
   2470 return result

TypeError: ufunc 'bitwise_and' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

temperature_int = [[18,20],[31,25]]
lon = [[-99.83, -99.32], [-99.79, -99.23]]
lat = [[42.25, 42.21], [42.63, 42.59]]

da_int = xr.DataArray(
    data=temperature_int,
    dims=["x", "y"],
    coords=dict(
        lon=(["x", "y"], lon),
        lat=(["x", "y"], lat)
    )
)

print(da_int & da_int)

<xarray.DataArray (x: 2, y: 2)>
array([[18, 20],
       [31, 25]])
Coordinates:
    lon      (x, y) float64 -99.83 -99.32 -99.79 -99.23
    lat      (x, y) float64 42.25 42.21 42.63 42.59
Dimensions without coordinates: x, y