Open SAR Toolkit - Tutorial 3, version 1.2, July 2020. Andreas Vollrath, ESA/ESRIN phi-lab

title

OST Tutorial III¶

Process the latest Sentinel-1 GRD product for a given point¶

Short description

This notebook demonstrates the interaction between the Sentinel1 class for data inventory and download, and the Sentinel1Scene class, together, for the generation of the latest Sentinel-1 product over a given point coordinate.

Requirements

a PC/Mac with at least 16GB of RAM
about 4GB of free disk space
a Copernicus Open Data Hub user account, ideally valid for at least 7 days (https://scihub.copernicus.eu)

NOTE: all cells that have an * after its number can be executed without changing any code.

0* - Install OST and dependencies¶

NOTE: Applies only if you haven’t fully installed OST and its dependencies yet, e.g. on Google Colab, so uncomment the lines in this case.

[ ]:

# !apt-get -y install wget
# !wget https://raw.githubusercontent.com/ESA-PhiLab/OST_Notebooks/master/install_ost.sh
# !bash install_ost.sh

1* - Import python libraries necessary for processing¶

[ ]:

# this is for the path handling and especially useful if you are on Windows
from pathlib import Path
from pprint import pprint

# we will need this for our time of period definition
from datetime import datetime, timedelta

# this is the s1Project class, that basically handles all the workflow from beginning to the end
from ost import Sentinel1, Sentinel1Scene

2 - Data selection parameters¶

NOTE: In case you want to process a different area, all you want to change is the lat and lon values in line 6

As already covered in OST Tutorial 2, we need a minimum of 3 basic parameters to initialise the Sentinel1 class.

1 Area of Interest:

In this case we only search for a specific spot on earth, i.e. Rome, Italy, that is defined by the Latitude and Longitude. We then create a Well-Known-Text formatted string.

2 Time of Interest:

In this example, the datetime class is used to set the start date to 30 days before today to assure we get any scene within our time of interest.

3 Project directory

Set this to anything you like if not happy with the default one.

[ ]:

# ----------------------------
# Area of interest
# ----------------------------

# Here we can either point to a shapefile or as well use
lat, lon = 41.8919, 12.5113
aoi = "POINT ({} {})".format(lon, lat)

# ----------------------------
# Time of interest
# ----------------------------
# we set only the start date to today - 30 days
start = (datetime.today() - timedelta(days=30)).strftime("%Y-%m-%d")

# ----------------------------
# Processing directory
# ----------------------------
# get home folder
home = Path.home()
# create a processing directory within the home folder
project_dir = home / "OST_Tutorials" / "Tutorial_3"

# ------------------------------
# Print out AOI and start date
# ------------------------------
print(
    "AOI: " + aoi,
)
print("TOI start: " + start)
print("Our project directory is located at: " + str(project_dir))

3* - Initialize the Sentinel1 project class¶

After initialisation of our class, where we explicitley add the GRD product type argument, we do a rough search over our AOI for the last 30 days. We print the first 5 entries and plot all images for visualization by using the Sentinel1 class attribute inventory and method plot_inventory.

[ ]:

# ---------------------------------------------------
# for plotting purposes we use this iPython magic
%matplotlib inline
%pylab inline
pylab.rcParams["figure.figsize"] = (19, 19)
# ---------------------------------------------------

# create s1Project class instance
s1_project = Sentinel1(project_dir=project_dir, aoi=aoi, start=start, product_type="GRD")

# search command
s1_project.search()
# uncomment in case you have issues with the registration procedure
# ost_s1.search(base_url='https://scihub.copernicus.eu/dhus')
print("We found {} products.".format(len(s1_project.inventory.identifier.unique())))
# combine OST class attribute with pandas head command to print out the first 5 rows of the
print(s1_project.inventory.head(5))

# we plot the full Inventory on a map
s1_project.plot_inventory(transparency=0.1)

4* - Select the latest scene found during the search¶

Here we use some python-panda syntax on our rough data inventory to filter out the latest scene and create store it in a new dataframe.

[ ]:

pylab.rcParams["figure.figsize"] = (13, 13)

# we give our inventory a shorter name iDf (for inventory Dataframe)
iDf = s1_project.inventory.copy()

# we select the latest scene based on the metadata entry endposition
latest_df = iDf[iDf.endposition == iDf.endposition.max()]

# we print out a little info on the date of the
print(" INFO: Latest scene found for {}".format(latest_df["acquisitiondate"].values[0]))

# we use the plotInventory method in combination with the newly
# created Geodataframe to see our scene footprint
s1_project.plot_inventory(latest_df, transparency=0.5)

7* Download scene¶

We use the build-in download method from the Sentinel1 class. Note that you can pass any Geodataframe generated by OST, and filtered by you (e.g. sort out rows that you do not need). In our case we are only interested in the latest scene, so we pass the newly generated latest_df Geodataframe object.

NOTE that you should use ESA’s scihub server in this case, since it is the place where the images arrive first. Other data mirrors might have slight delays, so that the scene found by the inventory might not be available.

[ ]:

s1_project.download(latest_df)

8* - Display some metadata of the latest scene¶

After use of the Sentinel1 class for finding and downloading the latest scene, we hand the scene identifier over to the Sentinel1Scene class for further processing as already demonstrated in OST Tutorial 1.

[ ]:

# create a S1Scene class instance based on the scene identifier coming from our "latest scene dataframe"
latest_scene = Sentinel1Scene(latest_df["identifier"].values[0])

# print summarising infos
latest_scene.info()

# print file location
file_location = latest_scene.get_path(s1_project.download_dir)

print(" File is located at: ")
print(" " + str(file_location))

9* - Produce a subsetted ARD product¶

The creation of the ARD product follows the same logic as presented in OST Tutorial 1. However, for this case we introduce the subset argument to eh create_ard function. Subsetting is adviced if only a small portion of the whole image is of interest. It will speed up processing and uses less storage.

In our case we use some helper functions within the OST package to create a squared buffer area fo 10000 meter around our point of interest defined as AOI.

[ ]:

# 10 km buffer around AOI Point
from shapely.wkt import loads
from ost.helpers import vector as vec

# turn WKT into shapely geometry object
shp_aoi = loads(s1_project.aoi)

# use OST helper's function to create a quadrant buffer for subset
subset_area = vec.geodesic_point_buffer(shp_aoi.x, shp_aoi.y, 10000, envelope=True)

print("-----------------------------------------------------------------------------")
latest_scene.create_ard(
    # we see our download path can be automatically generated by providing the Project's download directory
    infile=latest_scene.get_path(download_dir=s1_project.download_dir),
    # let's simply take our processing folder
    out_dir=s1_project.processing_dir,
    # define the subset
    subset=subset_area,
    # in case already processed, we will re-process
    overwrite=True,
)

print("-----------------------------------------------------------------------------")
print(" The path to our newly created ARD product can be obtained the following way:")
latest_scene.ard_dimap

10* - Create a RGB color composite¶

As already demonstrated in OST Tutorial 1, we create an RGB GeoTiff, and visualize it.

[ ]:

latest_scene.create_rgb(outfile=s1_project.processing_dir / f"{latest_scene.start_date}.tif")
latest_scene.visualise_rgb(shrink_factor=1)