Sentinel-1 SAR for Flood Mapping in the Philippines on CopPhil

What is going to be covered?

This notebook demonstrates how to use Sentinel-1 products to detect floods by comparing pre-event and post-event acquisitions.

We focus on major flood events in the Philippines (e.g., The Camarines Sur Province 30th november 2018 as shown in the code).

The workflow uses the CopPhil Infrastructure API and the /eodata S3 bucket to discover and access Sentinel-1 .SAFE products, extract VV and VH bands, convert to dB, apply thresholds, and map water/flood extent.

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

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-5P mission

Page Sentinel-5P mission shows basic information on Sentinel-5P 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.

import requests
import json
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os
import boto3
import rasterio

Connect with s3 storage - eodata

We first establish a connection to the CopPhil S3 storage provided

by CloudFerro.

This bucket (/eodata) hosts the Sentinel-1 archive, accessible

directly as .SAFE folders.

If the connection succeeds, we can later open Sentinel-1 products directly from /eodata paths without downloading them locally.

#Enter your eodata credentials here:
access_key = ''
secret_key = ''
host = 'https://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

Search Sentinel-1 Products (Camarines Sur, December 2018)

We create a search query for the CopPhil OData API to find

Sentinel-1 products that cover

Camarines Sur Province during the flood event in **late December

2018**.

Build a query

In order to find desired products it is needed to determinate some

specific filters:

  • collection_name: Sentinel-1 - product_type: GRD-COG - operationalMode = IW - aoi: extent (coordinates) of the area of interest (WGS84)

  • search_period_start: time range - start date

  • search_period_end: time range - end date

# OData search parameters: The Camarines Sur Province 30th november 2018
catalogue_odata_url = "https://catalogue.infra.copphil.philsa.gov.ph/odata/v1"
collection_name = "SENTINEL-1"
product_type = "GRD-COG"
operationalMode = "IW"

# (WKT, EPSG:4326). This small box straddles flood-affected river plains.
aoi = "POLYGON((123.05 14.00, 124.40 14.00, 124.40 13.00, 123.05 13.00, 123.05 14.00))"

# Flood window around peak inundation reports
search_period_start = "2018-12-18T00:00:00.000Z"
search_period_end   = "2018-12-30T23:59:59.000Z"

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

print(search_query)

https://catalogue.infra.copphil.philsa.gov.ph/odata/v1/Products?$filter=Collection/Name eq 'SENTINEL-1' and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' and att/OData.CSC.StringAttribute/Value eq 'GRD-COG') and Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'operationalMode' and att/OData.CSC.StringAttribute/Value eq 'IW') and OData.CSC.Intersects(area=geography'SRID=4326;POLYGON((123.05 14.00, 124.40 14.00, 124.40 13.00, 123.05 13.00, 123.05 14.00))') and ContentDate/Start gt 2018-12-18T00:00:00.000Z and ContentDate/Start lt 2018-12-30T23:59:59.000Z&$top=1000

Retrieve and Inspect Search Results

We send the query and parse the JSON response into a DataFrame for easier inspection. This gives us an overview of how many Sentinel-1 images from the search query are available.

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: 8

Create list of paths

Parse JSON data, extracts S3Path and Name values, and combines them into a list of strings formatted as name/path. It uses list comprehensions with zip to pair paths and names efficiently and then prints 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'])

path_list

['/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/29/S1B_IW_GRDH_1SSV_20181229T213005_20181229T213034_014262_01A86C_7391_COG.SAFE',
 '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/28/S1A_IW_GRDH_1SDV_20181228T213833_20181228T213902_025231_02CA19_F51D_COG.SAFE',
 '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/30/S1B_IW_GRDH_1SDV_20181230T095728_20181230T095753_014270_01A8AB_233C_COG.SAFE',
 '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/23/S1A_IW_GRDH_1SDV_20181223T213027_20181223T213056_025158_02C76E_B094_COG.SAFE',
 '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/23/S1A_IW_GRDH_1SDV_20181223T213056_20181223T213121_025158_02C76E_5D90_COG.SAFE',
 '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/28/S1A_IW_GRDH_1SDV_20181228T213902_20181228T213927_025231_02CA19_A477_COG.SAFE',
 '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/24/S1A_IW_GRDH_1SSV_20181224T095801_20181224T095832_025166_02C7B2_C662_COG.SAFE',
 '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/18/S1B_IW_GRDH_1SDV_20181218T095728_20181218T095753_014095_01A2D5_FF40_COG.SAFE']

And input the image folder path into output_dir, for extraction of vv and vh bands:

# Specify the directory where the files are located
output_dir = '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/18/S1B_IW_GRDH_1SDV_20181218T095728_20181218T095753_014095_01A2D5_FF40_COG.SAFE'
# We can either use the full path above or use the path_list from the search results
#output_dir = path_list[7]

# Function to find files containing specific substrings (e.g., "vv" and "vh") in their names
def find_files_with_names(folder_path, *names):
    vv_path = None
    vh_path = None

    # Traverse the directory tree starting at folder_path
    for root, dirs, files in os.walk(folder_path):
        # Look for 'measurement' folder in any subfolder
        if 'measurement' in dirs:
            measurement_folder = os.path.join(root, 'measurement')
            print(f"Found measurement folder: {measurement_folder}")

            # Look through the files in the "measurement" folder
            for file_name in os.listdir(measurement_folder):
                print(f"Checking file: {file_name}")  # Debugging line
                if "vv" in file_name.lower() and vv_path is None:  # Find the "vv" file (case-insensitive)
                    vv_path = os.path.join(measurement_folder, file_name)
                    print(f"Found vv file: {file_name}")  # Debugging line
                elif "vh" in file_name.lower() and vh_path is None:  # Find the "vh" file (case-insensitive)
                    vh_path = os.path.join(measurement_folder, file_name)
                    print(f"Found vh file: {file_name}")  # Debugging line
            break  # Exit after finding the files

    return vv_path, vh_path

# Find the "vv" and "vh" files in the specified directory
input_file_vv, input_file_vh = find_files_with_names(output_dir, "vv", "vh")

# Check if both vv and vh files were found
if input_file_vv and input_file_vh:
    # Open the "vv" file and read its data and profile
    with rasterio.open(input_file_vv) as src_vv:
        vv_band = src_vv.read(1)
        profile_vv = src_vv.profile

    # Open the "vh" file and read its data and profile
    with rasterio.open(input_file_vh) as src_vh:
        vh_band = src_vh.read(1)
        profile_vh = src_vh.profile

    print("VV file path:", input_file_vv)
    print("VH file path:", input_file_vh)
else:
    print("VV or VH files were not found.")

Found measurement folder: /eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/18/S1B_IW_GRDH_1SDV_20181218T095728_20181218T095753_014095_01A2D5_FF40_COG.SAFE/measurement
Checking file: s1b-iw-grd-vh-20181218t095728-20181218t095753-014095-01a2d5-002-cog.tiff
Found vh file: s1b-iw-grd-vh-20181218t095728-20181218t095753-014095-01a2d5-002-cog.tiff
Checking file: s1b-iw-grd-vv-20181218t095728-20181218t095753-014095-01a2d5-001-cog.tiff
Found vv file: s1b-iw-grd-vv-20181218t095728-20181218t095753-014095-01a2d5-001-cog.tiff
VV file path: /eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/18/S1B_IW_GRDH_1SDV_20181218T095728_20181218T095753_014095_01A2D5_FF40_COG.SAFE/measurement/s1b-iw-grd-vv-20181218t095728-20181218t095753-014095-01a2d5-001-cog.tiff
VH file path: /eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/18/S1B_IW_GRDH_1SDV_20181218T095728_20181218T095753_014095_01A2D5_FF40_COG.SAFE/measurement/s1b-iw-grd-vh-20181218t095728-20181218t095753-014095-01a2d5-002-cog.tiff

Process Sentinel-1 data

Convert images values to db scale

Now we convert the images to dB scale for easy thresholding:

def convert_to_db(image):
    # Replace zeros with a small positive value to avoid taking log of zero
    image[image == 0] = np.finfo(float).eps

    # Calculate dB values
    image_db = 10 * np.log10(image)

    # Replace infinite values with the maximum finite value
    max_value = np.nanmax(image_db[np.isfinite(image_db)])
    image_db[np.isinf(image_db)] = max_value

    return image_db

vv_band_db = convert_to_db(vv_band)
vh_band_db = convert_to_db(vh_band)

/tmp/ipykernel_1222422/220689011.py:6: RuntimeWarning: divide by zero encountered in log10
  image_db = 10 * np.log10(image)

Thresholding

Parse JSON data, extracts S3Path and Name values, and combines them into a list of strings formatted as name/path. It uses list comprehensions with zip to pair paths and names efficiently and then prints the first ten merged results.

With the images in dB scale we can make histograms of VV and VH backscatter to determine thresholds for separating water:

fig, axs = plt.subplots(1, 2, figsize=(18, 5))

# Histogram of VV Image in dB
axs[0].hist(vv_band_db.flatten(), bins=50, color='blue', alpha=0.7, label='VV dB')
axs[0].set_title('Histogram of VV Image in dB')
axs[0].set_xlabel('Radar Reflection Value (dB)')
axs[0].set_ylabel('Number of Pixels')
axs[0].legend()
axs[0].grid(True, linestyle='--', alpha=0.5)

# Histogram of VH Image in dB
axs[1].hist(vh_band_db.flatten(), bins=50, color='red', alpha=0.7, label='VH dB')
axs[1].set_title('Histogram of VH Image in dB')
axs[1].set_xlabel('Radar Reflection Value (dB)')
axs[1].set_ylabel('Number of Pixels')
axs[1].legend()
axs[1].grid(True, linestyle='--', alpha=0.5)

plt.show()
../_images/Flood-mapping-with-Sentinel-1-SAR-data-in-JupyterLab-on-CopPhil_18_0.png 
# Function to apply a threshold to an image
# Pixels greater than the threshold are set to 1, others to 0
def thresholding(image, threshold):
    return (image > threshold).astype(np.uint8)

# Define threshold values for VV and VH bands
threshold_value_vv = 21.5
threshold_value_vh = 18
# Apply thresholding to the VV and VH band
vv_band_threshold = thresholding(vv_band_db, threshold_value_vv)
vh_band_threshold = thresholding(vh_band_db, threshold_value_vh)

# Combine the thresholded images using a logical AND operation
# The result is 1 where both VV and VH are above their respective thresholds
combined_threshold = np.logical_and(vv_band_threshold, vh_band_threshold).astype(np.uint8)

Visualization

Before the flood

Now we can visualize original, dB scale, and images after thresholding:

# Create a figure with a 3x2 grid of subplots, setting the overall figure size
fig, ax = plt.subplots(3, 2, figsize=(18, 18))

# Display the original VV image
ax[0, 1].imshow(vv_band, cmap='viridis')
ax[0, 1].set_title('Original VV Image')

# Display the VV image in decibels
ax[1, 1].imshow(vv_band_db, cmap='viridis')
ax[1, 1].set_title('VV Image in dB')

# Display the VV image after thresholding
ax[2, 1].imshow(vv_band_threshold, cmap='inferno')
ax[2, 1].set_title('VV Image after Thresholding')

# Display the original VH image
ax[0, 0].imshow(vh_band, cmap='viridis')
ax[0, 0].set_title('Original VH Image')

# Display the VH image in decibels
ax[1, 0].imshow(vh_band_db, cmap='viridis')
ax[1, 0].set_title('VH Image in dB')

# Display the VH image after thresholding
ax[2, 0].imshow(vh_band_threshold, cmap='inferno')
ax[2, 0].set_title('VH Image after Thresholding')

# Turn off the axes for all subplots to improve the visual presentation
for i in range(2):
    for j in range(3):
        ax[j, i].axis('off')

# Show the complete figure
plt.show()
../_images/Flood-mapping-with-Sentinel-1-SAR-data-in-JupyterLab-on-CopPhil_22_0.png

We can take shape of combined_threshold image to determing the zoom area:

combined_threshold_flipped = np.flipud(combined_threshold)# Flippig the image to match the orientation of a normal map
print("Array shape:", combined_threshold_flipped.shape)
print("Max y index:", combined_threshold_flipped.shape[0] - 1)  # rows
print("Max x index:", combined_threshold_flipped.shape[1] - 1)  # cols

Array shape: (16798, 25189)
Max y index: 16797
Max x index: 25188

And plot the zoomed in image:

# Coordinates to zoom into
x_zoom = 11000
y_zoom = 8000

# Define a zoom window size (e.g., 100x100 pixels around the (x, y) coordinates)
zoom_window_size = 4000  # You can adjust this to zoom in/out more

# Crop the image around the specified coordinates (make sure not to exceed bounds)
cropped_image = combined_threshold_flipped[
    max(0, y_zoom - zoom_window_size): min(combined_threshold_flipped.shape[0], y_zoom + zoom_window_size),
    max(0, x_zoom - zoom_window_size): min(combined_threshold_flipped.shape[1], x_zoom + zoom_window_size)
]

# Display the zoomed-in region
plt.figure(figsize=(9, 9))
plt.imshow(cropped_image, cmap='viridis')

# Add axis with pixel coordinates
plt.title('Zoomed-In Region of Libmanan and Bicol Rivers (18.12.2018)')
plt.axis('on')  # Turn on the axis

plt.show()
../_images/Flood-mapping-with-Sentinel-1-SAR-data-in-JupyterLab-on-CopPhil_26_0.png

After the flood

Now we do the same steps for post flood image:

# Specify the directory where the files are located
output_dir_post = '/eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/30/S1B_IW_GRDH_1SDV_20181230T095728_20181230T095753_014270_01A8AB_233C_COG.SAFE'
# Find the "vv" and "vh" files in the specified directory
input_file_vv_post, input_file_vh_post = find_files_with_names(output_dir_post, "vv", "vh")

# Check if both vv and vh files were found
if input_file_vv_post and input_file_vh_post:
    # Open the "vv" file and read its data and profile
    with rasterio.open(input_file_vv_post) as src_vv_post:
        vv_band_post = src_vv_post.read(1)
        profile_vv_post = src_vv_post.profile
    # Open the "vh" file and read its data and profile
    with rasterio.open(input_file_vh_post) as src_vh_post:
        vh_band_post = src_vh_post.read(1)
        profile_vh_post = src_vh_post.profile
#Convert to dB scale
vv_band_db_post = convert_to_db(vv_band_post)
vh_band_db_post = convert_to_db(vh_band_post)
# Define threshold values for VV and VH bands
threshold_value_vv_post = 23
threshold_value_vh_post = 17.5
# Apply thresholding to the VV and VH band
vv_band_threshold_post = thresholding(vv_band_db_post, threshold_value_vv_post)
vh_band_threshold_post = thresholding(vh_band_db_post, threshold_value_vh_post)
# Combine the thresholded images
combined_threshold_post = np.logical_and(vv_band_threshold_post, vh_band_threshold_post).astype(np.uint8)
# Flippig the image to match the orientation of a normal map
combined_threshold_flipped_post = np.flipud(combined_threshold_post)

Found measurement folder: /eodata/Sentinel-1/SAR/IW_GRDH_1S-COG/2018/12/30/S1B_IW_GRDH_1SDV_20181230T095728_20181230T095753_014270_01A8AB_233C_COG.SAFE/measurement
Checking file: s1b-iw-grd-vh-20181230t095728-20181230t095753-014270-01a8ab-002-cog.tiff
Found vh file: s1b-iw-grd-vh-20181230t095728-20181230t095753-014270-01a8ab-002-cog.tiff
Checking file: s1b-iw-grd-vv-20181230t095728-20181230t095753-014270-01a8ab-001-cog.tiff
Found vv file: s1b-iw-grd-vv-20181230t095728-20181230t095753-014270-01a8ab-001-cog.tiff

/tmp/ipykernel_1222422/220689011.py:6: RuntimeWarning: divide by zero encountered in log10
  image_db = 10 * np.log10(image)

Now we can plot the same area in the image post flood:

# Coordinates to zoom into
x_zoom = 11000
y_zoom = 8000

# Define a zoom window size (e.g., 100x100 pixels around the (x, y) coordinates)
zoom_window_size = 4000  # You can adjust this to zoom in/out more

# Crop the image around the specified coordinates (make sure not to exceed bounds)
cropped_image_post = combined_threshold_flipped_post[
    max(0, y_zoom - zoom_window_size): min(combined_threshold_flipped_post.shape[0], y_zoom + zoom_window_size),
    max(0, x_zoom - zoom_window_size): min(combined_threshold_flipped_post.shape[1], x_zoom + zoom_window_size)
]

# Display the zoomed-in region
plt.figure(figsize=(9, 9))
plt.imshow(cropped_image_post, cmap='viridis')

# Add axis with pixel coordinates
plt.title('Zoomed-In Region of Libmanan River (30.12.2018)')
plt.axis('on')  # Turn on the axis

plt.show()
../_images/Flood-mapping-with-Sentinel-1-SAR-data-in-JupyterLab-on-CopPhil_30_0.png

And plot the thresholded images pre and post flood side by side to better visualize changes:

fig, ax = plt.subplots(1, 2, figsize=(18, 18))#setting grid for images pre and post flood
ax[0].imshow(cropped_image, cmap='viridis')
ax[0].set_title('Combined Thresholded Image Before Flood Event (Dec 18, 2018)')
ax[1].imshow(cropped_image_post, cmap='viridis')
ax[1].set_title('Combined Thresholded Image During The Flood Event (Dec 30, 2018)')
plt.show()
../_images/Flood-mapping-with-Sentinel-1-SAR-data-in-JupyterLab-on-CopPhil_32_0.png

Summary

This notebook demonstrated how to use Sentinel-1 data for flood monitoring. It showcased the potential of radar imagery, which can be used regardless of cloud cover, making it particularly valuable for tracking and monitoring natural phenomena and hazards.