Open SAR Toolkit - Tutorial 4a, version 1.2, June 2020. Andreas Vollrath, ESA/ESRIN phi-lab

title


OST Tutorial IV-A

How to create near-daily timeseries over Vienna. Introduction to GRD Batch Processing part I.

Open In Colab


Short description

This notebook provides an introduction to the batch processing of Sentinel-1 GRD data using OST’s Sentinel1Batch class. This is a subclass of the Sentinel1 class, and thus inherits all the functionalities of the Generic and the Sentinel1 classes for the generation of a project as well as data search and refinement as presented in the OST Tutorial II notebook. The Sentinel1Batch class holds functions for the batch processing of single calibrated backscatter products. Furthermore, time-series processing and the generation of multi-temporal statistics, referred to as timescans, are introduced.

Within the given example, time-series for 4 different overlapping tracks are going to be produced over the city of Vienna, Austria. The notebook demonstrates:

  1. the reduction of data processing by automatically subsetting the data,

  2. time-series processing and the corresponding ARD types,

  3. merging the track specific time-series into a unique time-series with almost daily coverage,

  4. creation of a timeseries animation for outreach purposes.


Requirements

  • a PC/Mac with at least 16GB of RAM

  • about 25 GB 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 yet, e.g. on Google Colab,

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

1* - Import of Libraries

In this step we import some standard python libraries for OS independent path handling as well as the Sentinel1_GRDBatch class thta handles the full workflow from search, download and processing of multiple GRD scenes. In addition, the OST helper module vector is loaded to create an AOI based on Point coordinates, and the raster module for creating a time-series animation.

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


# this is the s1Project class, that basically handles all the workflow from beginning to the end
from ost import Sentinel1Batch
from ost.helpers import vector, raster

2* - Set up the project

Here you going to initialize the Sentinel1Batch class by determining the project folder, the AOI and the start and end date. In addition we determine the image product type (i.e. GRD) and the ARD type that we use for processing. In this cas we choose the OST-RTC that produces Radiometrically Terrain Corrected products, i.e. the images are corrected for radiometric distortions along mountainous slopes. This type of ARD is advised when doing land cover and land use studiesover rugged terrain.

[ ]:
# define the project directory
project_dir = Path.home() / "OST_Tutorials" / "Tutorial_4a"

# define aoi with a 2 point coordinates and create a buffer with 20km radius
lat, lon = "48.25", "16.4"  # Vienna
aoi = vector.latlon_to_wkt(lat, lon, buffer_meter=17500, envelope=True)

# define the start and end date
start = "2020-05-01"
end = "2020-05-31"

# initialize the class to s1_grd instance
s1_grd = Sentinel1Batch(
    project_dir=project_dir,
    aoi=aoi,
    start=start,
    end=end,
    product_type="GRD",
    ard_type="OST-RTC",
)

# do the search
if not s1_grd.inventory_file:
    s1_grd.search()

3* - Plot refined data inventory

The resultant dataframe from the search inventory is visualised. We do not do a refinement step here, since all images are fully overlapping the AOI. This allows us to create a combined, almost daily time-series of all images.

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

# plot the inventory
s1_grd.plot_inventory(s1_grd.inventory, transparency=0.1)

4* - Download of GRD scenes

As already shown in Tutorial II, you will download the scenes based on the refined inventory dataframe for the respective produckt key.

[ ]:
s1_grd.download(s1_grd.inventory, concurrent=10)

5* - Set ARD parameters

Similar to the Sentinel1Scene class (Tutorial I and III), the Sentinel1Batch class handles the defintion of ARD types in a hierarchical dictionary structure. You can use the same types and steps to customize as for the Sentinel1Scene class. For our example, we already intialised the class instance with the OST-RTC ARD type, in order to remove the backscatter distortion on mountainous slopes. This applies to all single image processing in the first step. The subsequent time-series processing will bring all the imagery to a common grid and apply a multitemporal speckle filter, that is much more efficient than speckle filtering applied on a single image. Since speckle filters a conceptionalised to work on the raw power data of SAR imagery, the conversion to the dB scale is only applied after the multi-temporal speckle filtering. In addition, all images are clipped to the very same extent, that is defined by the common data coverage of all images per track as well as the AOI.

Note that it is possible to change the datatype of the output to either unsigned 16 or 8-bit integer. The backscatter data is therefore linearly stretched between -30 to 5 dB. This has the advantage to reduce the necessary storage by a factor of 2 for 16-bit uint and a factor of 4 for 8-bit uint data.

The exact processing steps are as follows and depend on the ARD type:

[ ]:
print(
    "-----------------------------------------------------------------------------------------------------------"
)
print("Time-series processing parameters hold in the configuration file:")
print(
    "-----------------------------------------------------------------------------------------------------------"
)
print("")
pprint(s1_grd.ard_parameters["time-series_ARD"])

# custimze some single scene ARD parameters
s1_grd.ard_parameters["single_ARD"]["resolution"] = 50  # reduce for processing time and disk space
s1_grd.ard_parameters["time-series_ARD"]["dtype_output"] = "uint8"

6* - Run the batch routine

To process all the data, including time-series and timescans is as easy as one command. All the complexity is handled by OST in the back, and you just have to wait, since processing can take quite a while. Note the keywords to aly higher level tim-series and timescan generation. Mosaicking refers toacross track mosaicking, which for this example is not the case. The overwrite argument tells OST if it should start processing from scratch (i.e. set to True), or process from where it stopped. The latter comes in handy when wokng on cloud instances that crash or automatically shutdown every once in a while.

Note that subsetting is set automatically to True when all tracks hold in the inventory fully overlap the AOI.

[ ]:
s1_grd.grds_to_ards(
    inventory_df=s1_grd.inventory,
    timeseries=True,
    timescan=True,
    mosaic=False,
    overwrite=False,
)

7* - Merge single per-track time-series to one single time-series

By using a little helper function from OST’s raster module, combining the 4 time-series to a unique one is as easy the following command. Within your processing directory, a new foler called combined is created. If multi-temporal statistics for this new time-series should be created, set the timescan argument to True.

[ ]:
raster.combine_timeseries(processing_dir=s1_grd.processing_dir, config_dict=s1_grd.config_dict, timescan=True)

8* - Create a time-series animation

Finally, a time-series animation is created. therefore we need to pass the time-series folder to the command. product_list expects a list of 1 to 3 elements. For GRD data this is either a single polarisation, or both bands. OST will calculate the power ratio of band 1 and 2 for a 3-band RGB composite. A shrink factor can be set to reduce image resolution and memory needs.

Note: This needs imagemagick installed, which is not a default requirement by OST. You can install it on e.g. Ubuntu by typing: sudo apt-get install magick

[ ]:
# define Time-series folder
ts_folder = s1_grd.processing_dir / "combined" / "Timeseries"

# create the animation
raster.create_timeseries_animation(
    timeseries_folder=ts_folder,
    product_list=["bs.VV", "bs.VH"],
    out_folder=s1_grd.processing_dir,
    shrink_factor=3,
    add_dates=True,
)