Pair

Basic examples of pairs, trying to extract fires. Here the reference product will be the pre-event image, the secondary being the post-event.

Initialization of the notebook

• Imports

• Loggers

• Paths

from cloudpathlib import AnyPath
from eosets.pair import Pair
from eoreader.bands import NBR
from matplotlib import pyplot as plt
import cartopy.crs as ccrs

# Create logger
import logging
from sertit import logs

logs.init_logger(logging.getLogger("eoreader"), log_lvl=logging.INFO)
logs.init_logger(logging.getLogger("eosets"))

# Get the base paths
data_path = AnyPath(r"/home/ds2_db3/CI/eosets/DATA")

# Get a dictionary of Sentinel-2, Sentinel-3 and Landsat-8 reference and secondary images
s2_paths = {
    "reference_paths": [data_path / "S2A_MSIL1C_20200824T110631_N0209_R137_T29TQE_20200824T150432.zip"],
    "secondary_paths": [data_path / "S2B_MSIL1C_20200908T110619_N0209_R137_T29TQE_20200908T132324.zip"],
}

s3_paths = {
    "reference_paths": [data_path / "S3B_SL_1_RBT____20200824T105515_20200824T105815_20200825T151744_0179_042_322_2340_LN2_O_NT_004.SEN3"],
    "secondary_paths": [data_path / "S3B_SL_1_RBT____20200909T104016_20200909T104316_20200910T161910_0179_043_165_2340_LN2_O_NT_004.SEN3"],
}

l8_paths = {
    "reference_paths": [data_path / "LC08_L1TP_202032_20200828_20200906_02_T1.tar"],
    "secondary_paths": [data_path / "LC08_L1TP_202032_20200929_20201006_02_T1.tar"],
}

all_paths = [s2_paths, s3_paths, l8_paths]

aoi_path = data_path / "Fire_Spain.geojson"

Creation of the Pairs and loading of dNBR

To create a pair, you just have to give the paths of your products. The reference and secondary products should overlap.

The goal here will be to extract the dNBR band of the pair and display it for every pair.

\[ dNBR = NBR_{pre} - NBR_{post}\]

\[ NBR = \frac{NIR - SWIR}{NIR + SWIR} \]

for paths in all_paths:
    # Create pair with reference and secondary paths
    pair = Pair(**paths)

    # Load dNBR
    _, _, dnbr_ds = pair.load(
        diff_bands=NBR,
        window=aoi_path,
        pixel_size=60
    )
    dnbr = dnbr_ds[NBR]
    dnbr.attrs["long_name"] = f"dNBR of {pair.condensed_name}"

    # Plot the pairs
    crs = ccrs.UTM(str(dnbr.rio.crs.to_epsg())[-2:])
    dnbr.plot(
        robust=True,
        transform=crs,
        x="x",
        y="y",
        cmap="gist_heat",
        subplot_kws={'projection':crs}
    )
    plt.show()

2026-01-06 16:43:03,722 - [DEBUG] - *** Loading ['NBR'] for 20200824T110631_S2_T29TQE_L1C_150432 ***

/opt/conda/lib/python3.11/site-packages/rasterio/warp.py:387: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
  dest = _reproject(

2026-01-06 16:43:08,628 - [DEBUG] - Collocating bands

2026-01-06 16:43:08,632 - [DEBUG] - *** Loading ['NBR'] for 20200908T110619_S2_T29TQE_L1C_132324 ***

2026-01-06 16:43:12,531 - [DEBUG] - Collocating bands

2026-01-06 16:43:14,889 - [DEBUG] - *** Loading ['NBR'] for 20200824T105515_S3_SLSTR_RBT ***

2026-01-06 16:43:59,281 - [DEBUG] - Collocating bands

2026-01-06 16:43:59,284 - [DEBUG] - *** Loading ['NBR'] for 20200909T104016_S3_SLSTR_RBT ***

2026-01-06 16:45:17,318 - [DEBUG] - Collocating bands

'_vectorize' function is not lazy yet. Computing the raster.

'_vectorize' function is not lazy yet. Computing the raster.

2026-01-06 16:45:26,786 - [DEBUG] - *** Loading ['NBR'] for 20200828T110153_L8_202032_OLI_TIRS ***

2026-01-06 16:45:29,287 - [DEBUG] - Collocating bands

2026-01-06 16:45:29,289 - [DEBUG] - *** Loading ['NBR'] for 20200929T110204_L8_202032_OLI_TIRS ***

2026-01-06 16:45:31,708 - [DEBUG] - Collocating bands

Stacking pairs

A pair can also be stacked.

pair = Pair(**s2_paths)
stack = pair.stack(
    reference_bands=NBR,
    secondary_bands=NBR,
    diff_bands=NBR,
    window=aoi_path,
    pixel_size=60
    )
stack

2026-01-06 16:45:44,344 - [DEBUG] - *** Loading ['NBR'] for 20200824T110631_S2_T29TQE_L1C_150432 ***

2026-01-06 16:45:47,470 - [DEBUG] - Collocating bands

2026-01-06 16:45:47,473 - [DEBUG] - *** Loading ['NBR'] for 20200908T110619_S2_T29TQE_L1C_132324 ***

2026-01-06 16:45:50,148 - [DEBUG] - Collocating bands

<xarray.DataArray 'Reference_NBR_Secondary_NBR_dNBR' (bands: 3, y: 272, x: 272)> Size: 888kB
dask.array<where, shape=(3, 272, 272), dtype=float32, chunksize=(1, 272, 272), chunktype=numpy.ndarray>
Coordinates:
  * x            (x) float64 2kB 7.705e+05 7.706e+05 ... 7.867e+05 7.868e+05
  * y            (y) float64 2kB 4.458e+06 4.458e+06 ... 4.442e+06 4.442e+06
    spatial_ref  int64 8B 0
  * bands        (bands) object 24B MultiIndex
  * variable     (bands) <U13 156B 'Reference_NBR' 'Secondary_NBR' 'dNBR'
  * band         (bands) int64 24B 1 1 1
Attributes:
    long_name:       Reference_NBR Secondary_NBR dNBR
    condensed_name:  20200824_S2_20200908_S2

xarray.DataArray

'Reference_NBR_Secondary_NBR_dNBR'

• bands: 3
• y: 272
• x: 272

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

  Array Chunk Bytes 867.00 kiB 289.00 kiB Shape (3, 272, 272) (1, 272, 272) Dask graph 3 chunks in 16 graph layers Data type float32 numpy.ndarray

• Coordinates: (6)
  • x
    (x)
    float64
    7.705e+05 7.706e+05 ... 7.868e+05

    array([770490.036765, 770550.110294, 770610.183824, ..., 786649.816176,
           786709.889706, 786769.963235], shape=(272,))

  • y
    (y)
    float64
    4.458e+06 4.458e+06 ... 4.442e+06

    array([4458110.01845 , 4458050.055351, 4457990.092251, ..., 4441979.944649,
           4441919.98155 , 4441860.01845 ], shape=(272,))

  • spatial_ref
    ()
    int64
    0

    spatial_ref :

    PROJCS["WGS 84 / UTM zone 29N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32629"]]

    crs_wkt :

    PROJCS["WGS 84 / UTM zone 29N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32629"]]

    GeoTransform :

    770460.0 60.07352941176471 0.0 4458140.0 0.0 -59.96309963099631

    array(0)

  • bands
    (bands)
    object
    MultiIndex

    array([('Reference_NBR', 1), ('Secondary_NBR', 1), ('dNBR', 1)], dtype=object)

  • variable
    (bands)
    <U13
    'Reference_NBR' ... 'dNBR'

    array(['Reference_NBR', 'Secondary_NBR', 'dNBR'], dtype='<U13')

  • band
    (bands)
    int64
    1 1 1

    array([1, 1, 1])

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Index([770490.0367647059, 770550.1102941177, 770610.1838235294,
           770670.2573529412,  770730.330882353, 770790.4044117647,
           770850.4779411765, 770910.5514705883,        770970.625,
           771030.6985294118,
           ...
           786229.3014705883,        786289.375, 786349.4485294118,
           786409.5220588235, 786469.5955882353, 786529.6691176471,
           786589.7426470588, 786649.8161764706, 786709.8897058824,
           786769.9632352941],
          dtype='float64', name='x', length=272))

  • y
    PandasIndex

    PandasIndex(Index([4458110.018450185, 4458050.055350554, 4457990.092250923,
           4457930.129151292, 4457870.166051661,  4457810.20295203,
           4457750.239852399, 4457690.276752768, 4457630.313653137,
           4457570.350553506,
           ...
           4442399.686346863, 4442339.723247233, 4442279.760147601,
           4442219.797047971, 4442159.833948339, 4442099.870848709,
           4442039.907749077, 4441979.944649447, 4441919.981549815,
           4441860.018450185],
          dtype='float64', name='y', length=272))

  • bands
    variable
    band
    PandasMultiIndex

    PandasIndex(MultiIndex([('Reference_NBR', 1),
                ('Secondary_NBR', 1),
                (         'dNBR', 1)],
               name='bands'))

• Attributes: (2)

  long_name :

  Reference_NBR Secondary_NBR dNBR

  condensed_name :

  20200824_S2_20200908_S2

stack.plot.imshow(vmin=0, vmax=1, robust=True)

<matplotlib.image.AxesImage at 0x7f7b04376190>