Create custom water or vegetation index#

Short description

This notebook demonstrates a general usecase for the SpaceSense library.

In this notebook, you will search for, select, and obtain a Sentinel-2 image. We want the selected data needs to be cloud-free in order to better analyze the index. Custom indices can then be calculated over the area. This example examines several water quality indices.

1 - Import spacesense object(s) and other dependencies#

[1]:

from spacesense import Client
import matplotlib.pyplot as plt
import numpy as np
import json

import os
if "SS_API_KEY" not in os.environ:
    from getpass import getpass
    api_key = getpass('Enter your api key : ')
    os.environ["SS_API_KEY"] = api_key

2 - Define AOI and output options#

The GeoJSON of this AOI, containing a small portion of Lake Matano and coastline in Indonesia, can be downloaded here

[2]:

# Define the AOI
with open("../resources/aois/lake_matano.geojson", mode="r") as file:
    aoi = json.load(file)

# Get an instance of the SpaceSense Client object
client = Client(id="s2_custom_index")

# Enable to save data in local files
client.enable_local_output()

3 - Search S2#

[3]:

start_date = "2021-05-01"
end_date = "2021-05-02"
# Scenes should be cloud-free in order to best analyze the resulting indices
res_S2 = client.s2_search(aoi=aoi, start_date=start_date, end_date=end_date, query_filters = {"valid_pixel_percentage" : {">=": 95}})
res_S2.dataframe

[3]:

	id	date	tile	valid_pixel_percentage	platform	relative_orbit_number	product_id	datetime	swath_coverage_percentage	no_data	cloud_shadows	vegetation	not_vegetated	water	cloud_medium_probability	cloud_high_probability	thin_cirrus	snow
0	S2B_51MUT_20210501_0_L2A	2021-05-01	51MUT	100.0	sentinel-2b	060	S2B_MSIL2A_20210501T021339_N0300_R060_T51MUT_2...	2021-05-01T02:18:15Z	100.0	0.0	0.0	26.02	0.0	72.95	0.0	0.0	0.0	0.0

4 - Specify bands#

Only selecting bands from S2 that we are interested in. The water quality indices we want to calculate use bands 1, 2, 3, 4, and 8, thus we only choose those bands plus the scene classification.

[4]:

res_S2.output_bands = ["B01","B02","B03","B04","B08","SCL"]

6 - Postprocessing: Calculating indices#

Various water quality indices are calculated from the resulting S2 band data.

The following quantities and formulas are used:

Concentration of chlorophyll a (Chl_a) (source)
Suspended particle matter (SPM) (source)

[6]:

ds = fusion.dataset
ds = ds.assign(CHL_A = (4.23 * ((ds.S2_B03/ds.S2_B01) ** 3.94)))
ds = ds.assign(SPM = 2887 * (ds.S2_B08)**1.223)
### Necessary to replace infinite values with nan
ds['CHL_A'] = ds['CHL_A'].where(np.isinf(ds['CHL_A'])==False)
ds

7 - Visualization#

Plot RGB truecolor image

[7]:

# Define a function to normalize reflectance from S2 bands
def norm(band):
    band_min, band_max = band.min(), band.max()
    return ((band - band_min)/(band_max - band_min))

# As we only have one date, we want to simplify the dimensions
ds = ds.isel(time=0)

# Normalize and brighten (if necessary) the R, G, and B bands
brightness_factor = 1
blue = norm(ds.S2_B02) * brightness_factor
green = norm(ds.S2_B03) * brightness_factor
red = norm(ds.S2_B04) * brightness_factor

# Using numpy, stack the RGB bands into a single object, and plot using matplotlib
rgb = np.dstack([red,green,blue])
plt.figure(figsize=(6,6))
# When plotting with imshow, it is sometimes necessary to directly assign the plot extent to the
# Min and max lat and lon values, as shown here. Otherwise, plots may be auto-scaled and skewed.
plt.imshow(rgb, extent=[float(ds.x.min().values),float(ds.x.max().values),float(ds.y.min().values),float(ds.y.max().values)])

[7]:

<matplotlib.image.AxesImage at 0x7f07ece81dd0>

../_images/notebooks_8_create_custom_index_and_masking_17_1.png

Plot the scene classifications according to the S2 SCL band. We can clearly see the yellow corresponding to 6 is water

[8]:

# Plot the values of the SCL band
ds.S2_SCL.plot()

[8]:

<matplotlib.collections.QuadMesh at 0x7f07ed1d71d0>

../_images/notebooks_8_create_custom_index_and_masking_19_1.png

Plot the newly calculated Chlophyll index

[9]:

ds.CHL_A.plot(vmin=0,vmax=5)

[9]:

<matplotlib.collections.QuadMesh at 0x7f07ec442790>

../_images/notebooks_8_create_custom_index_and_masking_21_1.png

Plot the newly calculated SPM index

[10]:

ds.SPM.plot(vmin=0,vmax=50)

[10]:

<matplotlib.collections.QuadMesh at 0x7f07ec32d450>

../_images/notebooks_8_create_custom_index_and_masking_23_1.png

Create custom water or vegetation index

Contents