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()

2024-04-25 13:54:46,683 - [DEBUG] - *** Loading ['NBR'] for 20200824T110631_S2_T29TQE_L1C_150432 ***

2024-04-25 13:54:50,725 - [DEBUG] - Merging bands NBR

2024-04-25 13:54:50,835 - [DEBUG] - Collocating bands

2024-04-25 13:54:50,839 - [DEBUG] - *** Loading ['NBR'] for 20200908T110619_S2_T29TQE_L1C_132324 ***

0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 13:54:53,442 - [DEBUG] - Merging bands NBR

2024-04-25 13:54:53,544 - [DEBUG] - Collocating bands

0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 13:55:10,929 - [DEBUG] - *** Loading ['NBR'] for 20200824T105515_S3_SLSTR_RBT ***

2024-04-25 13:55:33,274 - [DEBUG] - Merging bands NBR

2024-04-25 13:55:33,376 - [DEBUG] - Collocating bands

2024-04-25 13:55:33,380 - [DEBUG] - *** Loading ['NBR'] for 20200909T104016_S3_SLSTR_RBT ***

0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 13:56:00,657 - [DEBUG] - Merging bands NBR

2024-04-25 13:56:00,762 - [DEBUG] - Collocating bands

0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 13:56:08,953 - [DEBUG] - *** Loading ['NBR'] for 20200828T110153_L8_202032_OLI_TIRS ***

2024-04-25 13:56:10,953 - [DEBUG] - Merging bands NBR

2024-04-25 13:56:11,144 - [DEBUG] - Collocating bands

2024-04-25 13:56:11,150 - [DEBUG] - *** Loading ['NBR'] for 20200929T110204_L8_202032_OLI_TIRS ***

0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 13:56:12,783 - [DEBUG] - Merging bands NBR

2024-04-25 13:56:12,948 - [DEBUG] - Collocating bands

0...10...20...30...40...50...60...70...80...90...100 - done.

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

2024-04-25 13:56:16,780 - [DEBUG] - *** Loading ['NBR'] for 20200824T110631_S2_T29TQE_L1C_150432 ***

2024-04-25 13:56:20,859 - [DEBUG] - Merging bands NBR

2024-04-25 13:56:21,018 - [DEBUG] - Collocating bands

2024-04-25 13:56:21,026 - [DEBUG] - *** Loading ['NBR'] for 20200908T110619_S2_T29TQE_L1C_132324 ***

0...10...20...30...40...50...60...70...80...90...100 - done.

2024-04-25 13:56:24,587 - [DEBUG] - Merging bands NBR

2024-04-25 13:56:24,700 - [DEBUG] - Collocating bands

0...10...20...30...40...50...60...70...80...90...100 - done.

<xarray.DataArray 'Reference_NBR_Secondary_NBR_dNBR' (bands: 3, y: 271, x: 272)> Size: 885kB
dask.array<transpose, shape=(3, 271, 272), dtype=float32, chunksize=(1, 271, 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) object 24B '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: 271
• x: 272

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

  Array Chunk Bytes 863.81 kiB 287.94 kiB Shape (3, 271, 272) (1, 271, 272) Dask graph 3 chunks in 11 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])

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

    array([4458110.01845 , 4458050.055351, 4457990.092251, ..., 4442039.907749,
           4441979.944649, 4441919.98155 ])

  • 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.073529411764554 0.0 4458140.0 0.0 -59.96309963099797

    array(0)

  • bands
    (bands)
    object
    MultiIndex

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

  • variable
    (bands)
    object
    'Reference_NBR' ... 'dNBR'

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

  • 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.3014705882,        786289.375, 786349.4485294117,
           786409.5220588235, 786469.5955882353,  786529.669117647,
           786589.7426470588, 786649.8161764706, 786709.8897058823,
           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,
           ...
           4442459.649446494, 4442399.686346863, 4442339.723247232,
           4442279.760147601,  4442219.79704797, 4442159.833948339,
           4442099.870848708, 4442039.907749077, 4441979.944649446,
           4441919.981549815],
          dtype='float64', name='y', length=271))

  • 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 0x7f4ebc13b160>