Tracking typhoons over the Philippines using JupyterLab on CopPhil

Earth Observation data are coming from satellites and can amount to gigabytes and petabytes of information, captured as high resolution images. To speed up the search and visual recognition, each eodata high resolution image has its scaled-down version, called quicklook.

What Are We Going To Cover

This article demonstrates how to use Sentinel-3 data to track typhoons over the Philippines and its surrounding areas. It is divided into four main sections:

  • Searching for Sentinel-3 Level 1 products (OL_1_EFR___) using specific filters.

  • Accessing data via the S3 protocol.

  • Plotting quicklooks of the found images to identify where typhoons have been observed.

  • Creating maps based on the selected data.

Prerequisites

No. 1 Access to CopPhil site

You need a CopPhil hosting account, available at https://infra.copphil.philsa.gov.ph.

No. 2 Access to JupyterLab

JupyterLab used in this article is available here: https://jupyter.infra.copphil.philsa.gov.ph/hub/login?next=%2Fhub%2F.

No. 3 Working knowledge of JupyterLab

See article Introduction to JupyterLab on CopPhil

Also of interest: Python 3 kernel and Geo science kernel in JupyterLab on CopPhil

No. 4 Using API

This is an introductory example of obtaining EODATA using API calls, in Python:

Download satellite data using API on CopPhil

No. 5 Information on Sentinel-3 Level 1 mission

Page Sentinel-3 mission shows basic information on Sentinel-3 mission, which is used in this article as a source of information.

No. 6 Have credentials to access EODATA ready

To obtain the necessary credentials to access EODATA, follow this article:

How to get credentials to access EODATA on CopPhil

You can generate up to 200 different keys per account.

Prepare your environment

Upload necessary python libraries.

# HTTP requests
import requests

# JSON parser
import json

# data manipulation
import pandas as pd
import numpy as np
import xarray as xr
from netCDF4 import Dataset
import rasterio
from rasterio.errors import NotGeoreferencedWarning
import warnings

# image manipulation
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# HTTP requests
import requests

# File manipulation
import os
import tempfile

# Interact with s3
import boto3
from botocore.exceptions import ClientError
from io import BytesIO

Connect with s3 storage - eodata

To access eodata, you must have a pair of credentials called access key and secret key. Use Prerequisite No. 7 to obtain them, then enter them as values for variables access_key and secret_key, respectively.

If connection is successful, message

Connection to bucket eodata successful. Objects found.

will be printed.

access_key = ''
secret_key = ''
host = 'eodata.infra.copphil.philsa.gov.ph'

# Initialize S3 resource
s3 = boto3.resource(
    's3',
    aws_access_key_id=access_key,
    aws_secret_access_key=secret_key,
    endpoint_url=host,
)

# Verify connection by accessing the 'eodata' bucket
bucket_name = 'eodata'
try:
    bucket = s3.Bucket(bucket_name)
    # Check if the bucket is accessible by listing its contents
    objects = list(bucket.objects.limit(1))  # Limit to 1 object for testing
    if objects:
        print(f"Connection to bucket '{bucket_name}' successful. Objects found.")
    else:
        print(f"Connection to bucket '{bucket_name}' successful, but no objects found.")
except Exception as e:
    print(f"Connection to bucket '{bucket_name}' failed:", e)

Connection to bucket 'eodata' successful. Objects found.

Search for the Sentinel-3 products

Let us use principles from Prerequisite No. 4 to obtain data.

Build a query

In order to find desired products it is needed to specify the following filters:

collection_name

Sentinel-5

product_type

OL_1_EFR___

timeliness

NT (Non time critical)

aoi

extent (coordinates) of the area of interest (WGS84)

search_period_start

time range - start date

search_period_end

time range - end date

# Base URL of the product catalogue
catalogue_odata_url = "https://catalogue.infra.copphil.philsa.gov.ph/odata/v1"

# Search parameters
collection_name = "SENTINEL-3"
product_type = "OL_1_EFR___"
timeliness = "NT"
aoi = "POLYGON((116.87759 19.423629, 116.640126 4.793402, 149.362676 4.414686, 147.605441 18.570422, 116.87759 19.423629))"
search_period_start = "2024-11-13T00:00:00.000Z"
search_period_end = "2024-11-17T23:59:59.000Z"

# Construct the OData query
search_query = f"{catalogue_odata_url}/Products?$filter=Collection/Name eq '{collection_name}' and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' and att/Value eq '{product_type}') and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'timeliness' and att/Value eq '{timeliness}') and OData.CSC.Intersects(Footprint=geography'SRID=4326;{aoi}') and ContentDate/Start gt {search_period_start} and ContentDate/Start lt {search_period_end}&$orderby=ContentDate/Start asc&$top=1000"


# Print the formatted URL
print(search_query)

https://catalogue.infra.copphil.philsa.gov.ph/odata/v1/Products?$filter=Collection/Name eq 'SENTINEL-3' and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' and att/Value eq 'OL_1_EFR___') and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'timeliness' and att/Value eq 'NT') and OData.CSC.Intersects(Footprint=geography'SRID=4326;POLYGON((116.87759 19.423629, 116.640126 4.793402, 149.362676 4.414686, 147.605441 18.570422, 116.87759 19.423629))') and ContentDate/Start gt 2024-11-13T00:00:00.000Z and ContentDate/Start lt 2024-11-17T23:59:59.000Z&$orderby=ContentDate/Start asc&$top=1000

Inspect results of the request

Show the number of found products which fits into filters requirements and get a list of names and id’s of images which fits the query.

response = requests.get(search_query).json()
result = pd.DataFrame.from_dict(response["value"])

id_list = []
for element in response['value']:
    id_value = element['Id']
    id_list.append(id_value)


name_list = []
for element in response['value']:
    name_value = element['Name']
    name_list.append(name_value)


print(f'Number of products which meets query filters: {len(result)}')

Number of products which meets query filters: 26

Create list of paths

To create a list of paths, the following code will

  • parse JSON data, extract S3Path and Name values, then

  • combine them into a list of strings formatted as name/path.

It will use list comprehensions with zip to pair paths and names efficiently and then print the first ten merged results.

products = json.loads(requests.get(search_query).text)

path_list = []
for item in products['value']:
    path_list.append(item['S3Path'])

for path in path_list:
    print(path)

/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/13/S3B_OL_1_EFR____20241113T011826_20241113T012126_20241113T100413_0179_099_359_2700_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/13/S3B_OL_1_EFR____20241113T012126_20241113T012426_20241113T100429_0179_099_359_2880_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/13/S3A_OL_1_EFR____20241113T015411_20241113T015711_20241114T021208_0179_119_117_2520_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/13/S3A_OL_1_EFR____20241113T015711_20241113T020011_20241114T021344_0179_119_117_2700_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/13/S3A_OL_1_EFR____20241113T020011_20241113T020311_20241114T021343_0179_119_117_2880_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/14/S3B_OL_1_EFR____20241114T005515_20241114T005815_20241114T095715_0179_099_373_2880_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/14/S3A_OL_1_EFR____20241114T012800_20241114T013100_20241115T014852_0180_119_131_2520_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/14/S3A_OL_1_EFR____20241114T013100_20241114T013400_20241115T014852_0179_119_131_2700_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/14/S3A_OL_1_EFR____20241114T013400_20241114T013700_20241115T014854_0179_119_131_2880_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/14/S3B_OL_1_EFR____20241114T023315_20241114T023615_20241114T101201_0179_099_374_2700_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/15/S3A_OL_1_EFR____20241115T010449_20241115T010749_20241116T012820_0179_119_145_2700_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/15/S3A_OL_1_EFR____20241115T010749_20241115T011049_20241116T012822_0179_119_145_2880_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/15/S3B_OL_1_EFR____20241115T020404_20241115T020704_20241115T111213_0179_100_003_2520_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/15/S3B_OL_1_EFR____20241115T020704_20241115T021004_20241115T111224_0179_100_003_2700_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/15/S3B_OL_1_EFR____20241115T021004_20241115T021304_20241115T111235_0180_100_003_2880_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/16/S3B_OL_1_EFR____20241116T013754_20241116T014054_20241116T100657_0179_100_017_2520_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/16/S3B_OL_1_EFR____20241116T014054_20241116T014354_20241116T100708_0179_100_017_2700_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/16/S3B_OL_1_EFR____20241116T014354_20241116T014654_20241116T100720_0179_100_017_2880_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/16/S3A_OL_1_EFR____20241116T021638_20241116T021938_20241117T023828_0179_119_160_2520_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/16/S3A_OL_1_EFR____20241116T021938_20241116T022238_20241117T023828_0179_119_160_2700_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/16/S3A_OL_1_EFR____20241116T022238_20241116T022538_20241117T023835_0179_119_160_2880_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/17/S3B_OL_1_EFR____20241117T011443_20241117T011743_20241117T100707_0180_100_031_2700_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/17/S3B_OL_1_EFR____20241117T011743_20241117T012043_20241117T100717_0179_100_031_2880_PS2_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/17/S3A_OL_1_EFR____20241117T015027_20241117T015327_20241118T020958_0179_119_174_2520_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/17/S3A_OL_1_EFR____20241117T015327_20241117T015627_20241118T020958_0179_119_174_2700_PS1_O_NT_004.SEN3
/eodata/Sentinel-3/OLCI/OL_1_EFR___/2024/11/17/S3A_OL_1_EFR____20241117T015627_20241117T015927_20241118T021037_0179_119_174_2880_PS1_O_NT_004.SEN3

Verify whether the number of paths is the same as the number of products which met query filters.

if len(path_list) == len(path_list):
    print("Number of product which were found is the same as the number of paths.")
else:
    print("ERROR. Number of product which were found is not the same as the number of paths.")

Number of product which were found is the same as the number of paths.

Display and inspect quicklooks

Let’s visualize the quicklook images from Sentinel-3 data by displaying them alongside the directory names, making it easier to inspect the data visually.

Directory Iteration

Loop through the all found directories in path_list to locate the quicklook.jpg image file.

Quicklook Image Plotting

For each directory containing the quicklook.jpg, the image is loaded and displayed using matplotlib, with the directory path shown in the plot’s title.

All quicklook images will be plotted along with their position in path_list, helping to locate them. This will assist in identifying and visually pinpointing images where typhoons are visible.

# Loop through the directories in path_list (you can adjust this as needed)
for idx, path in enumerate(path_list[:10]):
    # Define the path to the quicklook.jpg file
    quicklook_path = os.path.join(path, "quicklook.jpg")

    # Check if the quicklook.jpg file exists in the directory
    if os.path.exists(quicklook_path):
        # Load and plot the quicklook image using matplotlib
        img = mpimg.imread(quicklook_path)

        # Create the plot
        plt.figure(figsize=(8, 6))
        plt.imshow(img)
        plt.axis('off')
        plt.title(f"path_list #{idx}: {os.path.basename(path)}")

        # Display the plot
        plt.show()
    else:
        print(f"No quicklook.jpg found in {path}")
../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_0.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_1.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_2.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_3.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_4.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_5.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_6.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_7.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_8.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_13_9.png

Four images were identified with typhoons: 0, 7, 10, 16, and 24. Let’s visualize them in a more geographic context, using their coordinates.

Visualize images where typhoons were observed

Visualize Sentinel-3 RGB composite images, with geographic coordinates, to track typhoons observed in the Philippines. For the typhoon images, the coordinates are extracted, and each image is overlaid with the country’s borders.

Key steps:

  • Extract the relevant bands (red, green, blue) from Sentinel-3 data.

  • Normalize and stack the bands to create an RGB composite.

  • Plot the composite image with geographic coordinates and overlay the Philippines boundary.

  • Display the RGB composites for specific indices where typhoons were observed (0, 7, 10, 16, and 24).

warnings.filterwarnings("ignore", category=NotGeoreferencedWarning)

def download_needed_sen3_files(s3_path):
    key = s3_path.strip().lstrip("/").replace("eodata/", "", 1)
    local_dir = tempfile.mkdtemp(suffix=".SEN3")

    needed_files = [
        "Oa08_radiance.nc",
        "Oa06_radiance.nc",
        "Oa04_radiance.nc",
        "geo_coordinates.nc",
    ]

    for fname in needed_files:
        bucket.download_file(
            f"{key}/{fname}",
            os.path.join(local_dir, fname)
        )

    return local_dir

def extract_band(sen3_path, band_index):
    fname = f"Oa{band_index:02d}_radiance.nc"
    subdataset_path = os.path.join(sen3_path, fname)
    with rasterio.open(subdataset_path) as dataset:
        return dataset.read(1)

def extract_coordinates(sen3_path):
    geo_path = os.path.join(sen3_path, "geo_coordinates.nc")
    with Dataset(geo_path, mode="r") as geo_nc:
        lat = geo_nc.variables["latitude"][:]
        lon = geo_nc.variables["longitude"][:]
    return lat, lon

def create_rgb_composite(sen3_path):
    red_band = extract_band(sen3_path, 8)
    green_band = extract_band(sen3_path, 6)
    blue_band = extract_band(sen3_path, 4)

    def normalize_band(band):
        band = np.ma.filled(band, np.nan)
        p2, p98 = np.nanpercentile(band, (2, 98))
        return np.clip((band - p2) / (p98 - p2) * 255, 0, 255).astype(np.uint8)

    return np.stack([
        normalize_band(red_band),
        normalize_band(green_band),
        normalize_band(blue_band)
    ], axis=-1)

def plot_rgb_with_coordinates_and_border(s3_path):
    local_sen3 = download_needed_sen3_files(s3_path)

    rgb_image = create_rgb_composite(local_sen3)
    lat, lon = extract_coordinates(local_sen3)

    fig = plt.figure(figsize=(10, 8))
    ax = plt.axes(projection=ccrs.PlateCarree())

    ax.imshow(
        rgb_image,
        extent=(lon.min(), lon.max(), lat.min(), lat.max()),
        transform=ccrs.PlateCarree(),
        origin="upper"
    )

    ax.add_feature(cfeature.COASTLINE, linewidth=1, edgecolor="red")
    ax.add_feature(cfeature.BORDERS, linewidth=1, edgecolor="red")
    ax.add_feature(cfeature.LAND, facecolor="none", edgecolor="black")

    ax.set_extent(
        [lon.min(), lon.max(), lat.min(), lat.max()],
        crs=ccrs.PlateCarree()
    )

    ax.set_title(f"RGB Composite with Cartopy borders\n{os.path.basename(s3_path)}")
    plt.show()

indices_to_plot = [0, 7, 16, 24]

for idx in indices_to_plot:
    plot_rgb_with_coordinates_and_border(path_list[idx])
../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_16_0.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_16_1.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_16_2.png ../_images/Disaster-preparedness-and-response-system-in-JupyterLab-on-Eumetsat-Elasticity_16_3.png

Summary

This article demonstrated the process of visualizing Sentinel-3 data to track typhoons over the Philippines. The data was accessed, and quicklook images were analyzed to identify typhoon events. Using specific filters, relevant Sentinel-3 Level 1 products were found and extracted, and geographic visualizations were created to show the typhoon locations.