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"]

5 - Obtain S2 data through Fuse function#

We pass the filtered S2 search result object to the “fuse” function in order to download the organize the S2 result

[5]:
fusion = client.fuse(
    catalogs_list=[res_S2],
    additional_info={"Info": "You can put any additional information, such as custom non-raster or non-vector data, in this parameter"}
    )
fusion.dataset
2022-09-14 14:45:02,374 INFO spacesense.core : created everything
[5]:
<xarray.Dataset>
Dimensions:  (time: 1, y: 2357, x: 327)
Coordinates:
  * time     (time) datetime64[ns] 2021-05-01
  * y        (y) float32 -2.487 -2.487 -2.487 -2.487 ... -2.504 -2.504 -2.504
  * x        (x) float64 121.3 121.3 121.3 121.3 ... 121.3 121.3 121.3 121.3
Data variables:
    s2_B01   (time, y, x) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0477 0.0477 0.0 0.0
    s2_B02   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0401 0.0442 0.0442 0.0
    s2_B03   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0257 0.0273 0.0263 0.0
    s2_B04   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0139 0.0128 0.0135 0.0
    s2_B08   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0071 0.0091 0.0096 0.0
    s2_SCL   (time, y, x) float32 0.0 0.0 0.0 0.0 0.0 ... 6.0 6.0 6.0 6.0 6.0
Attributes:
    crs:              +init=epsg:4326
    transform:        [ 8.98315284e-05  0.00000000e+00  1.21303645e+02  0.000...
    res:              [8.98315284e-05 7.00728646e-06]
    ulx, uly:         [121.30364507  -2.48706115]
    additional_info:  {'Info': 'You can put any additional information, such ...

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
[6]:
<xarray.Dataset>
Dimensions:  (time: 1, y: 2357, x: 327)
Coordinates:
  * time     (time) datetime64[ns] 2021-05-01
  * y        (y) float32 -2.487 -2.487 -2.487 -2.487 ... -2.504 -2.504 -2.504
  * x        (x) float64 121.3 121.3 121.3 121.3 ... 121.3 121.3 121.3 121.3
Data variables:
    s2_B01   (time, y, x) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0477 0.0477 0.0 0.0
    s2_B02   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0401 0.0442 0.0442 0.0
    s2_B03   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0257 0.0273 0.0263 0.0
    s2_B04   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0139 0.0128 0.0135 0.0
    s2_B08   (time, y, x) float32 0.0 0.0 0.0 0.0 ... 0.0071 0.0091 0.0096 0.0
    s2_SCL   (time, y, x) float32 0.0 0.0 0.0 0.0 0.0 ... 6.0 6.0 6.0 6.0 6.0
    Chl_a    (time, y, x) float32 nan nan nan nan nan ... 0.3699 0.4693 nan nan
    SPM      (time, y, x) float32 0.0 0.0 0.0 0.0 0.0 ... 6.8 9.212 9.835 0.0
Attributes:
    crs:              +init=epsg:4326
    transform:        [ 8.98315284e-05  0.00000000e+00  1.21303645e+02  0.000...
    res:              [8.98315284e-05 7.00728646e-06]
    ulx, uly:         [121.30364507  -2.48706115]
    additional_info:  {'Info': 'You can put any additional information, such ...

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 0x7f1e440a7a50>
../_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 0x7f1e440e9150>
../_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 0x7f1e3ecb54d0>
../_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 0x7f1e3eccea90>
../_images/notebooks_8_create_custom_index_and_masking_23_1.png