Automatic ice damage detection from Sentinel-1 radar imagery¶

Step 2: Get Sentinel-1 imagery data.¶

For more info about this science case, please check the documentation at this link.

In this notebook SentinelHub is exploited to query, preprocess and download Sentinel-1 data.
The imagery data is then saved into a datacube, with the same resolution and extent as the damage model.

NB: A SentinelHub account is required in order to use its functionalities. By default the sh.SHConfig() instance will look for the SH_CLIENT_ID and SH_CLIENT_SECRET environment variables.
If they are not set, they can be passed to the sh.SHConfig() instance after its creation by running sh_config.sh_client_id = <ID> and sh_config.sh_client_secret = <SECRET>.

This notebook runs with the python environment users-science-case-polar, please checkout the documentation for help on changing the environment.

In [1]:

# import the necessary packages

import os
from datetime import datetime, timedelta
from pathlib import Path

import sentinelhub as sh
import xarray as xr
from xcube.core.store import new_data_store

# import the necessary packages import os from datetime import datetime, timedelta from pathlib import Path import sentinelhub as sh import xarray as xr from xcube.core.store import new_data_store

Define the parameters needed to access the S3 storage. They are saved as environment variables.

In [2]:

S3_USER_STORAGE_KEY = os.environ["S3_USER_STORAGE_KEY"]
S3_USER_STORAGE_SECRET = os.environ["S3_USER_STORAGE_SECRET"]
S3_USER_STORAGE_BUCKET = os.environ["S3_USER_STORAGE_BUCKET"]

team_store = new_data_store(
    data_store_id="s3",
    root=S3_USER_STORAGE_BUCKET,
    storage_options=dict(
        anon=False, key=S3_USER_STORAGE_KEY,
        secret=S3_USER_STORAGE_SECRET
    )
)

S3_USER_STORAGE_KEY = os.environ["S3_USER_STORAGE_KEY"] S3_USER_STORAGE_SECRET = os.environ["S3_USER_STORAGE_SECRET"] S3_USER_STORAGE_BUCKET = os.environ["S3_USER_STORAGE_BUCKET"] team_store = new_data_store( data_store_id="s3", root=S3_USER_STORAGE_BUCKET, storage_options=dict( anon=False, key=S3_USER_STORAGE_KEY, secret=S3_USER_STORAGE_SECRET ) )

Set the directory where the Sentinel-1 scenes will be downloaded.

In [3]:

iceshelf_name = "amery"

download_directory = Path.cwd() / "sh_download_directory" / iceshelf_name
download_directory.mkdir(parents=True, exist_ok=True)

iceshelf_name = "amery" download_directory = Path.cwd() / "sh_download_directory" / iceshelf_name download_directory.mkdir(parents=True, exist_ok=True)

Query setup¶

Coordinate Reference Systems

• EPSG 3031 is Antarctic Polar Stereographic projected coordinates
• UTM depends on the longitude of the area of interest and is needed to specify the resolution / size of the imagery data, due to how SentinelHub is implemented

Bounding box of the area of interest

In [4]:

iceshelf_bbox = sh.BBox((1566250, 533250, 2266750, 933750), crs=sh.CRS(3031))

# resolution in [m], in EPSG:3031 projected coordinates
resolution_x, resolution_y = (500, 500)

# size of the imagery array given the resolution specified above
imagery_width = round(abs(iceshelf_bbox.max_x - iceshelf_bbox.min_x) / resolution_x)
imagery_height = round(abs(iceshelf_bbox.max_y - iceshelf_bbox.min_y) / resolution_y)

iceshelf_bbox = sh.BBox((1566250, 533250, 2266750, 933750), crs=sh.CRS(3031)) # resolution in [m], in EPSG:3031 projected coordinates resolution_x, resolution_y = (500, 500) # size of the imagery array given the resolution specified above imagery_width = round(abs(iceshelf_bbox.max_x - iceshelf_bbox.min_x) / resolution_x) imagery_height = round(abs(iceshelf_bbox.max_y - iceshelf_bbox.min_y) / resolution_y)

Time interval of interest

As an example, let's download the data with acquisition date within the first 5 days of October 2016.
This use case employs data from the last three months (Q4) of both 2015 and 2016.

In [5]:

time_start = datetime(2016, 10, 1)
time_end = datetime(2016, 10, 5)

# set the time interval in the format required by the SentinelHub API
time_interval = time_start.strftime("%Y-%m-%d"), time_end.strftime("%Y-%m-%d")

# this is useful to distinguish the files in case data is downloaded for multiple time intervals
time_period_label = "test_Oct2016"

time_start = datetime(2016, 10, 1) time_end = datetime(2016, 10, 5) # set the time interval in the format required by the SentinelHub API time_interval = time_start.strftime("%Y-%m-%d"), time_end.strftime("%Y-%m-%d") # this is useful to distinguish the files in case data is downloaded for multiple time intervals time_period_label = "test_Oct2016"

Sentinel-1 categories

• Acquisition mode
  • IW: Interferometric Wide swath mode
  • EW: Extra Wide swath mode
• Polarisation
  • SH: only HH co-polarisation
  • SV: only VV co-polarisation
  • DH: double polarisation, both co- (HH) and cross- (HV)
  • DV: double polarisation, both co- (VV) and cross- (VH)
• Orbit
  • Ascending
  • Descending

Evalscripts custom scripts

Sentinel-1 queries can be performed with custom Javascript scripts (also called evalscripts).
Check out the SentinelHub documentation at this link for a thorough explanation of what the evalscripts look like and what they can be used for.

What follows are the v3 evalscripts used to query Sentinel-1 data over the polar regions, where only DH and SH polarisation modes are available.
Only the latter is used in this notebook, to get all the scenes in single polarisation mode containing only co-polarised data.

In [6]:

# this is for the scenes with double polarisation

evalscript_s1_pol_dh = """
    //VERSION=3

    function setup() {
        return {
            input: [{
                bands: ["HH", "HV"]
            }],
            output: {
                bands: 2, sampleType: "FLOAT32"
            }
        };
    }

    function evaluatePixel(sample) {
        return [sample.HH, sample.HV];
    }
"""


# this is for the scenes with single polarisation

evalscript_s1_pol_sh = """
    //VERSION=3

    function setup() {
        return {
            input: [{
                bands: ["HH"]
            }],
            output: {
                bands: 1, sampleType: "FLOAT32"
            }
        };
    }

    function evaluatePixel(sample) {
        return [sample.HH];
    }
"""

# this is for the scenes with double polarisation evalscript_s1_pol_dh = """ //VERSION=3 function setup() { return { input: [{ bands: ["HH", "HV"] }], output: { bands: 2, sampleType: "FLOAT32" } }; } function evaluatePixel(sample) { return [sample.HH, sample.HV]; } """ # this is for the scenes with single polarisation evalscript_s1_pol_sh = """ //VERSION=3 function setup() { return { input: [{ bands: ["HH"] }], output: { bands: 1, sampleType: "FLOAT32" } }; } function evaluatePixel(sample) { return [sample.HH]; } """

Query and download the data¶

Find the scenes that overlap the area of interest.¶

In [7]:

# set up the Sentinel Hub Config object
sh_config = sh.SHConfig()
sh_config.sh_base_url = sh.DataCollection.SENTINEL1.service_url
sh_catalog = sh.SentinelHubCatalog(config=sh_config)

# the search is performed only based on the bounding box and the time interval
search_iterator = sh_catalog.search(
    sh.DataCollection.SENTINEL1,
    bbox=iceshelf_bbox,
    time=time_interval,
    fields={"include": ["id", "properties.datetime"],
            "exclude": []},
)
results = list(search_iterator)
print("Total number of results:", len(results))

# set up the Sentinel Hub Config object sh_config = sh.SHConfig() sh_config.sh_base_url = sh.DataCollection.SENTINEL1.service_url sh_catalog = sh.SentinelHubCatalog(config=sh_config) # the search is performed only based on the bounding box and the time interval search_iterator = sh_catalog.search( sh.DataCollection.SENTINEL1, bbox=iceshelf_bbox, time=time_interval, fields={"include": ["id", "properties.datetime"], "exclude": []}, ) results = list(search_iterator) print("Total number of results:", len(results))

Total number of results: 8

Have a look at the IDs of the scenes.

In [8]:

search_iterator.get_ids()

search_iterator.get_ids()

Out[8]:

['S1B_EW_GRDM_1SSH_20161004T152819_20161004T152923_002358_003FBD_384F',
 'S1A_EW_GRDM_1SSH_20161004T143928_20161004T144032_013341_015466_DA25',
 'S1A_IW_GRDH_1SSH_20161003T220432_20161003T220502_013331_015416_78EA',
 'S1A_IW_GRDH_1SSH_20161003T220407_20161003T220432_013331_015416_12CA',
 'S1A_IW_GRDH_1SSH_20161003T220342_20161003T220407_013331_015416_1A83',
 'S1A_IW_GRDH_1SSH_20161003T220317_20161003T220342_013331_015416_4FF6',
 'S1A_IW_GRDH_1SSH_20161003T220252_20161003T220317_013331_015416_47F5',
 'S1B_EW_GRDM_1SSH_20161003T144654_20161003T144758_002343_003F5C_63D4']

And their acquisition timestamps.

In [9]:

search_iterator.get_timestamps()

search_iterator.get_timestamps()

Out[9]:

[datetime.datetime(2016, 10, 4, 15, 28, 19, tzinfo=tzlocal()),
 datetime.datetime(2016, 10, 4, 14, 39, 28, tzinfo=tzlocal()),
 datetime.datetime(2016, 10, 3, 22, 4, 32, tzinfo=tzlocal()),
 datetime.datetime(2016, 10, 3, 22, 4, 7, tzinfo=tzlocal()),
 datetime.datetime(2016, 10, 3, 22, 3, 42, tzinfo=tzlocal()),
 datetime.datetime(2016, 10, 3, 22, 3, 17, tzinfo=tzlocal()),
 datetime.datetime(2016, 10, 3, 22, 2, 52, tzinfo=tzlocal()),
 datetime.datetime(2016, 10, 3, 14, 46, 54, tzinfo=tzlocal())]

Preprocess and download the scenes and save them as datacubes.¶

In [10]:

# define a (narrow) time window (in this case of 1 minute) within which subsequent scenes can be joined into the same mosaic
max_time_difference_for_mosaicing = timedelta(seconds=60)

for scene_id, scene_timestamp in zip(search_iterator.get_ids(), search_iterator.get_timestamps()):

    # get the polarisation mode and orbit type from the scene metadata
    polarisation = sh_catalog.get_feature(sh.DataCollection.SENTINEL1, scene_id)["properties"]["s1:polarization"]
    orbit = sh_catalog.get_feature(sh.DataCollection.SENTINEL1, scene_id)["properties"]["sat:orbit_state"]

    # set the directory where to download the data and the name of the zarr file that will be created
    scene_download_directory = download_directory / orbit / f"{polarisation.lower()}_data" / time_period_label / f"{scene_timestamp.strftime('%Y%m%dT%H%M%S')}"
    zarr_output_file_path = download_directory / orbit / f"{polarisation.lower()}_data" / time_period_label / f"{scene_id}.zarr"

    # if a scene is DH don't make any request
    if polarisation == "DH":
        continue

    # now we can expect that all the scenes will be SH
    assert polarisation == "SH"
    request = sh.SentinelHubRequest(
        evalscript=evalscript_s1_pol_sh,
        input_data=[
            sh.SentinelHubRequest.input_data(
                data_collection=sh.DataCollection.SENTINEL1,
                time_interval=(
                    scene_timestamp - max_time_difference_for_mosaicing,
                    scene_timestamp + max_time_difference_for_mosaicing),
                other_args={"processing": {
                    "orthorectify": True,
                    "demInstance": "COPERNICUS",
                    "backCoeff": "GAMMA0_TERRAIN"}}
            )
        ],
        responses=[sh.SentinelHubRequest.output_response("default", sh.MimeType.TIFF)],
        bbox=iceshelf_bbox,
        size=(imagery_width, imagery_height),
        config=sh_config,
        data_folder=scene_download_directory.as_posix(),
    )

    # download the scene and save into a datacube
    downloaded_data = sh.SentinelHubDownloadClient(config=sh_config).download(request.download_list, max_threads=1)
    tiff_path = scene_download_directory / request.get_filename_list()[0]
    tiff_file = xr.open_mfdataset(tiff_path, engine="rasterio")
    imagery_datacube = tiff_file.isel(band=0).drop_vars("band").expand_dims(
        dim={"time": [int(scene_timestamp.timestamp())]},
        axis=0, create_index_for_new_dim=True)
    imagery_datacube = imagery_datacube.reindex(y=imagery_datacube.y[::-1])
    imagery_datacube = imagery_datacube.rename({"band_data": "s1_imagery"})
    imagery_datacube = imagery_datacube.chunk(dict(time=1, x=128, y=128))
    imagery_datacube.to_zarr(zarr_output_file_path)

# define a (narrow) time window (in this case of 1 minute) within which subsequent scenes can be joined into the same mosaic max_time_difference_for_mosaicing = timedelta(seconds=60) for scene_id, scene_timestamp in zip(search_iterator.get_ids(), search_iterator.get_timestamps()): # get the polarisation mode and orbit type from the scene metadata polarisation = sh_catalog.get_feature(sh.DataCollection.SENTINEL1, scene_id)["properties"]["s1:polarization"] orbit = sh_catalog.get_feature(sh.DataCollection.SENTINEL1, scene_id)["properties"]["sat:orbit_state"] # set the directory where to download the data and the name of the zarr file that will be created scene_download_directory = download_directory / orbit / f"{polarisation.lower()}_data" / time_period_label / f"{scene_timestamp.strftime('%Y%m%dT%H%M%S')}" zarr_output_file_path = download_directory / orbit / f"{polarisation.lower()}_data" / time_period_label / f"{scene_id}.zarr" # if a scene is DH don't make any request if polarisation == "DH": continue # now we can expect that all the scenes will be SH assert polarisation == "SH" request = sh.SentinelHubRequest( evalscript=evalscript_s1_pol_sh, input_data=[ sh.SentinelHubRequest.input_data( data_collection=sh.DataCollection.SENTINEL1, time_interval=( scene_timestamp - max_time_difference_for_mosaicing, scene_timestamp + max_time_difference_for_mosaicing), other_args={"processing": { "orthorectify": True, "demInstance": "COPERNICUS", "backCoeff": "GAMMA0_TERRAIN"}} ) ], responses=[sh.SentinelHubRequest.output_response("default", sh.MimeType.TIFF)], bbox=iceshelf_bbox, size=(imagery_width, imagery_height), config=sh_config, data_folder=scene_download_directory.as_posix(), ) # download the scene and save into a datacube downloaded_data = sh.SentinelHubDownloadClient(config=sh_config).download(request.download_list, max_threads=1) tiff_path = scene_download_directory / request.get_filename_list()[0] tiff_file = xr.open_mfdataset(tiff_path, engine="rasterio") imagery_datacube = tiff_file.isel(band=0).drop_vars("band").expand_dims( dim={"time": [int(scene_timestamp.timestamp())]}, axis=0, create_index_for_new_dim=True) imagery_datacube = imagery_datacube.reindex(y=imagery_datacube.y[::-1]) imagery_datacube = imagery_datacube.rename({"band_data": "s1_imagery"}) imagery_datacube = imagery_datacube.chunk(dict(time=1, x=128, y=128)) imagery_datacube.to_zarr(zarr_output_file_path)

Combine the data cubes and save them to S3 storage.¶

At this point we should have one folder containing the scenes with ascending orbit and one with descending orbit.
All scenes have "SH" polarisation mode, hence containing only "HH" co-polarised data.

In [11]:

polarisation_mode = "SH"
polarisation_type = "HH"

for orbit_type in ["ascending", "descending"]:

    data_folder = download_directory / orbit_type / f"{polarisation_mode.lower()}_data" / time_period_label

    scenes_list = list(data_folder.glob("*.zarr"))
    scenes_list.sort(key=os.path.getmtime, reverse=True)
    imagery_datacube = xr.open_mfdataset(
        scenes_list, drop_variables=["spatial_ref"])

    del imagery_datacube.s1_imagery.attrs["AREA_OR_POINT"]
    del imagery_datacube.s1_imagery.attrs["TIFFTAG_RESOLUTIONUNIT"]
    del imagery_datacube.s1_imagery.attrs["TIFFTAG_XRESOLUTION"]
    del imagery_datacube.s1_imagery.attrs["TIFFTAG_YRESOLUTION"]
    del imagery_datacube.s1_imagery.attrs["grid_mapping"]
    imagery_datacube["s1_imagery"].attrs = {
        "description": "Imagery scenes",
        "units": "Normalised backscatter (linear power)."}
    imagery_datacube = imagery_datacube.rio.write_crs("epsg:3031", grid_mapping_name="crs")
    imagery_datacube = imagery_datacube.assign_attrs(
        description=f"Sentinel-1 imagery over Amery Ice Shelf for the period '{time_period_label}'. \
    Polarisation mode: {polarisation_mode}. Polarisation: {polarisation_type}. Orbit: {orbit_type}.")

    storage_path = f"datacubes/{iceshelf_name}/s1_imagery_res_500m_{polarisation_mode}_{polarisation_type}_{orbit_type.capitalize()}_{time_period_label}.zarr"
    team_store.write_data(
        imagery_datacube, storage_path, replace=False)

polarisation_mode = "SH" polarisation_type = "HH" for orbit_type in ["ascending", "descending"]: data_folder = download_directory / orbit_type / f"{polarisation_mode.lower()}_data" / time_period_label scenes_list = list(data_folder.glob("*.zarr")) scenes_list.sort(key=os.path.getmtime, reverse=True) imagery_datacube = xr.open_mfdataset( scenes_list, drop_variables=["spatial_ref"]) del imagery_datacube.s1_imagery.attrs["AREA_OR_POINT"] del imagery_datacube.s1_imagery.attrs["TIFFTAG_RESOLUTIONUNIT"] del imagery_datacube.s1_imagery.attrs["TIFFTAG_XRESOLUTION"] del imagery_datacube.s1_imagery.attrs["TIFFTAG_YRESOLUTION"] del imagery_datacube.s1_imagery.attrs["grid_mapping"] imagery_datacube["s1_imagery"].attrs = { "description": "Imagery scenes", "units": "Normalised backscatter (linear power)."} imagery_datacube = imagery_datacube.rio.write_crs("epsg:3031", grid_mapping_name="crs") imagery_datacube = imagery_datacube.assign_attrs( description=f"Sentinel-1 imagery over Amery Ice Shelf for the period '{time_period_label}'. \ Polarisation mode: {polarisation_mode}. Polarisation: {polarisation_type}. Orbit: {orbit_type}.") storage_path = f"datacubes/{iceshelf_name}/s1_imagery_res_500m_{polarisation_mode}_{polarisation_type}_{orbit_type.capitalize()}_{time_period_label}.zarr" team_store.write_data( imagery_datacube, storage_path, replace=False)

In [ ]: