Search and fuse S1 and S2 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 and Sentinel-2 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#

[2]:
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#

[1]:
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

[5]:
# 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_s1_s2_with_custom_data")

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

4 - Search S1#

[6]:
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.dataframe
[6]:
title date relativeorbitnumber lastrelativeorbitnumber producttype sensoroperationalmode acquisitiontype polarisationmode beginposition platformname missiondatatakeid orbitdirection orbitnumber instrumentname lastorbitnumber endposition ingestiondate slicenumber platformidentifier
0 S1A_IW_GRDH_1SDV_20210101T052857_20210101T0529... 2021-01-01 168.0 168.0 GRD IW NOMINAL VV VH 2021-01-01T05:28:57.811000 Sentinel-1 275881.0 DESCENDING 35940.0 Synthetic Aperture Radar (C-band) 35940.0 2021-01-01T05:29:22.811000 2021-01-01T09:42:44.752000 24.0 2014-016A

5 - Search S2#

[16]:
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
[16]:
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 - Filter search results#

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

[17]:
res_S2.dataframe = res_S2.dataframe[res_S2.dataframe["valid_pixel_percentage"] > 80]
res_S2.dataframe
[17]:
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

7 - Provide custom data#

[18]:
# 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")

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

[19]:
fusion = client.fuse(
    catalogs_list=[res_S1,res_S2],
    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-09-14 14:48:13,930 INFO spacesense.core : created everything

9 - 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.

[20]:
fusion.dataset
[20]:
<xarray.Dataset>
Dimensions:                           (time: 2, y: 43, x: 88)
Coordinates:
  * time                              (time) datetime64[ns] 2021-01-01 2021-0...
  * y                                 (y) float32 39.83 39.83 ... 39.83 39.83
  * x                                 (x) float64 8.623 8.623 ... 8.63 8.63
Data variables:
    s1_vh                             (time, y, x) float32 0.0 0.02966 ... nan
    s1_vv                             (time, y, x) float32 0.0 0.07535 ... nan
    s1_lia                            (time, y, x) float32 0.0 39.92 ... nan nan
    s1_mask                           (time, y, x) float32 0.0 255.0 ... nan nan
    s2_B02                            (time, y, x) float32 nan nan ... 0.001136
    s2_B03                            (time, y, x) float32 nan nan ... 0.001654
    s2_B04                            (time, y, x) float32 nan nan ... 0.001859
    s2_B08                            (time, y, x) float32 nan nan ... 0.009711
    s2_SCL                            (time, y, x) float32 nan nan ... 0.1781
    Sardinia_land_use_band_IDFEATURE  (y, x) float32 1.285e+05 ... 1.361e+05
    ndvi_crop_sardinia_band_1         (y, x) float32 nan nan nan ... nan nan nan
Attributes:
    crs:              +init=epsg:4326
    transform:        [ 8.97434772e-05  0.00000000e+00  8.62246354e+00  0.000...
    res:              [8.97434772e-05 9.03483104e-05]
    ulx, uly:         [ 8.62246354 39.83104648]
    additional_info:  {'Info': 'You can put any additional information, such ...

10 - Plot the fused dataset with matplotlib#

[21]:
ds_p = fusion.dataset

# Configure matplotlib to show the 11 layers of fused dataset in a 6x2 grid
fig, axs = plt.subplots(nrows=6,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]
    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_24_0.png
[22]:
# Select and plot individual layer and time
fusion.dataset.s2_B04.sel(time="2021-01-16").plot()
[22]:
<matplotlib.collections.QuadMesh at 0x7f6f7c3bd8d0>
../_images/notebooks_5_fuse_raster_vector_25_1.png