Create custom water or vegetation index
Contents
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-11-07 16:16:54,385 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] s2_data_lineage: {"Data origin": "S3 bucket (ARN=arn:aws:s3:::sentinel-c... 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:
[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] s2_data_lineage: {"Data origin": "S3 bucket (ARN=arn:aws:s3:::sentinel-c... 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 0x7f07ece81dd0>

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>

Plot the newly calculated Chlophyll index
[9]:
ds.CHL_A.plot(vmin=0,vmax=5)
[9]:
<matplotlib.collections.QuadMesh at 0x7f07ec442790>

Plot the newly calculated SPM index
[10]:
ds.SPM.plot(vmin=0,vmax=50)
[10]:
<matplotlib.collections.QuadMesh at 0x7f07ec32d450>
