Search and fuse S1, S2 and weather data with custom raster or vector#

Short description

This notebook introduces the Fuse functionality in the SpaceSense library.

In this notebook, you will search for, filter, and fuse Sentinel-1, Sentinel-2 and ERA5 weather data with user-provided raster and vector files.

0 - Install Spacesense Client Library and dependencies#

Please follow the installation process and use the virtual environment in this notebook.

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

[1]:

from spacesense import Client, Raster, Vector
import matplotlib.pyplot as plt
import numpy as np
import json
import datetime
iso = datetime.date.fromisoformat

2 - Configure the API Key by setting the SS_API_KEY environment variable#

[2]:

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

3 - Define AOI and output options#

The GeoJSON of this AOI, containing an agricultural field in Sardinia, can be downloaded here

[3]:

# 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="fuse_with_custom_data")

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

4 - Search S1#

[4]:

start_date = "2021-01-01"
end_date = "2021-01-02"
# Retrieves all S1 images corresponding to the aoi, start date, and end date
res_S1 = client.s1_search(aoi=aoi, start_date=start_date, end_date=end_date)
res_S1.filter_duplicate_dates()
res_S1.dataframe

[4]:
title date swath_coverage_percentage relativeorbitnumber lastrelativeorbitnumber producttype sensoroperationalmode acquisitiontype polarisationmode beginposition platformname missiondatatakeid orbitdirection orbitnumber instrumentname lastorbitnumber endposition ingestiondate slicenumber platformidentifier
0 S1B_IW_GRDH_1SDV_20210101T172105_20210101T1721... 2021-01-01 100.0 88.0 88.0 GRD IW NOMINAL VV VH 2021-01-01T17:21:05.249000 Sentinel-1 194718.0 ASCENDING 24964.0 Synthetic Aperture Radar (C-band) 24964.0 2021-01-01T17:21:30.248000 2021-01-01T21:41:44.626000 10.0 2016-025A

5 - Search S2#

[5]:

start_date = "2021-01-01"
end_date = "2021-01-17"
# Retrieves all S2 images corresponding to the aoi, start date, and end date
res_S2 = client.s2_search(aoi=aoi, start_date=start_date, end_date=end_date)
res_S2.dataframe

[5]:
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
6 S2B_32TMK_20210101_0_L2A 2021-01-01 32TMK 0.00 sentinel-2b 065 S2B_MSIL2A_20210101T102329_N0214_R065_T32TMK_2... 2021-01-01T10:29:42Z 100.0 0.0 0.0 0.00 0.0 0.0 0.00 100.00 0.0 0.0
5 S2A_32TMK_20210103_0_L2A 2021-01-03 32TMK 0.00 sentinel-2a 022 S2A_MSIL2A_20210103T101411_N0214_R022_T32TMK_2... 2021-01-03T10:19:48Z 100.0 0.0 0.0 0.00 0.0 0.0 0.00 100.00 0.0 0.0
4 S2A_32TMK_20210106_0_L2A 2021-01-06 32TMK 0.00 sentinel-2a 065 S2A_MSIL2A_20210106T102411_N0214_R065_T32TMK_2... 2021-01-06T10:29:44Z 100.0 0.0 0.0 0.00 0.0 0.0 0.00 100.00 0.0 0.0
3 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.0 0.00 0.0 0.0 49.27 0.00 0.0 0.0
2 S2B_32TMK_20210111_0_L2A 2021-01-11 32TMK 5.08 sentinel-2b 065 S2B_MSIL2A_20210111T102309_N0214_R065_T32TMK_2... 2021-01-11T10:29:43Z 100.0 0.0 0.0 0.00 0.0 0.0 16.69 77.53 0.7 0.0
1 S2A_32TMK_20210113_0_L2A 2021-01-13 32TMK 0.00 sentinel-2a 022 S2A_MSIL2A_20210113T101401_N0214_R022_T32TMK_2... 2021-01-13T10:19:48Z 100.0 0.0 0.0 0.00 0.0 0.0 0.00 100.00 0.0 0.0
0 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.0 97.56 0.0 0.0 0.00 0.00 0.0 0.0

6 - Search Weather#

[6]:

start_date = "2021-01-01"
end_date = "2021-01-17"
# Retrieves temperature extrema corresponding to the aoi, start date, and end date
res_weather = client.weather_search(aoi=aoi, start_date=start_date, end_date=end_date, variables = ["AVGTEMP", "PREC"])
res_weather.dataframe

[6]:
date avgtemp prec
0 2021-01-01 9.388605 14.856473
1 2021-01-02 9.047998 12.129927
2 2021-01-03 7.030481 3.415852
3 2021-01-04 7.352747 9.735805
4 2021-01-05 7.531793 13.345673
5 2021-01-06 6.930902 0.578632
6 2021-01-07 6.775079 4.149154
7 2021-01-08 9.184137 9.411143
8 2021-01-09 12.055780 9.186649
9 2021-01-10 11.520410 10.954764
10 2021-01-11 8.974512 10.341069
11 2021-01-12 9.047998 0.744074
12 2021-01-13 12.357477 4.554079
13 2021-01-14 12.422113 1.272354
14 2021-01-15 11.482599 8.809325
15 2021-01-16 9.194818 2.330621
16 2021-01-17 11.172235 1.816961

7 - Filter search results#

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

[7]:

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

[7]:
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 S2A_32TMK_20210116_0_L2A 2021-01-16 32TMK 100.0 sentinel-2a 065 S2A_MSIL2A_20210116T102351_N0214_R065_T32TMK_2... 2021-01-16T10:29:44Z 100.0 0.0 0.0 97.56 0.0 0.0 0.0 0.0 0.0 0.0 
[8]:

res_weather.dataframe = res_weather.dataframe[[date in list(res_S2.dataframe["date"])+list(res_S1.dataframe["date"]) for date in list(res_weather.dataframe["date"])]]

8 - Select S2 bands#

[9]:

res_S2.output_bands = ["RED", "GREEN", "BLUE", "NIR"]

9 - Provide custom data#

The custom data used in this notebook can be found here. Please note, the reading in the Vector shapefile (“.shp”) needs all the files with “.cpg”, “.dbf”, “.prj”, “.qmd”, and “.shx” to be present in the same directory.

[10]:

# Example vector of land use over the AOI
# The "attributes" parameter can be used to select individual layers to fuse
land_use = Vector("../resources/example_files/Sardinia_land_use.shp", attributes=["IDFEATURE"])

# Example raster of NDVI (say for comparison) over the AOI
ndvi = Raster("../resources/example_files/ndvi_crop_sardinia.tif")

10 - Obtain all S1 and S2 images, and Fuse with custom files#

[11]:

fusion = client.fuse(
    catalogs_list=[res_S1,res_S2,res_weather],
    to_fuse=[land_use, ndvi],
    additional_info={"Info": "You can put any additional information, such as custom non-raster or non-vector data, in this parameter"}
    )

2022-11-28 09:53:21,611 INFO spacesense.core : created everything

11 - See the fused dataset#

The result of the Fusion is an Xarray Dataset. This is a 3 dimensional array with a structure and in-memory representation of a NetCDF file.

Here is a link to several tutorials from Xarray. And this link gives some introductory tutorials including basic data structure organization, selection, and computation.

For this example, the Dataset has dimensions of latitude, longitude, and time, with the Data Variables being the data layers of S1, S2, and the custom raster and vector that we fused together. As custom rasters and vectors do not currently support association with a date, they are fused as time-independant variables, thus only having an x (latitude) and y (longitude) dimension. S1 and S2 bands were not changed from their default values, which can be found here for S1, and here for S2.

[12]:

fusion.dataset

[12]:

<xarray.Dataset>
Dimensions:                      (time: 2, y: 43, x: 88)
Coordinates:
  * time                         (time) datetime64[ns] 2021-01-01 2021-01-16
  * y                            (y) float32 39.83 39.83 39.83 ... 39.83 39.83
  * x                            (x) float64 8.623 8.623 8.623 ... 8.63 8.63
Data variables:
    S1_VH                        (time, y, x) float32 0.0 0.0206 ... nan nan
    S1_VV                        (time, y, x) float32 0.0 0.1091 ... nan nan
    S1_DV                        (time, y, x) float32 nan 5.294 ... nan nan
    S1_LIA                       (time, y, x) float32 0.0 38.73 ... nan nan
    S2_RED                       (time, y, x) float32 nan nan nan ... 0.0417 0.0
    S2_GREEN                     (time, y, x) float32 nan nan nan ... 0.0371 0.0
    S2_BLUE                      (time, y, x) float32 nan nan nan ... 0.0255 0.0
    S2_NIR                       (time, y, x) float32 nan nan nan ... 0.2182 0.0
    WEATHER_AVGTEMP              (time) float32 9.389 9.195
    WEATHER_PREC                 (time) float32 14.86 2.331
    Sardinia_land_use_IDFEATURE  (y, x) float32 1.361e+05 ... 1.361e+05
    ndvi_crop_sardinia_BAND1     (y, x) float32 nan nan nan nan ... nan nan nan
Attributes:
    crs:                   +init=epsg:4326
    transform:             [ 8.97434772e-05  0.00000000e+00  8.62246354e+00  ...
    res:                   [8.97434772e-05 9.03483104e-05]
    ulx, uly:              [ 8.62246354 39.83104648]
    weather_data_lineage:  {"Data origin": "ERA5 reanalysis by ECMWF download...
    s1_data_lineage:       {"Data origin": "S3 bucket (ARN=arn:aws:s3:::senti...
    s2_data_lineage:       {"Data origin": "S3 bucket (ARN=arn:aws:s3:::senti...
    additional_info:       {'Info': 'You can put any additional information, ...

12 - Plot the fused dataset with matplotlib#

[13]:

ds_p = fusion.dataset.drop(["WEATHER_AVGTEMP", "WEATHER_PREC"])

# Configure matplotlib to show the 11 layers of fused dataset in a 6x2 grid
fig, axs = plt.subplots(nrows=5,ncols=2, figsize=(16, 30))

# Loop over the 11 variables in the fused dataset, select a date, and plot the layer
for idx, data_var in enumerate(ds_p.data_vars):
    # Set the variable to plot, and take the first date
    var = ds_p.isel(time=0)[data_var]
    # If the first date is all nans (no data for that variable), attempt to plot the next date
    for i, time in enumerate(ds_p.time.values):
        if np.isnan(var.mean().values) == True:
            var = ds_p.isel(time=i+1)[data_var]
    if "WEATHER" not in data_var:
        ax = fig.axes[idx]
        var.plot(ax=ax,cbar_kwargs={'label':""})
        ax.set_title(data_var)
        ax.axis('off')
        plt.draw()
../_images/notebooks_5_fuse_raster_vector_29_0.png 
[14]:

# Select and plot individual layer and time
fusion.dataset.S2_RED.sel(time="2021-01-16").plot()

[14]:

<matplotlib.collections.QuadMesh at 0x7f806c194590>
../_images/notebooks_5_fuse_raster_vector_30_1.png 
[16]:

rgb = fusion.plot_rgb(all_dates=False)
../_images/notebooks_5_fuse_raster_vector_31_0.png