Interactive dashboard for Sentinel-2 satellite imagery

---

Overview

In this notebook, we will take a look at how to retrieve Sentinel-2 L2A satellite imagery from the Microsoft Planetary Computer Data Catalog (MSPC). We will go over how to interact with the Data Catalog, which exposes a SpatioTemporal Asset Catalog (STAC) interface for querying, searching and retrieving data. We will use the odc-stac package to load the data lazily, which means data is not actually read unless required (say, for plotting). Once loaded, we will process the data and make a simple interactive dashboard to look at the satellite imagery over a location for different seasons. We will use the HoloViz ecosystem for the interactive dashboard.

Prerequisites

Concepts	Importance	Notes
Xarray	Helpful	Background
Dask + Xarray	Necessary	Background
About the Microsoft Planetary Computer (MSPC)	Helpful	Background
Documentation of pystac-client	Helpful	Consult as needed
Landsat ML Cookbook	Helpful	Similar cookbook, illustrates accessing Landsat data from MSPC
About the HoloViz ecosystem	Helpful	How different HoloViz packages work with each other
Sentinel-2 L2A User Guide	Necessary	Background about the satellite data
Sentinel-2 L2A data definitions	Helpful	Section 2.3.10 has some useful information about the data we access in this cookbook

• Time to learn: 15 minutes

Imports

import os
import pandas as pd
import numpy as np
import xarray as xr
import pystac_client
import planetary_computer
import panel as pn
import panel.widgets as pnw
import hvplot.xarray
import holoviews as hv
import geoviews as gv
from pystac.extensions.eo import EOExtension as eo
import datetime
from cartopy import crs
from dask.distributed import Client, LocalCluster
import odc.stac

xr.set_options(keep_attrs=True)
hv.extension('bokeh')
gv.extension('bokeh')

Since we will use dask to distribute our computation, we can create a dask cluster and connect to it via a dask client. For this recipe, we will create a LocalCluster with default values.

cluster = LocalCluster()
client = Client(cluster)
client

Client

Client-f8e44681-122e-11ef-8976-6045bd854962

Connection method: Cluster object	Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

LocalCluster

2bddd65c

Dashboard: http://127.0.0.1:8787/status	Workers: 4
Total threads: 4	Total memory: 15.61 GiB
Status: running	Using processes: True

Scheduler Info

Scheduler

Scheduler-d05ea11b-edb2-4860-aa4c-af25c0e30f8a

Comm: tcp://127.0.0.1:42743	Workers: 4
Dashboard: http://127.0.0.1:8787/status	Total threads: 4
Started: Just now	Total memory: 15.61 GiB

Workers

Worker: 0

Comm: tcp://127.0.0.1:34403	Total threads: 1
Dashboard: http://127.0.0.1:46707/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:46415
Local directory: /tmp/dask-scratch-space/worker-ifecgyrk

Worker: 1

Comm: tcp://127.0.0.1:38251	Total threads: 1
Dashboard: http://127.0.0.1:34881/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:39967
Local directory: /tmp/dask-scratch-space/worker-ye9fst9k

Worker: 2

Comm: tcp://127.0.0.1:36149	Total threads: 1
Dashboard: http://127.0.0.1:40137/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:34365
Local directory: /tmp/dask-scratch-space/worker-0l5lwuey

Worker: 3

Comm: tcp://127.0.0.1:46269	Total threads: 1
Dashboard: http://127.0.0.1:34987/status	Memory: 3.90 GiB
Nanny: tcp://127.0.0.1:42301
Local directory: /tmp/dask-scratch-space/worker-d98xynn5

Open the catalog

The root of Microsoft Planetary Computer’s STAC API endpoint is located at https://planetarycomputer.microsoft.com/api/stac/v1. We will load in the catalog using the pystac_client.Client.open method. Even though the STAC metadata in MSPC is publicly accessible, authentication is required to access the actual data. The modifier keyword can be used to explicitly “sign” an item, which basically means we can access the privately stored data (more information here).

stac_root = 'https://planetarycomputer.microsoft.com/api/stac/v1'
catalog = pystac_client.Client.open(
    stac_root,
    modifier=planetary_computer.sign_inplace
)
print(f"{catalog.title} - {catalog.description}")

Microsoft Planetary Computer STAC API - Searchable spatiotemporal metadata describing Earth science datasets hosted by the Microsoft Planetary Computer

Let’s search for any collections that have the substring “sentinel-2” in it to discover the sentinel-2 data.

sentinel2_collections = [collection for collection in catalog.get_collections() if "sentinel-2" in collection.id]
sentinel2_collections

[<CollectionClient id=sentinel-2-l2a>]

Looks like there is only one collection (sentinel-2-l2a) available in the catalog - which is the dataset we want to use.

Query, Filter and Load the collection

Now that we have the ID of the collection of interest, we can specify certain filters to narrow down to exactly the data we want to look at. The final visualization in the recipe will look at how the NCAR Mesa Lab, Boulder CO looks like throughout the year as seen from space. To narrow down our search, we will use the following criteria -

• Bounding box: We will limit our spatial extent to the bounding box of the NCAR Mesa Lab region.

• Date range: We will look at the year 2022

• Collection: sentinel-2-l2a (from previous cell)

• Cloud threshold: Since cloud can block the satellite from making an observation of the ground, we will limit our search to satellite images where the cloud cover is less than a certain threshold, 40% in this case.

Feel free to change these filtering paramters to suit your needs when running in an interactive session.

bbox = [-105.283263,39.972809,-105.266569,39.987640] # NCAR, boulder, CO. bbox from http://bboxfinder.com/
date_range = "2022-01-01/2022-12-31"
collection = "sentinel-2-l2a"                        # full id of collection
cloud_thresh = 30

search = catalog.search(
    collections = sentinel2_collections,
    bbox = bbox,
    datetime = date_range,
    query={"eo:cloud_cover": {"lt": cloud_thresh}}
)
items = search.item_collection()
print(f"Found {len(items)} items in the {collection}")

Found 64 items in the sentinel-2-l2a

We now have an ItemCollection with the data that we need. Let’s look at one of the items in the collection and explore what assets it has.

first_item = items.items[0]
all_bands = list(first_item.assets.keys())
print("Assets available:")
print(*all_bands, sep=', ')

Assets available:
AOT, B01, B02, B03, B04, B05, B06, B07, B08, B09, B11, B12, B8A, SCL, WVP, visual, preview, safe-manifest, granule-metadata, inspire-metadata, product-metadata, datastrip-metadata, tilejson, rendered_preview

Seems like there are a lot of assets associated with the item - you can read about them here. We are interested in the assets that start with a ‘B’, which are the bands associated with the different spectral bands in which the MultiSpectral Instrument (MSI) of the satellite captures observations. Specifically, the RGB – or Red, Green and Blue - bands that we need to create a “True color” image are as follows:

band	corresponds to
B04	Red
B03	Green
B02	Blue

We will use the odc.stac.stac_load function to load in the assets that start with the alphabet ‘B’. This function will return a lazily-loaded xr.DataSet (using dask). For plotting purposes it is better if we have the data as a xr.DataArray instead with the bands as a dimension. We can do that using .to_array(dim=<dim_name>) method of a dataset.

bands_of_interest = [b for b in all_bands if b.startswith('B')]

da = odc.stac.stac_load(
    items,
    bands=bands_of_interest,
    bbox=bbox,
    chunks={},  # <-- use Dask
).to_array(dim='band')
da

<xarray.DataArray (band: 12, time: 64, y: 166, x: 144)> Size: 73MB
dask.array<stack, shape=(12, 64, 166, 144), dtype=float32, chunksize=(1, 1, 166, 144), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 1kB 4.426e+06 4.426e+06 ... 4.425e+06 4.425e+06
  * x            (x) float64 1kB 4.758e+05 4.758e+05 ... 4.772e+05 4.772e+05
    spatial_ref  int32 4B 32613
  * time         (time) datetime64[ns] 512B 2022-01-26T17:56:21.024000 ... 20...
  * band         (band) object 96B 'B01' 'B02' 'B03' 'B04' ... 'B11' 'B12' 'B8A'

xarray.DataArray

• band: 12
• time: 64
• y: 166
• x: 144

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

  Array Chunk Bytes 70.03 MiB 93.38 kiB Shape (12, 64, 166, 144) (1, 1, 166, 144) Dask graph 768 chunks in 13 graph layers Data type float32 numpy.ndarray

• Coordinates: (5)
  • y
    (y)
    float64
    4.426e+06 4.426e+06 ... 4.425e+06

    units :

    metre

    resolution :

    -10.0

    crs :

    EPSG:32613

    array([4426425., 4426415., 4426405., 4426395., 4426385., 4426375., 4426365.,
           4426355., 4426345., 4426335., 4426325., 4426315., 4426305., 4426295.,
           4426285., 4426275., 4426265., 4426255., 4426245., 4426235., 4426225.,
           4426215., 4426205., 4426195., 4426185., 4426175., 4426165., 4426155.,
           4426145., 4426135., 4426125., 4426115., 4426105., 4426095., 4426085.,
           4426075., 4426065., 4426055., 4426045., 4426035., 4426025., 4426015.,
           4426005., 4425995., 4425985., 4425975., 4425965., 4425955., 4425945.,
           4425935., 4425925., 4425915., 4425905., 4425895., 4425885., 4425875.,
           4425865., 4425855., 4425845., 4425835., 4425825., 4425815., 4425805.,
           4425795., 4425785., 4425775., 4425765., 4425755., 4425745., 4425735.,
           4425725., 4425715., 4425705., 4425695., 4425685., 4425675., 4425665.,
           4425655., 4425645., 4425635., 4425625., 4425615., 4425605., 4425595.,
           4425585., 4425575., 4425565., 4425555., 4425545., 4425535., 4425525.,
           4425515., 4425505., 4425495., 4425485., 4425475., 4425465., 4425455.,
           4425445., 4425435., 4425425., 4425415., 4425405., 4425395., 4425385.,
           4425375., 4425365., 4425355., 4425345., 4425335., 4425325., 4425315.,
           4425305., 4425295., 4425285., 4425275., 4425265., 4425255., 4425245.,
           4425235., 4425225., 4425215., 4425205., 4425195., 4425185., 4425175.,
           4425165., 4425155., 4425145., 4425135., 4425125., 4425115., 4425105.,
           4425095., 4425085., 4425075., 4425065., 4425055., 4425045., 4425035.,
           4425025., 4425015., 4425005., 4424995., 4424985., 4424975., 4424965.,
           4424955., 4424945., 4424935., 4424925., 4424915., 4424905., 4424895.,
           4424885., 4424875., 4424865., 4424855., 4424845., 4424835., 4424825.,
           4424815., 4424805., 4424795., 4424785., 4424775.])

  • x
    (x)
    float64
    4.758e+05 4.758e+05 ... 4.772e+05

    units :

    metre

    resolution :

    10.0

    crs :

    EPSG:32613

    array([475815., 475825., 475835., 475845., 475855., 475865., 475875., 475885.,
           475895., 475905., 475915., 475925., 475935., 475945., 475955., 475965.,
           475975., 475985., 475995., 476005., 476015., 476025., 476035., 476045.,
           476055., 476065., 476075., 476085., 476095., 476105., 476115., 476125.,
           476135., 476145., 476155., 476165., 476175., 476185., 476195., 476205.,
           476215., 476225., 476235., 476245., 476255., 476265., 476275., 476285.,
           476295., 476305., 476315., 476325., 476335., 476345., 476355., 476365.,
           476375., 476385., 476395., 476405., 476415., 476425., 476435., 476445.,
           476455., 476465., 476475., 476485., 476495., 476505., 476515., 476525.,
           476535., 476545., 476555., 476565., 476575., 476585., 476595., 476605.,
           476615., 476625., 476635., 476645., 476655., 476665., 476675., 476685.,
           476695., 476705., 476715., 476725., 476735., 476745., 476755., 476765.,
           476775., 476785., 476795., 476805., 476815., 476825., 476835., 476845.,
           476855., 476865., 476875., 476885., 476895., 476905., 476915., 476925.,
           476935., 476945., 476955., 476965., 476975., 476985., 476995., 477005.,
           477015., 477025., 477035., 477045., 477055., 477065., 477075., 477085.,
           477095., 477105., 477115., 477125., 477135., 477145., 477155., 477165.,
           477175., 477185., 477195., 477205., 477215., 477225., 477235., 477245.])

  • spatial_ref
    ()
    int32
    32613

    spatial_ref :

    PROJCRS["WGS 84 / UTM zone 13N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 13N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-105,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 108°W and 102°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],BBOX[0,-108,84,-102]],ID["EPSG",32613]]

    crs_wkt :

    PROJCRS["WGS 84 / UTM zone 13N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 13N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-105,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 108°W and 102°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],BBOX[0,-108,84,-102]],ID["EPSG",32613]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    projected_crs_name :

    WGS 84 / UTM zone 13N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -105.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    GeoTransform :

    475810 10 0 4426430 0 -10

    array(32613, dtype=int32)

  • time
    (time)
    datetime64[ns]
    2022-01-26T17:56:21.024000 ... 2...

    array(['2022-01-26T17:56:21.024000000', '2022-01-31T17:55:49.024000000',
           '2022-02-05T17:55:31.024000000', '2022-02-07T17:45:09.024000000',
           '2022-02-15T17:54:31.024000000', '2022-02-20T17:53:49.025000000',
           '2022-02-27T17:42:49.024000000', '2022-03-07T17:52:01.024000000',
           '2022-03-14T17:41:11.024000000', '2022-03-19T17:39:39.024000000',
           '2022-04-01T17:49:09.024000000', '2022-04-08T17:38:59.024000000',
           '2022-04-18T17:38:59.024000000', '2022-04-21T17:48:59.024000000',
           '2022-04-26T17:49:11.024000000', '2022-05-11T17:49:09.024000000',
           '2022-05-18T17:38:59.024000000', '2022-05-26T17:49:21.024000000',
           '2022-06-02T17:39:11.024000000', '2022-06-07T17:39:09.024000000',
           '2022-06-15T17:49:21.024000000', '2022-06-17T17:39:09.024000000',
           '2022-06-20T17:49:09.024000000', '2022-06-27T17:39:09.024000000',
           '2022-07-02T17:39:21.024000000', '2022-07-07T17:39:09.024000000',
           '2022-07-12T17:39:21.024000000', '2022-07-17T17:39:09.024000000',
           '2022-07-22T17:39:11.024000000', '2022-07-25T17:49:21.024000000',
           '2022-07-27T17:39:09.024000000', '2022-07-30T17:49:09.024000000',
           '2022-08-04T17:49:21.024000000', '2022-08-06T17:39:09.024000000',
           '2022-08-09T17:49:09.024000000', '2022-08-11T17:39:21.024000000',
           '2022-08-21T17:39:21.024000000', '2022-08-24T17:49:21.024000000',
           '2022-08-29T17:49:09.024000000', '2022-08-31T17:39:21.024000000',
           '2022-09-03T17:49:21.024000000', '2022-09-05T17:39:09.024000000',
           '2022-09-08T17:49:09.024000000', '2022-09-18T17:50:09.024000000',
           '2022-09-23T17:50:51.024000000', '2022-09-25T17:40:59.024000000',
           '2022-10-05T17:41:59.024000000', '2022-10-08T17:52:29.024000000',
           '2022-10-10T17:42:51.024000000', '2022-10-13T17:53:11.024000000',
           '2022-10-18T17:53:29.024000000', '2022-10-25T17:44:19.024000000',
           '2022-10-28T17:54:39.024000000', '2022-10-30T17:44:51.024000000',
           '2022-11-02T17:55:11.024000000', '2022-11-07T17:55:39.024000000',
           '2022-11-09T17:45:51.024000000', '2022-11-19T17:46:31.024000000',
           '2022-11-22T17:56:41.024000000', '2022-12-04T17:47:19.024000000',
           '2022-12-09T17:47:31.024000000', '2022-12-17T17:57:39.024000000',
           '2022-12-22T17:57:41.024000000', '2022-12-29T17:47:41.024000000'],
          dtype='datetime64[ns]')

  • band
    (band)
    object
    'B01' 'B02' 'B03' ... 'B12' 'B8A'

    array(['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B11',
           'B12', 'B8A'], dtype=object)

• Indexes: (4)
  • y
    PandasIndex

    PandasIndex(Index([4426425.0, 4426415.0, 4426405.0, 4426395.0, 4426385.0, 4426375.0,
           4426365.0, 4426355.0, 4426345.0, 4426335.0,
           ...
           4424865.0, 4424855.0, 4424845.0, 4424835.0, 4424825.0, 4424815.0,
           4424805.0, 4424795.0, 4424785.0, 4424775.0],
          dtype='float64', name='y', length=166))

  • x
    PandasIndex

    PandasIndex(Index([475815.0, 475825.0, 475835.0, 475845.0, 475855.0, 475865.0, 475875.0,
           475885.0, 475895.0, 475905.0,
           ...
           477155.0, 477165.0, 477175.0, 477185.0, 477195.0, 477205.0, 477215.0,
           477225.0, 477235.0, 477245.0],
          dtype='float64', name='x', length=144))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2022-01-26 17:56:21.024000', '2022-01-31 17:55:49.024000',
                   '2022-02-05 17:55:31.024000', '2022-02-07 17:45:09.024000',
                   '2022-02-15 17:54:31.024000', '2022-02-20 17:53:49.025000',
                   '2022-02-27 17:42:49.024000', '2022-03-07 17:52:01.024000',
                   '2022-03-14 17:41:11.024000', '2022-03-19 17:39:39.024000',
                   '2022-04-01 17:49:09.024000', '2022-04-08 17:38:59.024000',
                   '2022-04-18 17:38:59.024000', '2022-04-21 17:48:59.024000',
                   '2022-04-26 17:49:11.024000', '2022-05-11 17:49:09.024000',
                   '2022-05-18 17:38:59.024000', '2022-05-26 17:49:21.024000',
                   '2022-06-02 17:39:11.024000', '2022-06-07 17:39:09.024000',
                   '2022-06-15 17:49:21.024000', '2022-06-17 17:39:09.024000',
                   '2022-06-20 17:49:09.024000', '2022-06-27 17:39:09.024000',
                   '2022-07-02 17:39:21.024000', '2022-07-07 17:39:09.024000',
                   '2022-07-12 17:39:21.024000', '2022-07-17 17:39:09.024000',
                   '2022-07-22 17:39:11.024000', '2022-07-25 17:49:21.024000',
                   '2022-07-27 17:39:09.024000', '2022-07-30 17:49:09.024000',
                   '2022-08-04 17:49:21.024000', '2022-08-06 17:39:09.024000',
                   '2022-08-09 17:49:09.024000', '2022-08-11 17:39:21.024000',
                   '2022-08-21 17:39:21.024000', '2022-08-24 17:49:21.024000',
                   '2022-08-29 17:49:09.024000', '2022-08-31 17:39:21.024000',
                   '2022-09-03 17:49:21.024000', '2022-09-05 17:39:09.024000',
                   '2022-09-08 17:49:09.024000', '2022-09-18 17:50:09.024000',
                   '2022-09-23 17:50:51.024000', '2022-09-25 17:40:59.024000',
                   '2022-10-05 17:41:59.024000', '2022-10-08 17:52:29.024000',
                   '2022-10-10 17:42:51.024000', '2022-10-13 17:53:11.024000',
                   '2022-10-18 17:53:29.024000', '2022-10-25 17:44:19.024000',
                   '2022-10-28 17:54:39.024000', '2022-10-30 17:44:51.024000',
                   '2022-11-02 17:55:11.024000', '2022-11-07 17:55:39.024000',
                   '2022-11-09 17:45:51.024000', '2022-11-19 17:46:31.024000',
                   '2022-11-22 17:56:41.024000', '2022-12-04 17:47:19.024000',
                   '2022-12-09 17:47:31.024000', '2022-12-17 17:57:39.024000',
                   '2022-12-22 17:57:41.024000', '2022-12-29 17:47:41.024000'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • band
    PandasIndex

    PandasIndex(Index(['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B11',
           'B12', 'B8A'],
          dtype='object', name='band'))

• Attributes: (0)

Prepare the data for visualization

On January 25th, 2022, the European Space Agency (data provider for Sentinel-2 satellite) made a change in their processing pipeline to address some issues that you can read about here if interested. For the purpose of this notebook, we will process the newer dataset such that it is harmonized with the old processing pipeline - in simple words, we will make sure that the data has the same statistical properties so that we can visualize them seamlessly.

# from https://planetarycomputer.microsoft.com/dataset/sentinel-2-l2a#Baseline-Change
def harmonize_to_old(data):  
    """
    Harmonize new Sentinel-2 data to the old baseline.

    Parameters
    ----------
    data: xarray.DataArray
        A DataArray with four dimensions: time, band, y, x

    Returns
    -------
    harmonized: xarray.DataArray
        A DataArray with all values harmonized to the old
        processing baseline.
    """
    cutoff = datetime.datetime(2022, 1, 25)
    offset = 1000
    bands = ["B01","B02","B03","B04","B05","B06","B07","B08","B8A","B09","B10","B11","B12"]

    old = data.sel(time=slice(cutoff))

    to_process = list(set(bands) & set(data.band.data.tolist()))
    new = data.sel(time=slice(cutoff, None)).drop_sel(band=to_process)

    new_harmonized = data.sel(time=slice(cutoff, None), band=to_process).clip(offset)
    new_harmonized -= offset

    new = xr.concat([new, new_harmonized], "band").sel(band=data.band.data.tolist())
    return xr.concat([old, new], dim="time")

da = harmonize_to_old(da)
da

<xarray.DataArray (band: 12, time: 64, y: 166, x: 144)> Size: 73MB
dask.array<getitem, shape=(12, 64, 166, 144), dtype=float32, chunksize=(1, 1, 166, 144), chunktype=numpy.ndarray>
Coordinates:
  * y            (y) float64 1kB 4.426e+06 4.426e+06 ... 4.425e+06 4.425e+06
  * x            (x) float64 1kB 4.758e+05 4.758e+05 ... 4.772e+05 4.772e+05
    spatial_ref  int32 4B 32613
  * time         (time) datetime64[ns] 512B 2022-01-26T17:56:21.024000 ... 20...
  * band         (band) object 96B 'B01' 'B02' 'B03' 'B04' ... 'B11' 'B12' 'B8A'

xarray.DataArray

• band: 12
• time: 64
• y: 166
• x: 144

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

  Array Chunk Bytes 70.03 MiB 93.38 kiB Shape (12, 64, 166, 144) (1, 1, 166, 144) Dask graph 768 chunks in 17 graph layers Data type float32 numpy.ndarray

• Coordinates: (5)
  • y
    (y)
    float64
    4.426e+06 4.426e+06 ... 4.425e+06

    units :

    metre

    resolution :

    -10.0

    crs :

    EPSG:32613

    array([4426425., 4426415., 4426405., 4426395., 4426385., 4426375., 4426365.,
           4426355., 4426345., 4426335., 4426325., 4426315., 4426305., 4426295.,
           4426285., 4426275., 4426265., 4426255., 4426245., 4426235., 4426225.,
           4426215., 4426205., 4426195., 4426185., 4426175., 4426165., 4426155.,
           4426145., 4426135., 4426125., 4426115., 4426105., 4426095., 4426085.,
           4426075., 4426065., 4426055., 4426045., 4426035., 4426025., 4426015.,
           4426005., 4425995., 4425985., 4425975., 4425965., 4425955., 4425945.,
           4425935., 4425925., 4425915., 4425905., 4425895., 4425885., 4425875.,
           4425865., 4425855., 4425845., 4425835., 4425825., 4425815., 4425805.,
           4425795., 4425785., 4425775., 4425765., 4425755., 4425745., 4425735.,
           4425725., 4425715., 4425705., 4425695., 4425685., 4425675., 4425665.,
           4425655., 4425645., 4425635., 4425625., 4425615., 4425605., 4425595.,
           4425585., 4425575., 4425565., 4425555., 4425545., 4425535., 4425525.,
           4425515., 4425505., 4425495., 4425485., 4425475., 4425465., 4425455.,
           4425445., 4425435., 4425425., 4425415., 4425405., 4425395., 4425385.,
           4425375., 4425365., 4425355., 4425345., 4425335., 4425325., 4425315.,
           4425305., 4425295., 4425285., 4425275., 4425265., 4425255., 4425245.,
           4425235., 4425225., 4425215., 4425205., 4425195., 4425185., 4425175.,
           4425165., 4425155., 4425145., 4425135., 4425125., 4425115., 4425105.,
           4425095., 4425085., 4425075., 4425065., 4425055., 4425045., 4425035.,
           4425025., 4425015., 4425005., 4424995., 4424985., 4424975., 4424965.,
           4424955., 4424945., 4424935., 4424925., 4424915., 4424905., 4424895.,
           4424885., 4424875., 4424865., 4424855., 4424845., 4424835., 4424825.,
           4424815., 4424805., 4424795., 4424785., 4424775.])

  • x
    (x)
    float64
    4.758e+05 4.758e+05 ... 4.772e+05

    units :

    metre

    resolution :

    10.0

    crs :

    EPSG:32613

    array([475815., 475825., 475835., 475845., 475855., 475865., 475875., 475885.,
           475895., 475905., 475915., 475925., 475935., 475945., 475955., 475965.,
           475975., 475985., 475995., 476005., 476015., 476025., 476035., 476045.,
           476055., 476065., 476075., 476085., 476095., 476105., 476115., 476125.,
           476135., 476145., 476155., 476165., 476175., 476185., 476195., 476205.,
           476215., 476225., 476235., 476245., 476255., 476265., 476275., 476285.,
           476295., 476305., 476315., 476325., 476335., 476345., 476355., 476365.,
           476375., 476385., 476395., 476405., 476415., 476425., 476435., 476445.,
           476455., 476465., 476475., 476485., 476495., 476505., 476515., 476525.,
           476535., 476545., 476555., 476565., 476575., 476585., 476595., 476605.,
           476615., 476625., 476635., 476645., 476655., 476665., 476675., 476685.,
           476695., 476705., 476715., 476725., 476735., 476745., 476755., 476765.,
           476775., 476785., 476795., 476805., 476815., 476825., 476835., 476845.,
           476855., 476865., 476875., 476885., 476895., 476905., 476915., 476925.,
           476935., 476945., 476955., 476965., 476975., 476985., 476995., 477005.,
           477015., 477025., 477035., 477045., 477055., 477065., 477075., 477085.,
           477095., 477105., 477115., 477125., 477135., 477145., 477155., 477165.,
           477175., 477185., 477195., 477205., 477215., 477225., 477235., 477245.])

  • spatial_ref
    ()
    int32
    32613

    spatial_ref :

    PROJCRS["WGS 84 / UTM zone 13N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 13N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-105,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 108°W and 102°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],BBOX[0,-108,84,-102]],ID["EPSG",32613]]

    crs_wkt :

    PROJCRS["WGS 84 / UTM zone 13N",BASEGEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],ID["EPSG",4326]],CONVERSION["UTM zone 13N",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-105,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9996,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",500000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1]],USAGE[SCOPE["Navigation and medium accuracy spatial referencing."],AREA["Between 108°W and 102°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - Northwest Territories (NWT); Nunavut; Saskatchewan. Mexico. United States (USA)."],BBOX[0,-108,84,-102]],ID["EPSG",32613]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984 ensemble

    projected_crs_name :

    WGS 84 / UTM zone 13N

    grid_mapping_name :

    transverse_mercator

    latitude_of_projection_origin :

    0.0

    longitude_of_central_meridian :

    -105.0

    false_easting :

    500000.0

    false_northing :

    0.0

    scale_factor_at_central_meridian :

    0.9996

    GeoTransform :

    475810 10 0 4426430 0 -10

    array(32613, dtype=int32)

  • time
    (time)
    datetime64[ns]
    2022-01-26T17:56:21.024000 ... 2...

    array(['2022-01-26T17:56:21.024000000', '2022-01-31T17:55:49.024000000',
           '2022-02-05T17:55:31.024000000', '2022-02-07T17:45:09.024000000',
           '2022-02-15T17:54:31.024000000', '2022-02-20T17:53:49.025000000',
           '2022-02-27T17:42:49.024000000', '2022-03-07T17:52:01.024000000',
           '2022-03-14T17:41:11.024000000', '2022-03-19T17:39:39.024000000',
           '2022-04-01T17:49:09.024000000', '2022-04-08T17:38:59.024000000',
           '2022-04-18T17:38:59.024000000', '2022-04-21T17:48:59.024000000',
           '2022-04-26T17:49:11.024000000', '2022-05-11T17:49:09.024000000',
           '2022-05-18T17:38:59.024000000', '2022-05-26T17:49:21.024000000',
           '2022-06-02T17:39:11.024000000', '2022-06-07T17:39:09.024000000',
           '2022-06-15T17:49:21.024000000', '2022-06-17T17:39:09.024000000',
           '2022-06-20T17:49:09.024000000', '2022-06-27T17:39:09.024000000',
           '2022-07-02T17:39:21.024000000', '2022-07-07T17:39:09.024000000',
           '2022-07-12T17:39:21.024000000', '2022-07-17T17:39:09.024000000',
           '2022-07-22T17:39:11.024000000', '2022-07-25T17:49:21.024000000',
           '2022-07-27T17:39:09.024000000', '2022-07-30T17:49:09.024000000',
           '2022-08-04T17:49:21.024000000', '2022-08-06T17:39:09.024000000',
           '2022-08-09T17:49:09.024000000', '2022-08-11T17:39:21.024000000',
           '2022-08-21T17:39:21.024000000', '2022-08-24T17:49:21.024000000',
           '2022-08-29T17:49:09.024000000', '2022-08-31T17:39:21.024000000',
           '2022-09-03T17:49:21.024000000', '2022-09-05T17:39:09.024000000',
           '2022-09-08T17:49:09.024000000', '2022-09-18T17:50:09.024000000',
           '2022-09-23T17:50:51.024000000', '2022-09-25T17:40:59.024000000',
           '2022-10-05T17:41:59.024000000', '2022-10-08T17:52:29.024000000',
           '2022-10-10T17:42:51.024000000', '2022-10-13T17:53:11.024000000',
           '2022-10-18T17:53:29.024000000', '2022-10-25T17:44:19.024000000',
           '2022-10-28T17:54:39.024000000', '2022-10-30T17:44:51.024000000',
           '2022-11-02T17:55:11.024000000', '2022-11-07T17:55:39.024000000',
           '2022-11-09T17:45:51.024000000', '2022-11-19T17:46:31.024000000',
           '2022-11-22T17:56:41.024000000', '2022-12-04T17:47:19.024000000',
           '2022-12-09T17:47:31.024000000', '2022-12-17T17:57:39.024000000',
           '2022-12-22T17:57:41.024000000', '2022-12-29T17:47:41.024000000'],
          dtype='datetime64[ns]')

  • band
    (band)
    object
    'B01' 'B02' 'B03' ... 'B12' 'B8A'

    array(['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B11',
           'B12', 'B8A'], dtype=object)

• Indexes: (4)
  • y
    PandasIndex

    PandasIndex(Index([4426425.0, 4426415.0, 4426405.0, 4426395.0, 4426385.0, 4426375.0,
           4426365.0, 4426355.0, 4426345.0, 4426335.0,
           ...
           4424865.0, 4424855.0, 4424845.0, 4424835.0, 4424825.0, 4424815.0,
           4424805.0, 4424795.0, 4424785.0, 4424775.0],
          dtype='float64', name='y', length=166))

  • x
    PandasIndex

    PandasIndex(Index([475815.0, 475825.0, 475835.0, 475845.0, 475855.0, 475865.0, 475875.0,
           475885.0, 475895.0, 475905.0,
           ...
           477155.0, 477165.0, 477175.0, 477185.0, 477195.0, 477205.0, 477215.0,
           477225.0, 477235.0, 477245.0],
          dtype='float64', name='x', length=144))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2022-01-26 17:56:21.024000', '2022-01-31 17:55:49.024000',
                   '2022-02-05 17:55:31.024000', '2022-02-07 17:45:09.024000',
                   '2022-02-15 17:54:31.024000', '2022-02-20 17:53:49.025000',
                   '2022-02-27 17:42:49.024000', '2022-03-07 17:52:01.024000',
                   '2022-03-14 17:41:11.024000', '2022-03-19 17:39:39.024000',
                   '2022-04-01 17:49:09.024000', '2022-04-08 17:38:59.024000',
                   '2022-04-18 17:38:59.024000', '2022-04-21 17:48:59.024000',
                   '2022-04-26 17:49:11.024000', '2022-05-11 17:49:09.024000',
                   '2022-05-18 17:38:59.024000', '2022-05-26 17:49:21.024000',
                   '2022-06-02 17:39:11.024000', '2022-06-07 17:39:09.024000',
                   '2022-06-15 17:49:21.024000', '2022-06-17 17:39:09.024000',
                   '2022-06-20 17:49:09.024000', '2022-06-27 17:39:09.024000',
                   '2022-07-02 17:39:21.024000', '2022-07-07 17:39:09.024000',
                   '2022-07-12 17:39:21.024000', '2022-07-17 17:39:09.024000',
                   '2022-07-22 17:39:11.024000', '2022-07-25 17:49:21.024000',
                   '2022-07-27 17:39:09.024000', '2022-07-30 17:49:09.024000',
                   '2022-08-04 17:49:21.024000', '2022-08-06 17:39:09.024000',
                   '2022-08-09 17:49:09.024000', '2022-08-11 17:39:21.024000',
                   '2022-08-21 17:39:21.024000', '2022-08-24 17:49:21.024000',
                   '2022-08-29 17:49:09.024000', '2022-08-31 17:39:21.024000',
                   '2022-09-03 17:49:21.024000', '2022-09-05 17:39:09.024000',
                   '2022-09-08 17:49:09.024000', '2022-09-18 17:50:09.024000',
                   '2022-09-23 17:50:51.024000', '2022-09-25 17:40:59.024000',
                   '2022-10-05 17:41:59.024000', '2022-10-08 17:52:29.024000',
                   '2022-10-10 17:42:51.024000', '2022-10-13 17:53:11.024000',
                   '2022-10-18 17:53:29.024000', '2022-10-25 17:44:19.024000',
                   '2022-10-28 17:54:39.024000', '2022-10-30 17:44:51.024000',
                   '2022-11-02 17:55:11.024000', '2022-11-07 17:55:39.024000',
                   '2022-11-09 17:45:51.024000', '2022-11-19 17:46:31.024000',
                   '2022-11-22 17:56:41.024000', '2022-12-04 17:47:19.024000',
                   '2022-12-09 17:47:31.024000', '2022-12-17 17:57:39.024000',
                   '2022-12-22 17:57:41.024000', '2022-12-29 17:47:41.024000'],
                  dtype='datetime64[ns]', name='time', freq=None))

  • band
    PandasIndex

    PandasIndex(Index(['B01', 'B02', 'B03', 'B04', 'B05', 'B06', 'B07', 'B08', 'B09', 'B11',
           'B12', 'B8A'],
          dtype='object', name='band'))

• Attributes: (0)

Now that we have a harmonized dataset, we still need to process the data as follows:

• Sentinel-2 L2A provides the Surface Reflectance (SR) data, which usually ranges from 0 (no reflection) to 1.0 (complete reflection). However, the actual values in the loaded dataset ranges from 0 to ~10,000. These data values need to be scaled to 0.0-1.0 by dividing the data by 10,000. More details can be found in section 2.3.10 of this document.

We will then explicitly trigger the dask computation using the compute() method and load the result into memory. This is to reduce repeated calls to retrieve data from MSPC. By loading the processed This wouldn’t have been possible if the dataset was large.

da = da / 1e4   # Scale data values from 0:10000 to 0:1.0
da = da / da.max(dim='band')  # additionally scale from 0-max -> 0-1 for visual quality
da = da.compute()

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

/home/runner/miniconda3/envs/interactive-sentinel-2-cookbook-dev/lib/python3.12/site-packages/dask/core.py:127: RuntimeWarning: invalid value encountered in divide
  return func(*(_execute_task(a, cache) for a in args))

We have now processed the data so that we can visualize it!

Let’s look at the Blue, Green and Red bands.

def plot_band(band, cmap):
    return da.sel(band=band).isel(time=0).hvplot(
        x='x', y='y', data_aspect=1, 
        cmap=cmap, geo=True, tiles='ESRI', 
        crs=crs.epsg(items[0].properties['proj:epsg']), rasterize=True,
        title=f"band: {band}, cmap: {cmap}",
        clabel='surface reflectance [0.0-1.0]'
    ).opts(
        frame_width=300,
        xlabel='longitude',
        ylabel='latitude',
        
    )

(plot_band('B04', 'Blues') + plot_band('B03', 'Greens') + plot_band('B02', 'Reds')).cols(2)

Let us make a dashboard composed of 4 different interactive plots showing the RGB view of the satellite observation for four different seasons. We need a function that will take a time input and does the following tasks:

1. plot an interactive RGB image of the data and overlay it on a map of the world.

2. provide a date slider widget which can be used to interact with the plot.

3. only set the default value of the date slider to the time, but allow the user to slide through the length of the entire dataset.

Using this function, we will be able to compose the dashboard.

def rgb_during(time):
    season_names = {
        1: 'Winter',
        2: 'Spring',
        3: 'Summer',
        4: 'Fall'
    }
    da_rgb = da.sel(band=['B04', 'B03', 'B02'])
    start_date = pd.to_datetime(da_rgb['time'].min().data).to_pydatetime()
    end_date = pd.to_datetime(da_rgb['time'].max().data).to_pydatetime()
    closest_date = pd.to_datetime(da_rgb.sel(time=time, method='nearest').time.data).to_pydatetime()
    dt_slider = pnw.DateSlider(name='Date', start=start_date, end=end_date, value=closest_date)
    
    def get_obs_on(t):
        season_key = [month%12 // 3 + 1 for month in range(1, 13)][t.month-1]
        season = season_names[season_key]
        return da.sel(band=['B04', 'B03', 'B02']).sel(time=t, method='nearest').transpose('y', 'x', 'band').hvplot.rgb(
            x='x', y='y', bands='band', 
            geo=True, tiles='ESRI', crs=crs.epsg(items[0].properties['proj:epsg']), 
            rasterize=True, title=f"{season}: {t.strftime('%Y-%m-%d')}"
        ).opts(
            frame_width=300,
            xlabel='longitude',
            ylabel='latitude',
        )
        
    
    return pn.panel(pn.Column(
                pn.bind(get_obs_on, t=dt_slider), 
                pn.Row(
                    pn.Spacer(width=60),
                    dt_slider,
                )
            ))

rgb_during('2022-02-01')

Let’s now compose a dashboard using panel.

winter = '2022-01-15'
spring = '2022-04-30'
summer = '2022-08-01'
fall = '2022-09-15'

winter_plot = rgb_during(winter)
spring_plot = rgb_during(spring)
summer_plot = rgb_during(summer)
fall_plot = rgb_during(fall)

pn.Column(
    pn.Row(winter_plot, spring_plot),
    pn.Row(summer_plot, fall_plot)
)

---

Summary

In this recipe, we looked at how to access Sentinel-2 satellite data over a region of interest and create an interactive dashboard to visualize the data.

What’s next?

We plotted the RGB or True color image of our region of interest using a subset of all the bands available. We can further calculate useful indices, such as the Normalized Difference Snow Index (NDSI) or the Normalized Difference Vegetation Index.

Resources and references

• Authored by Pritam Das (@pritamd47), June 2023 during Project Pythia cookoff 2023.

• This notebook takes a lot of inspiration from the Landsat ML Cookbook by Demetris Roumis.

• This notebook uses concepts and code illustrated in the Reading Data from the STAC API.