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

---

OST Tutorial IV-B

How to create a timeseries animation of Iceberg A-68. Introduction to GRD Batch Processing part II.

---

Short description

This notebook continues to introduce you to the general workflow of OST for the batch processing of GRD data using the Sentinel1Batch class. In this example:

1. across-track mosaicking based on a refined inventory and

2. processing in polar regions is shown.

---

Requirements

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

• about 75 GB of free disk space

• a Copernicus Open Data Hub user account, 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

from pathlib import Path
from pprint import pprint

from ost import Sentinel1Batch
from ost.helpers import vector, raster

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

2 - Set up the project

This follows the logic of the prior OST Tutorial notebooks.

# define a project directory
home = Path.home()
# create a processing directory
project_dir = home / "OST_Tutorials" / "4b"

# define aoi with helper function, i.e. get a buffered area around point coordinates
lat, lon = "-67", "-61"
aoi = vector.latlon_to_wkt(lat, lon, buffer_degree=1.5, envelope=True)

# define the start and end date
start = "2017-06-30"
end = "2017-08-31"

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

# trigger the search
s1_grd.search()
s1_grd.plot_inventory()

3 - Search Refinement

In order to create a time-series of multiple tracks, a pre-condition is that all tracks feature the same amount of acquistions within our Time of Interest. Let’s use some pandas syntax to see if this is the case:

df = s1_grd.inventory.pivot_table(index=["relativeorbit", "acquisitiondate"], aggfunc="size").reset_index()
df.pivot_table(index="relativeorbit", aggfunc="size").reset_index()

As in most cases, we do not fulfill this pre-condition. By considering all tracks, our time-series would need to be reduced to 5 acquisitions. However, images taken over track 9 are not necessary, since our AOI is fully covered by the other 2 tracks.

As already mentioned in OST Tutorial 2, the refine_inventory method takes care of those issues and prepares the inventory in a way that it is suitable for across-track mosaic time-series. This includes the splitting of images by orbit direction and polarzation mode in the first place. In addition, it checks if some tracks can be excluded because all the others fully overlap the AOI. In this way we reduce the amount of images to process, while optimising for our later time-series processing. See OST Tutorial 2 for full explanation and arguments.

# do the refinement
s1_grd.refine_inventory()

The output of the refinement procedure gives some infos, e.g. the exclusion of track 9. At the very end it summarises the information. Since in our case we only have imagery acquired in descending orbit and HH polarization, we see that 10 mosaics in time can be created. Another important infomation is the key defiend by orbit direction and polarisation, i.e. DESCENDING_HH. We will need this to select the refined inventory stored in the refined_inventory_dict attribute of our class instance as follows:

# select the key from output of refinement command
key = "DESCENDING_HH"

# we wrap the information of the length of our refined inventory in a print statement
print(f"The refined inventory holds {len(s1_grd.refined_inventory_dict[key])} acquisitions to process.")

# we plot the full Inventory on a map
s1_grd.plot_inventory(s1_grd.refined_inventory_dict[key], transparency=0.05)

3 - Data download

Here we download the data. It is best to use the ASF data mirror.

s1_grd.download(s1_grd.refined_inventory_dict[key], concurrent=10)

5* - Customise ARD parameters

1. For data reduction we lower the resolution to 100 meters.

2. This will also reduce speckle, so we do not need it neither.

3. We do not care about Layover and Shadow for this example, since there is anyway no high-resolution DEM for Antarctica that could provide sufficient inforamtion.

4. Our time-series output will be scaled to dB scale for better stretching and visualisation

5. We further reduce the amount of data by converting from 32-bit float to unsigned 8-bit integer data type

6. Our AOI is only a rough seletion for the Area of Interest. We do not want to cut it to the boundaries to see the full data provided by the selected imagery.

s1_grd.ard_parameters["single_ARD"]["resolution"] = 100
s1_grd.ard_parameters["single_ARD"]["remove_mt_speckle"] = False
s1_grd.ard_parameters["single_ARD"]["create_ls_mask"] = False
s1_grd.ard_parameters["single_ARD"]["dem"]["dem_name"] = "GETASSE30"
s1_grd.ard_parameters["time-series_ARD"]["to_db"] = True
s1_grd.ard_parameters["time-series_ARD"]["dtype_output"] = "uint8"
s1_grd.ard_parameters["mosaic"]["cut_to_aoi"] = False

pprint(s1_grd.ard_parameters)

6* - Produce Timeseries Mosaics

Now the creation of our mosaic time-series is just a single command away.

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

7* - Creation of an animated GIF of a timeseries

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

from ost.helpers import raster

# define the timeseries folder for which the animation should be created
ts_folder = s1_grd.processing_dir / "Mosaic" / "Timeseries"

raster.create_timeseries_animation(
    ts_folder,
    ["bs.HH"],
    s1_grd.processing_dir,
    shrink_factor=15,
    add_dates=True,
    duration=0.25,
)