Data processing using the apply Operator with OpenEO

Apply Operator

The apply operator employ a process on the datacube that calculates new pixel values for each pixel, based on n other pixels.

Let’s start again with the same sample data from the Sentinel-2 STAC Collection, applying the filters directly in the load_stac call, to reduce the amount of data.

import openeo
from openeo.local import LocalConnection
local_conn = LocalConnection('')

url = "https://earth-search.aws.element84.com/v1/collections/sentinel-2-l2a"
spatial_extent = {"west": 11.259613, "east": 11.406212, "south": 46.461019, "north": 46.522237}
temporal_extent = ['2022-07-10T00:00:00Z','2022-07-13T00:00:00Z']
bands = ["red","green","blue"]
properties = {"eo:cloud_cover": dict(lt=50)}
datacube = local_conn.load_stac(url=url,
                    spatial_extent=spatial_extent,
                    temporal_extent = temporal_extent,
                    bands=bands,
                    properties=properties
)

datacube.execute()

Did not load machine learning processes due to missing dependencies: Install them like this: `pip install openeo-processes-dask[implementations, ml]`
/home/749c4c89-469c-4d4a-9ce4-9feffea0504a/.local/lib/python3.11/site-packages/stackstac/prepare.py:408: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.
  times = pd.to_datetime(
/home/749c4c89-469c-4d4a-9ce4-9feffea0504a/.local/lib/python3.11/site-packages/stackstac/prepare.py:408: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.
  times = pd.to_datetime(

<xarray.DataArray 'stackstac-b814a604da3396a4e43965fd3693671a' (time: 1,
                                                                band: 3,
                                                                y: 713, x: 1145)>
dask.array<getitem, shape=(1, 3, 713, 1145), dtype=float64, chunksize=(1, 1, 606, 860), chunktype=numpy.ndarray>
Coordinates: (12/54)
  * time                                     (time) datetime64[ns] 2022-07-12...
    id                                       (time) <U24 'S2B_32TPS_20220712_...
  * band                                     (band) <U5 'red' 'green' 'blue'
  * x                                        (x) float64 6.733e+05 ... 6.848e+05
  * y                                        (y) float64 5.155e+06 ... 5.148e+06
    instruments                              <U3 'msi'
    ...                                       ...
    proj:shape                               object {10980}
    gsd                                      int64 10
    common_name                              (band) <U5 'red' 'green' 'blue'
    center_wavelength                        (band) float64 0.665 0.56 0.49
    full_width_half_max                      (band) float64 0.038 0.045 0.098
    epsg                                     int64 32632
Attributes:
    spec:        RasterSpec(epsg=32632, bounds=(600000.0, 5090220.0, 709800.0...
    crs:         epsg:32632
    transform:   | 10.00, 0.00, 600000.00|\n| 0.00,-10.00, 5200020.00|\n| 0.0...
    resolution:  10.0

xarray.DataArray

'stackstac-b814a604da3396a4e43965fd3693671a'

• time: 1
• band: 3
• y: 713
• x: 1145

• dask.array<chunksize=(1, 1, 606, 860), meta=np.ndarray>

  Array Chunk Bytes 18.69 MiB 3.98 MiB Shape (1, 3, 713, 1145) (1, 1, 606, 860) Dask graph 12 chunks in 4 graph layers Data type float64 numpy.ndarray

• Coordinates: (54)
  • time
    (time)
    datetime64[ns]
    2022-07-12T10:18:03.747000

    array(['2022-07-12T10:18:03.747000000'], dtype='datetime64[ns]')

  • id
    (time)
    <U24
    'S2B_32TPS_20220712_0_L2A'

    array(['S2B_32TPS_20220712_0_L2A'], dtype='<U24')

  • band
    (band)
    <U5
    'red' 'green' 'blue'

    array(['red', 'green', 'blue'], dtype='<U5')

  • x
    (x)
    float64
    6.733e+05 6.733e+05 ... 6.848e+05

    array([673320., 673330., 673340., ..., 684740., 684750., 684760.])

  • y
    (y)
    float64
    5.155e+06 5.155e+06 ... 5.148e+06

    array([5154880., 5154870., 5154860., ..., 5147780., 5147770., 5147760.])

  • instruments
    ()
    <U3
    'msi'

    array('msi', dtype='<U3')

  • s2:degraded_msi_data_percentage
    ()
    float64
    0.0106

    array(0.0106)

  • s2:datatake_id
    ()
    <U34
    'GS2B_20220712T100559_027933_N04...

    array('GS2B_20220712T100559_027933_N04.00', dtype='<U34')

  • s2:sequence
    ()
    <U1
    '0'

    array('0', dtype='<U1')

  • eo:cloud_cover
    ()
    float64
    21.61

    array(21.614978)

  • view:sun_azimuth
    ()
    float64
    146.9

    array(146.8735106)

  • s2:product_type
    ()
    <U7
    'S2MSI2A'

    array('S2MSI2A', dtype='<U7')

  • platform
    ()
    <U11
    'sentinel-2b'

    array('sentinel-2b', dtype='<U11')

  • s2:product_uri
    ()
    <U65
    'S2B_MSIL2A_20220712T100559_N040...

    array('S2B_MSIL2A_20220712T100559_N0400_R022_T32TPS_20220712T125540.SAFE',
          dtype='<U65')

  • s2:generation_time
    ()
    <U27
    '2022-07-12T12:55:40.000000Z'

    array('2022-07-12T12:55:40.000000Z', dtype='<U27')

  • s2:reflectance_conversion_factor
    ()
    float64
    0.9674

    array(0.96742114)

  • s2:unclassified_percentage
    ()
    float64
    0.3591

    array(0.35912)

  • earthsearch:boa_offset_applied
    ()
    bool
    True

    array(True)

  • s2:datastrip_id
    ()
    <U64
    'S2B_OPER_MSI_L2A_DS_2BPS_202207...

    array('S2B_OPER_MSI_L2A_DS_2BPS_20220712T125540_S20220712T101515_N04.00',
          dtype='<U64')

  • earthsearch:payload_id
    ()
    <U74
    'roda-sentinel2/workflow-sentine...

    array('roda-sentinel2/workflow-sentinel2-to-stac/167cebdcc5ebde4c6e505d54c546cc02',
          dtype='<U74')

  • s2:thin_cirrus_percentage
    ()
    float64
    0.09388

    array(0.093882)

  • updated
    ()
    <U24
    '2022-11-06T09:30:00.290Z'

    array('2022-11-06T09:30:00.290Z', dtype='<U24')

  • s2:not_vegetated_percentage
    ()
    float64
    13.17

    array(13.165817)

  • s2:granule_id
    ()
    <U62
    'S2B_OPER_MSI_L2A_TL_2BPS_202207...

    array('S2B_OPER_MSI_L2A_TL_2BPS_20220712T125540_A027933_T32TPS_N04.00',
          dtype='<U62')

  • s2:processing_baseline
    ()
    <U5
    '04.00'

    array('04.00', dtype='<U5')

  • s2:cloud_shadow_percentage
    ()
    float64
    3.223

    array(3.222707)

  • proj:epsg
    ()
    int64
    32632

    array(32632)

  • s2:water_percentage
    ()
    float64
    0.3914

    array(0.391429)

  • processing:software
    ()
    object
    {'sentinel2-to-stac': '0.1.0'}

    array({'sentinel2-to-stac': '0.1.0'}, dtype=object)

  • s2:high_proba_clouds_percentage
    ()
    float64
    14.17

    array(14.174789)

  • created
    ()
    <U24
    '2022-11-06T09:30:00.290Z'

    array('2022-11-06T09:30:00.290Z', dtype='<U24')

  • s2:snow_ice_percentage
    ()
    float64
    1.195

    array(1.194618)

  • mgrs:latitude_band
    ()
    <U1
    'T'

    array('T', dtype='<U1')

  • mgrs:grid_square
    ()
    <U2
    'PS'

    array('PS', dtype='<U2')

  • s2:saturated_defective_pixel_percentage
    ()
    int64
    0

    array(0)

  • earthsearch:s3_path
    ()
    <U79
    's3://sentinel-cogs/sentinel-s2-...

    array('s3://sentinel-cogs/sentinel-s2-l2a-cogs/32/T/PS/2022/7/S2B_32TPS_20220712_0_L2A',
          dtype='<U79')

  • constellation
    ()
    <U10
    'sentinel-2'

    array('sentinel-2', dtype='<U10')

  • s2:nodata_pixel_percentage
    ()
    float64
    1e-05

    array(1.e-05)

  • grid:code
    ()
    <U10
    'MGRS-32TPS'

    array('MGRS-32TPS', dtype='<U10')

  • s2:vegetation_percentage
    ()
    float64
    58.37

    array(58.367532)

  • mgrs:utm_zone
    ()
    int64
    32

    array(32)

  • s2:datatake_type
    ()
    <U8
    'INS-NOBS'

    array('INS-NOBS', dtype='<U8')

  • s2:dark_features_percentage
    ()
    float64
    1.684

    array(1.683803)

  • view:sun_elevation
    ()
    float64
    62.31

    array(62.31315147)

  • s2:medium_proba_clouds_percentage
    ()
    float64
    7.346

    array(7.346307)

  • title
    (band)
    <U20
    'Red (band 4) - 10m' ... 'Blue (...

    array(['Red (band 4) - 10m', 'Green (band 3) - 10m',
           'Blue (band 2) - 10m'], dtype='<U20')

  • proj:transform
    ()
    object
    {0, 600000, 10, 5200020, -10}

    array({0, 600000, 10, 5200020, -10}, dtype=object)

  • raster:bands
    ()
    object
    {'nodata': 0, 'data_type': 'uint...

    array({'nodata': 0, 'data_type': 'uint16', 'bits_per_sample': 15, 'spatial_resolution': 10, 'scale': 0.0001, 'offset': -0.1},
          dtype=object)

  • proj:shape
    ()
    object
    {10980}

    array({10980}, dtype=object)

  • gsd
    ()
    int64
    10

    array(10)

  • common_name
    (band)
    <U5
    'red' 'green' 'blue'

    array(['red', 'green', 'blue'], dtype='<U5')

  • center_wavelength
    (band)
    float64
    0.665 0.56 0.49

    array([0.665, 0.56 , 0.49 ])

  • full_width_half_max
    (band)
    float64
    0.038 0.045 0.098

    array([0.038, 0.045, 0.098])

  • epsg
    ()
    int64
    32632

    array(32632)

• Indexes: (4)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2022-07-12 10:18:03.747000'], dtype='datetime64[ns]', name='time', freq=None))

  • band
    PandasIndex

    PandasIndex(Index(['red', 'green', 'blue'], dtype='object', name='band'))

  • x
    PandasIndex

    PandasIndex(Index([673320.0, 673330.0, 673340.0, 673350.0, 673360.0, 673370.0, 673380.0,
           673390.0, 673400.0, 673410.0,
           ...
           684670.0, 684680.0, 684690.0, 684700.0, 684710.0, 684720.0, 684730.0,
           684740.0, 684750.0, 684760.0],
          dtype='float64', name='x', length=1145))

  • y
    PandasIndex

    PandasIndex(Index([5154880.0, 5154870.0, 5154860.0, 5154850.0, 5154840.0, 5154830.0,
           5154820.0, 5154810.0, 5154800.0, 5154790.0,
           ...
           5147850.0, 5147840.0, 5147830.0, 5147820.0, 5147810.0, 5147800.0,
           5147790.0, 5147780.0, 5147770.0, 5147760.0],
          dtype='float64', name='y', length=713))

• Attributes: (4)

  spec :

  RasterSpec(epsg=32632, bounds=(600000.0, 5090220.0, 709800.0, 5200020.0), resolutions_xy=(10.0, 10.0))

  crs :

  epsg:32632

  transform :

  | 10.00, 0.00, 600000.00| | 0.00,-10.00, 5200020.00| | 0.00, 0.00, 1.00|

  resolution :

  10.0

Visualize the RGB bands of our sample dataset:

data = datacube.execute()
data[0].plot.imshow()

/home/749c4c89-469c-4d4a-9ce4-9feffea0504a/.local/lib/python3.11/site-packages/stackstac/prepare.py:408: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.
  times = pd.to_datetime(
Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

<matplotlib.image.AxesImage at 0x7fee562919d0>

Apply an arithmetic formula

We would like to improve the previous visualization, rescaling all the pixels between 0 and 1.

We can use apply in combination with other math processes.

from openeo.processes import linear_scale_range
input_min = -0.1
input_max = 0.2
output_min = 0
output_max = 1

def rescale(x):
    return linear_scale_range(x,input_min,input_max,output_min,output_max)

scaled_data = datacube.apply(rescale)
scaled_data

Visualize the result and see how apply scaled the data resulting in a more meaningful visualization:

scaled_data_xr = scaled_data.execute()
scaled_data_xr[0].plot.imshow()

/home/749c4c89-469c-4d4a-9ce4-9feffea0504a/.local/lib/python3.11/site-packages/stackstac/prepare.py:408: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.
  times = pd.to_datetime(

<matplotlib.image.AxesImage at 0x7fee57f68210>