Time series of vegetation indices#

Short description

This notebook demonstrates a general usecase for the SpaceSense library.

In this notebook, you will search for, select, and obtain Sentinel-2 data and vegetation indices, and performa simple time series analysis.

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

[1]:

from spacesense import Client
import matplotlib.pyplot as plt
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 an agricultural field in Sardinia, can be downloaded here

[2]:

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

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

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

3 - Search S2#

[3]:

start_date = "2021-01-01"
end_date = "2022-01-01"
res_S2 = client.s2_search(aoi=aoi, start_date=start_date, end_date=end_date, query_filters = {"valid_pixel_percentage" : {">=": 50}})
res_S2.filter_duplicate_dates()
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
80 S2B_32TMK_20210108_0_L2A 2021-01-08 32TMK 50.73 sentinel-2b 022 S2B_MSIL2A_20210108T101319_N0214_R022_T32TMK_2... 2021-01-08T10:19:47Z 100.0 0.0 0.00 0.00 0.00 0.0 49.27 0.0 0.0 0.0
79 S2A_32TMK_20210116_0_L2A 2021-01-16 32TMK 100.00 sentinel-2a 065 S2A_MSIL2A_20210116T102351_N0214_R065_T32TMK_2... 2021-01-16T10:29:44Z 100.0 0.0 0.00 97.56 0.00 0.0 0.00 0.0 0.0 0.0
78 S2B_32TMK_20210118_0_L2A 2021-01-18 32TMK 100.00 sentinel-2b 022 S2B_MSIL2A_20210118T101249_N0214_R022_T32TMK_2... 2021-01-18T10:19:47Z 100.0 0.0 0.00 96.59 0.49 0.0 0.00 0.0 0.0 0.0
77 S2A_32TMK_20210123_0_L2A 2021-01-23 32TMK 100.00 sentinel-2a 022 S2A_MSIL2A_20210123T101321_N0214_R022_T32TMK_2... 2021-01-23T10:19:48Z 100.0 0.0 0.00 97.40 0.00 0.0 0.00 0.0 0.0 0.0
76 S2A_32TMK_20210126_0_L2A 2021-01-26 32TMK 100.00 sentinel-2a 065 S2A_MSIL2A_20210126T102311_N0214_R065_T32TMK_2... 2021-01-26T10:29:44Z 100.0 0.0 0.00 97.40 0.00 0.0 0.00 0.0 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4 S2B_32TMK_20211127_0_L2A 2021-11-27 32TMK 89.81 sentinel-2b 065 S2B_MSIL2A_20211127T102259_N0301_R065_T32TMK_2... 2021-11-27T10:29:39Z 100.0 0.0 10.19 49.81 7.91 0.0 0.00 0.0 0.0 0.0
3 S2A_32TMK_20211212_0_L2A 2021-12-12 32TMK 100.00 sentinel-2a 065 S2A_MSIL2A_20211212T102431_N0301_R065_T32TMK_2... 2021-12-12T10:29:43Z 100.0 0.0 0.00 64.09 7.70 0.0 0.00 0.0 0.0 0.0
2 S2B_32TMK_20211217_0_L2A 2021-12-17 32TMK 100.00 sentinel-2b 065 S2B_MSIL2A_20211217T102329_N0301_R065_T32TMK_2... 2021-12-17T10:29:37Z 100.0 0.0 0.00 65.58 18.54 0.0 0.00 0.0 0.0 0.0
1 S2A_32TMK_20211219_0_L2A 2021-12-19 32TMK 100.00 sentinel-2a 022 S2A_MSIL2A_20211219T101431_N0301_R022_T32TMK_2... 2021-12-19T10:19:48Z 100.0 0.0 0.00 67.48 10.16 0.0 0.00 0.0 0.0 0.0
0 S2B_32TMK_20211224_0_L2A 2021-12-24 32TMK 94.91 sentinel-2b 022 S2B_MSIL2A_20211224T101329_N0301_R022_T32TMK_2... 2021-12-24T10:19:43Z 100.0 0.0 5.09 44.01 0.00 0.0 0.00 0.0 0.0 0.0

81 rows × 18 columns

[4]:

print(res_S2.item_collection.calendar())

                              2021

      January               February               March
Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su
             1  2  3   1  2  3  4  5  6  7   1  2  3  4  5  6  7
 4  5  6  7  8  9 10   8  9 10 11 12 13 14   8  9 10 11 12 13 14
11 12 13 14 15 16 17  15 16 17 18 19 20 21  15 16 17 18 19 20 21
18 19 20 21 22 23 24  22 23 24 25 26 27 28  22 23 24 25 26 27 28
25 26 27 28 29 30 31                        29 30 31

       April                  May                   June
Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su
          1  2  3  4                  1  2      1  2  3  4  5  6
 5  6  7  8  9 10 11   3  4  5  6  7  8  9   7  8  9 10 11 12 13
12 13 14 15 16 17 18  10 11 12 13 14 15 16  14 15 16 17 18 19 20
19 20 21 22 23 24 25  17 18 19 20 21 22 23  21 22 23 24 25 26 27
26 27 28 29 30        24 25 26 27 28 29 30  28 29 30

        July                 August              September
Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su
          1  2  3  4                     1         1  2  3  4  5
 5  6  7  8  9 10 11   2  3  4  5  6  7  8   6  7  8  9 10 11 12
12 13 14 15 16 17 18   9 10 11 12 13 14 15  13 14 15 16 17 18 19
19 20 21 22 23 24 25  16 17 18 19 20 21 22  20 21 22 23 24 25 26
26 27 28 29 30 31     23 24 25 26 27 28 29  27 28 29 30

      October               November              December
Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su  Mo Tu We Th Fr Sa Su
             1  2  3   1  2  3  4  5  6  7         1  2  3  4  5
 4  5  6  7  8  9 10   8  9 10 11 12 13 14   6  7  8  9 10 11 12
11 12 13 14 15 16 17  15 16 17 18 19 20 21  13 14 15 16 17 18 19
18 19 20 21 22 23 24  22 23 24 25 26 27 28  20 21 22 23 24 25 26
25 26 27 28 29 30 31  29 30                 27 28 29 30 31

sentinel-2a (1)
sentinel-2b (1)
81 total dates

4 - Specify bands and pre-calculated indices#

Here, we select only the blue, red, and near infrared (NIR) bands of Sentinel-2, as well as the vegetation indices of NDVI and LAI

[5]:

res_S2.output_bands = ["B02","B04","B08","NDVI","LAI"]

5 - Filter search results#

Select only the results which you require. Learn more about how to filter and manipulate search results with this example.

In this example, we wish to use data with little cloud interference to see the most accurate representation of vegetation indices. The initial search took only those S2 images with 50% or more “valid_pixel_percentage”. This initial search can be used to quickly see how many images, in general, are available for the AOI and TOI. Further filtering and refinement can be done later depending on those results, like the one performed below. Here we see there are plenty of available images, so we can further filter them based on “valid_pixel_percentage”, to attempt to get more cloud-free data.

It is important to note that the Sentinel cloud detection is not perfect, so some exploration and visualziation of the data may be needed to validate the images.

[6]:

res_S2.dataframe = res_S2.dataframe[res_S2.dataframe["valid_pixel_percentage"] > 80]
res_S2.dataframe

[6]:
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
79 S2A_32TMK_20210116_0_L2A 2021-01-16 32TMK 100.00 sentinel-2a 065 S2A_MSIL2A_20210116T102351_N0214_R065_T32TMK_2... 2021-01-16T10:29:44Z 100.0 0.0 0.00 97.56 0.00 0.0 0.0 0.0 0.0 0.0
78 S2B_32TMK_20210118_0_L2A 2021-01-18 32TMK 100.00 sentinel-2b 022 S2B_MSIL2A_20210118T101249_N0214_R022_T32TMK_2... 2021-01-18T10:19:47Z 100.0 0.0 0.00 96.59 0.49 0.0 0.0 0.0 0.0 0.0
77 S2A_32TMK_20210123_0_L2A 2021-01-23 32TMK 100.00 sentinel-2a 022 S2A_MSIL2A_20210123T101321_N0214_R022_T32TMK_2... 2021-01-23T10:19:48Z 100.0 0.0 0.00 97.40 0.00 0.0 0.0 0.0 0.0 0.0
76 S2A_32TMK_20210126_0_L2A 2021-01-26 32TMK 100.00 sentinel-2a 065 S2A_MSIL2A_20210126T102311_N0214_R065_T32TMK_2... 2021-01-26T10:29:44Z 100.0 0.0 0.00 97.40 0.00 0.0 0.0 0.0 0.0 0.0
75 S2A_32TMK_20210205_0_L2A 2021-02-05 32TMK 100.00 sentinel-2a 065 S2A_MSIL2A_20210205T102221_N0214_R065_T32TMK_2... 2021-02-05T10:29:43Z 100.0 0.0 0.00 97.40 0.16 0.0 0.0 0.0 0.0 0.0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
4 S2B_32TMK_20211127_0_L2A 2021-11-27 32TMK 89.81 sentinel-2b 065 S2B_MSIL2A_20211127T102259_N0301_R065_T32TMK_2... 2021-11-27T10:29:39Z 100.0 0.0 10.19 49.81 7.91 0.0 0.0 0.0 0.0 0.0
3 S2A_32TMK_20211212_0_L2A 2021-12-12 32TMK 100.00 sentinel-2a 065 S2A_MSIL2A_20211212T102431_N0301_R065_T32TMK_2... 2021-12-12T10:29:43Z 100.0 0.0 0.00 64.09 7.70 0.0 0.0 0.0 0.0 0.0
2 S2B_32TMK_20211217_0_L2A 2021-12-17 32TMK 100.00 sentinel-2b 065 S2B_MSIL2A_20211217T102329_N0301_R065_T32TMK_2... 2021-12-17T10:29:37Z 100.0 0.0 0.00 65.58 18.54 0.0 0.0 0.0 0.0 0.0
1 S2A_32TMK_20211219_0_L2A 2021-12-19 32TMK 100.00 sentinel-2a 022 S2A_MSIL2A_20211219T101431_N0301_R022_T32TMK_2... 2021-12-19T10:19:48Z 100.0 0.0 0.00 67.48 10.16 0.0 0.0 0.0 0.0 0.0
0 S2B_32TMK_20211224_0_L2A 2021-12-24 32TMK 94.91 sentinel-2b 022 S2B_MSIL2A_20211224T101329_N0301_R022_T32TMK_2... 2021-12-24T10:19:43Z 100.0 0.0 5.09 44.01 0.00 0.0 0.0 0.0 0.0 0.0

79 rows × 18 columns

6 - 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

[7]:

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"}
    )

2022-11-07 16:17:37,288 INFO spacesense.core : created everything

7 - See the fused dataset#

All the S2 images that we selected are now fused into the same file, and can be easily accessed via date or variable. `Click here <>`__ to learn more about how to use the resulting Xarray Datasets

[8]:

ds = fusion.dataset
ds

[8]:

<xarray.Dataset>
Dimensions:  (time: 79, y: 43, x: 88)
Coordinates:
  * time     (time) datetime64[ns] 2021-01-16 2021-01-18 ... 2021-12-24
  * y        (y) float32 39.83 39.83 39.83 39.83 ... 39.83 39.83 39.83 39.83
  * x        (x) float64 8.623 8.623 8.623 8.623 8.623 ... 8.63 8.63 8.63 8.63
Data variables:
    S2_B02   (time, y, x) float32 0.0 0.0123 0.0144 0.0126 ... 0.0247 0.0265 0.0
    S2_B04   (time, y, x) float32 0.0 0.0208 0.0214 0.021 ... 0.0455 0.0449 0.0
    S2_B08   (time, y, x) float32 0.0 0.1198 0.1394 0.1282 ... 0.1478 0.1556 0.0
    S2_NDVI  (time, y, x) float32 nan 0.7041 0.7338 0.7185 ... 0.5292 0.5521 nan
    S2_LAI   (time, y, x) float32 0.4768 1.603 1.801 ... 1.773 1.856 0.4768
Attributes:
    crs:              +init=epsg:4326
    transform:        [ 8.98360017e-05  0.00000000e+00  8.62245945e+00  0.000...
    res:              [8.98360017e-05 9.04414549e-05]
    ulx, uly:         [ 8.62245945 39.83104984]
    s2_data_lineage:  {"Data origin": "S3 bucket (ARN=arn:aws:s3:::sentinel-c...
    additional_info:  {'Info': 'You can put any additional information, such ...

8 - Plot timeseries of vegetation indices with matplotlib#

First, process the data into time series format

[9]:

ds = fusion.dataset

# Take the field level mean (average over latitude and longitude) of the vegetation indices and remove N/A
ndvi_ts = ds.S2_NDVI.mean(dim={'y','x'}).dropna(dim='time')
lai_ts = ds.S2_LAI.mean(dim={'y','x'}).dropna(dim='time')

Then plot

[10]:

# Set up a time series plot
fig, ax = plt.subplots(figsize=(12,6))
ax2 = ax.twinx()

# Plot both NDVI and LAI on the same plot, with different y axes
ndvi_ts.plot(color='blue',label="NDVI",ax=ax)
lai_ts.plot(color='green',label="LAI",ax=ax2)
fig.legend()
plt.title("Veg Indices Time Series")

[10]:

Text(0.5, 1.0, 'Veg Indices Time Series')
../_images/notebooks_6_timeseries_veg_index_23_1.png

As you can see, this time series gives a good depiction of the changes in vegetation over a year. However, we do still see some errors, likely due to cloud coverage not correctly classified from the Sentinel scene classification.