Create a spectral plot for Sentinel 2a pixels

What does this notebook do? This notebook allows you to plot up pixel spectra at the bands’ native resolutions. The user loads in some data, then selects a pixel for interrogation using a plot widget. The spectra for this pixel is extracted and plotted. This notebook takes advantage of the different native resolutions of bands from the Sentinel 2 satellite to explore the variability associated with bands being recorded at different spatial resolutions. The code uses the 60m resolution of the nbar_coastal_aerosol band to determine a bounding box for the higher resolution bands. The pixels that fall within that bounding box are selected, and the range of spectral values across the additional pixels are plotted to characterise the spectral heterogeneity of the selected area.

Before you run this notebook: You need to run the following commands from the command line prior to launching jupyter notebook from the same terminal so that the required libraries and paths are set.

module use /g/data/v10/public/modules/modulefiles

module load dea

Date: September 2018.

Author: Claire Krause

Tags: Sentinel2, products, DEAPlotting, three_band_image, query, spectra, ipywidgets, Scripts, dc.load, on_click

About Sentinel 2 bands

The spatial resolution of Sentinel 2 is dependent on the band, with different bands being collected at 10m, 20m or 60m resolution.

** 10m spatial resolution ** 10m spatial resolution

Sentinel 2 bands

DEA band name

Band number

Central wavelength (nm)

Resolution (m)

Bandwidth (nm)

Blue

nbar(t)_blue

2

490

10

65

Green

nbar(t)_green

3

560

10

35

Red

nbar(t)_red

4

665

10

30

NIR

nbar(t)_nir_1

8

842

10

115

** 20m spatial resolution ** 20m spatial resolution

Sentinel 2 bands

DEA band name

Band number

Central wavelength (nm)

Resolution (m)

Bandwidth (nm)

Vegetation red edge

nbar(t)_red_edge_1

5

705

20

15

Vegetation red edge

nbar(t)_red_edge_2

6

740

20

15

Vegetation red edge

nbar(t)_red_edge_3

7

783

20

20

Narrow NIR

nbar(t)_nir_2

8A

865

20

20

SWIR

nbar(t)_swir_2

11

1610

20

90

SWIR

nbar(t)_swir_3

12

2190

20

180

** 60m spatial resolution ** 60m spatial resolution

| Sentinel 2 bands | DEA band name | Band number | Central wavelength (nm) | Resolution (m) | Bandwidth (nm) | | —————–|—————|————-|————————-|—————-|—————-| | Coastal aerosol | nbar(t)_coastal_aerosol | 1 | 443 | 60 | 20 | N.B. Bands 9 and 10 are not available within DEA

Images from ESA

Import the required libraries

[1]:
%pylab notebook

from datacube import Datacube

# Import widgets for interactive notebook
from ipywidgets import interact, fixed
import ipywidgets as widgets

# Import the custom script for plotting. This script can be found in the dea-notebooks repository.
import sys
import os.path
sys.path.append(os.path.expanduser('../10_Scripts'))
import DEAPlotting

# The path to the Sentinel 2 development DEA instance configuration file
dc = Datacube(app='Sentinel2 spectra')
Populating the interactive namespace from numpy and matplotlib

Define a function that does all of the work

We are putting this all into a single widget function so that once a pixel is selected, the user does not need to run additional cells to get to the desired outputs

[2]:
def onclick(event):
    '''
    This function performs all of the work in this notebook. We have put all of the
    processing within this function so that all the calculations are done following
    a 'click' and therefore do not need to be in separate cells. This makes the notebook
    smoother and minimises the cells that need to be manually run by the user once a
    location has been selected.

    This particular widget function uses the selected to location to plot two images:
    - spectra for the pixel closest to the chosen location
    - spectra for the pixel/s closest to the chosen location using the 60m pixel as
    a bounding box for the higher resolution pixels.
    These two images are run in the subsequent cells.
    '''
    global pixelx, pixely, spectra, spectramin, spectramax, spectramean
    pixelx, pixely = int(event.xdata), int(event.ydata)
    w.value = 'pixelx : {}, pixely : {}'.format(pixelx, pixely)
    plt.plot(pixelx, pixely, 'ro', markersize=5)

    # Find the pixel closest to the chosen pixel at each resolution
    Pixel10m = Resolution10m.sel(y=pixely, x=pixelx, method='nearest')
    Pixel20m = Resolution20m.sel(y=pixely, x=pixelx, method='nearest')
    Pixel60m = Resolution60m.sel(y=pixely, x=pixelx, method='nearest')

    # Grab the pixel spectral values for each band
    spectra = [Pixel60m.nbar_coastal_aerosol.isel(time=mytime).values,
               Pixel10m.nbar_blue.isel(time=mytime).values,
               Pixel10m.nbar_green.isel(time=mytime).values,
               Pixel10m.nbar_red.isel(time=mytime).values,
               Pixel20m.nbar_red_edge_1.isel(time=mytime).values,
               Pixel20m.nbar_red_edge_2.isel(time=mytime).values,
               Pixel20m.nbar_red_edge_3.isel(time=mytime).values,
               Pixel10m.nbar_nir_1.isel(time=mytime).values,
               Pixel20m.nbar_nir_2.isel(time=mytime).values,
               Pixel20m.nbar_swir_2.isel(time=mytime).values,
               Pixel20m.nbar_swir_3.isel(time=mytime).values,
               ]

    # Get the location of the selected pixel at the coursest resolution
    Pixel60m = Resolution60m.sel(y=pixely, x=pixelx, method='nearest')

    # Find the index locations of the lat/lon of that pixel
    xindex = Resolution60m.indexes['x'].get_loc(
        Pixel60m.x.values.item(), method='nearest')
    yindex = Resolution60m.indexes['y'].get_loc(
        Pixel60m.y.values.item(), method='nearest')

    # Get the index for the pixels next to the chosen pixel
    xmax = Resolution60m.x.isel(x=xindex+1)
    xmin = Resolution60m.x.isel(x=xindex-1)
    ymax = Resolution60m.y.isel(y=yindex-1)
    ymin = Resolution60m.y.isel(y=yindex+1)

    # Now work out what the lat/lon is for halfway between pixel +-1 (to keep our resolution at 60 x 60 m)
    latmin = mean((xmin, pixelx))
    latmax = mean((pixelx, xmax))
    lonmin = mean((ymin, pixely))
    lonmax = mean((pixely, ymax))

    # Grab all of the pixels that fall within the 60m pixel bounds
    bluepixels = Resolution10m.nbar_blue.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    greenpixels = Resolution10m.nbar_green.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    redpixels = Resolution10m.nbar_red.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    rededge1pixels = Resolution20m.nbar_red_edge_1.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    rededge2pixels = Resolution20m.nbar_red_edge_2.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    rededge3pixels = Resolution20m.nbar_red_edge_3.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    nir1pixels = Resolution10m.nbar_nir_1.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    nir2pixels = Resolution20m.nbar_nir_2.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    swir2pixels = Resolution20m.nbar_swir_2.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))
    swir3pixels = Resolution20m.nbar_swir_3.sel(
        x=slice(latmin, latmax), y=slice(lonmax, lonmin))

    # Grab the min, max and mean of the pixels within the 60m bounding box
    spectramin = [Pixel60m.nbar_coastal_aerosol.isel(time=mytime).min().item(),
                  bluepixels.isel(time=mytime).min().item(),
                  greenpixels.isel(time=mytime).min().item(),
                  redpixels.isel(time=mytime).min().item(),
                  rededge1pixels.isel(time=mytime).min().item(),
                  rededge2pixels.isel(time=mytime).min().item(),
                  rededge3pixels.isel(time=mytime).min().item(),
                  nir1pixels.isel(time=mytime).min().item(),
                  nir2pixels.isel(time=mytime).min().item(),
                  swir2pixels.isel(time=mytime).min().item(),
                  swir3pixels.isel(time=mytime).min().item(),
                  ]

    spectramax = [Pixel60m.nbar_coastal_aerosol.isel(time=mytime).max().item(),
                  bluepixels.isel(time=mytime).max().item(),
                  greenpixels.isel(time=mytime).max().item(),
                  redpixels.isel(time=mytime).max().item(),
                  rededge1pixels.isel(time=mytime).max().item(),
                  rededge2pixels.isel(time=mytime).max().item(),
                  rededge3pixels.isel(time=mytime).max().item(),
                  nir1pixels.isel(time=mytime).max().item(),
                  nir2pixels.isel(time=mytime).max().item(),
                  swir2pixels.isel(time=mytime).max().item(),
                  swir3pixels.isel(time=mytime).max().item(),
                  ]

    spectramean = [Pixel60m.nbar_coastal_aerosol.isel(time=mytime).mean().item(),
                   bluepixels.isel(time=mytime).mean().item(),
                   greenpixels.isel(time=mytime).mean().item(),
                   redpixels.isel(time=mytime).mean().item(),
                   rededge1pixels.isel(time=mytime).mean().item(),
                   rededge2pixels.isel(time=mytime).mean().item(),
                   rededge3pixels.isel(time=mytime).mean().item(),
                   nir1pixels.isel(time=mytime).mean().item(),
                   nir2pixels.isel(time=mytime).mean().item(),
                   swir2pixels.isel(time=mytime).mean().item(),
                   swir3pixels.isel(time=mytime).mean().item(),
                   ]

Set up the wavelengths and corresponding labels for later on

[3]:
wavelengths = [443, 490, 560, 665, 705, 740, 783, 842, 865, 1610, 2190]
wavelengthlabels = ['aerosol', 'blue', 'green', 'red', 'rededge1', 'rededge2',
                    'rededge3', 'nir1', 'nir2', 'swir2', 'swir3']

Set up the extraction queries and load the corresponding data

Grab the 10m resolution bands

[4]:
query = {
        'lat': (-35.25, -35.35),
        'lon': (149.05, 149.17),
        'output_crs': 'EPSG:3577',
        'resolution': (-10, 10),
        'time':('2017-01-01', '2017-02-15'),
        }

Resolution10m = dc.load(product='s2a_ard_granule', group_by='solar_day',
                        measurements = ['nbar_blue', 'nbar_green', 'nbar_red', 'nbar_nir_1'], **query)

Grab the 20m resolution bands

[5]:
query = {
        'lat': (-35.25, -35.35),
        'lon': (149.05, 149.17),
        'output_crs': 'EPSG:3577',
        'resolution': (-20, 20),
        'time':('2017-01-01', '2017-02-15'),
        }

Resolution20m = dc.load(product='s2a_ard_granule', group_by='solar_day',
                        measurements = ['nbar_red_edge_1', 'nbar_red_edge_2', 'nbar_red_edge_3',
                                        'nbar_nir_2', 'nbar_swir_2', 'nbar_swir_3'], **query)

Grab the 60m resolution aerosol band

[6]:
query = {
        'lat': (-35.25, -35.35),
        'lon': (149.05, 149.17),
        'output_crs': 'EPSG:3577',
        'resolution': (-60, 60),
        'time':('2017-01-01', '2017-02-15'),
        }

Resolution60m = dc.load(product='s2a_ard_granule', group_by='solar_day',
                        measurements = ['nbar_coastal_aerosol'], **query)

Find a time where there is clear image without clouds etc

This needs to be done by trial and error by changing the time here, and checking the output in the next plot cell.

[7]:
mytime = 1

Plot the extracted region, and select a pixel for spectral interrogation

[10]:
print('\033[1m' + 'Click on the pixel you would like to interrogate' + '\033[0m')

DEAPlotting.three_band_image(Resolution10m, bands = ['nbar_red', 'nbar_green', 'nbar_blue'],
                             time = mytime, contrast_enhance=True)

fig = plt.gcf()
w = widgets.HTML("Click on the pixel you would like to interrogate")

cid = fig.canvas.mpl_connect('button_press_event', onclick)
display(w)
Click on the pixel you would like to interrogate