Monitoring Urban Sprawl with Sentinel-1 SAR Data Using Pixel Value and Polarization Thresholding on Creodias

In this article you will utilize thresholding techniques and compare the values of VV and VH polarizations to analyze urban areas. As an illustration, you will attempt to download Sentinel-1 images containing data of the urban area of Warsaw, Poland.

The analysis of Synthetic Aperture Radar (SAR) images enables the detection of variations in radar wave reflections, which is crucial for identifying urbanized areas. SAR images exhibit distinct reflection patterns that allow for the delineation of urban areas, including agglomerations and smaller urban settlements. Urbanized areas stand out in these images primarily due to the intense reflection of radar waves by man-made structures, roads, and other artificial elements of infrastructure, facilitating their clear identification contrasted with natural terrain. Such analysis of areas like the Warsaw agglomeration can serve for

  • monitoring city growth and development,

  • urban planning,

  • assessment of natural and anthropogenic threats, and

  • environmental management.

Example result from the original workflow:

Combined VV and VH threshold mask for the Warsaw area

It clearly depicts the Warsaw agglomeration area, adjacent to the region between the Warsaw and Łódź agglomerations. This space is currently planned for development to connect both agglomerations. In the central part of the frame, there are plans for establishing a large airport complex.

What Are We Going To Do

  • Access data from Creodias

  • Extract urban areas by thresholding

  • Provide quantitative analysis by analyzing individual pixels

Prerequisites

No. 1 EODATA access

Obtain credentials that can list and download Sentinel-1 measurements from the EODATA S3 service available to your account. Use the access method for your execution environment:

No. 2 Catalogue and S3 connection settings

The query requires an OData catalogue with Collection, Attributes, ContentDate and S3Path fields and the OData.CSC.Intersects spatial filter. The selected product’s S3Path must refer to the bucket accessible through your S3 endpoint. A STAC URL cannot be used in this OData workflow.

Use these CREODIAS settings when prompted, or set the environment variables before starting Jupyter:

Connection settings

Environment variable

Value

EODATA_CATALOGUE_URL

https://datahub.creodias.eu/odata/v1/Products

EODATA_S3_ENDPOINT

https://eodata.cloudferro.com

EODATA_S3_BUCKET

eodata

EODATA_S3_REGION

Leave unset; press Enter at the signing-region prompt

These values follow EOData Catalogue API Manual on Creodias and How to access EODATA using boto3 on Creodias. The documented boto3 configuration does not set a separate signing region. Use EODATA credentials authorized for this S3 endpoint.

No. 3 Jupyter Notebook or JupyterLab with Python 3.9 or newer

Use an environment with network access to the selected catalogue and S3 endpoint. It needs enough memory for two complete Sentinel-1 measurement arrays and their derived arrays, plus disk space for downloaded TIFFs and results. Follow the Jupyter installation instructions if Jupyter is not installed.

No. 4 Python packages

Install the following packages in the Python environment used by your notebook:

python -m pip install requests boto3 numpy rasterio matplotlib

Download urban_sprawl.ipynb and run its cells in order. The notebook and the code below use the same connection settings. The S3 setup reads EODATA_ACCESS_KEY and EODATA_SECRET_KEY from the environment, or requests them at hidden prompts if they are unset. Keep real keys out of saved cells.

The stored figures and numerical output illustrate an earlier run; your selected product and results may differ.

Step 1 Accessing data from the cloud

1. Imports and Variable Declarations

The query selects an IW GRDH product with VV/VH polarization in the date range and Warsaw area of interest.

import os
from pathlib import Path
from getpass import getpass

import boto3
from botocore.config import Config
import requests
import numpy as np
import rasterio
import matplotlib.pyplot as plt

# Use the OData Products URL documented by your data provider.
CATALOGUE_URL = (
    os.environ.get("EODATA_CATALOGUE_URL")
    or input("OData Products URL (ending in /odata/v1/Products): ")
).strip().rstrip("/")
if not CATALOGUE_URL.startswith("https://") or not CATALOGUE_URL.endswith("/odata/v1/Products"):
    raise ValueError("Provide an HTTPS OData Products URL ending in /odata/v1/Products.")

# Warsaw area of interest.
collection = "SENTINEL-1"
product_type = "GRD"
start_date = "2024-05-15T00:00:00.000Z"
end_date = "2024-05-17T00:05:00.000Z"
polygon_coords = (
    "20.851215 52.261859, 21.271151 52.261859, "
    "21.271151 52.017798, 20.851215 52.017798, 20.851215 52.261859"
)
query_params = {
    "$filter": (
        f"Collection/Name eq '{collection}' and "
        "contains(Name,'_IW_GRDH_1SDV_') and "
        "Attributes/OData.CSC.StringAttribute/any(att:att/Name eq 'productType' "
        f"and att/OData.CSC.StringAttribute/Value eq '{product_type}') and "
        f"ContentDate/Start gt {start_date} and ContentDate/End lt {end_date} and "
        f"OData.CSC.Intersects(area=geography'SRID=4326;POLYGON(({polygon_coords}))')"
    ),
    "$orderby": "ContentDate/Start asc",
    "$top": 1,
}

2. S3 Configuration and Function to List TIFF Files

S3_ENDPOINT = (
    os.environ.get("EODATA_S3_ENDPOINT") or input("EODATA S3 endpoint (https://...): ")
).strip().rstrip("/")
S3_BUCKET = (os.environ.get("EODATA_S3_BUCKET") or input("EODATA bucket [eodata]: ").strip() or "eodata").strip()
S3_REGION = (os.environ.get("EODATA_S3_REGION") or input("S3 signing region (blank for SDK default): ")).strip() or None
if not S3_ENDPOINT.startswith("https://"):
    raise ValueError("Provide the HTTPS endpoint for your EODATA S3 service.")

# Use credentials authorized to list and download the selected EODATA objects.
# Environment variables are optional; otherwise enter the keys at hidden prompts.
access_key = os.environ.get("EODATA_ACCESS_KEY") or getpass("EODATA access key: ")
secret_key = os.environ.get("EODATA_SECRET_KEY") or getpass("EODATA secret key: ")
if not access_key or not secret_key:
    raise ValueError("Provide both EODATA access keys authorized for the configured S3 service.")

s3 = boto3.client(
    "s3",
    aws_access_key_id=access_key,
    aws_secret_access_key=secret_key,
    endpoint_url=S3_ENDPOINT,
    region_name=S3_REGION,
    config=Config(s3={"addressing_style": "path"}),
)

def list_tiff_files_in_s3(prefix):
    tiff_files = []
    paginator = s3.get_paginator("list_objects_v2")
    for page in paginator.paginate(Bucket=S3_BUCKET, Prefix=prefix):
        for obj in page.get("Contents", []):
            if obj["Key"].lower().endswith((".tif", ".tiff")):
                tiff_files.append(obj["Key"])
    return sorted(tiff_files)

3. Fetching Product Path from API Response

response = requests.get(CATALOGUE_URL, params=query_params, timeout=90)
response.raise_for_status()
products = response.json().get("value", [])
if not products:
    raise RuntimeError("No Sentinel-1 products matched the date range and area of interest.")

product = products[0]
product_name = product["Name"]
s3_path = product["S3Path"]
bucket_prefix = f"/{S3_BUCKET}/"
if not s3_path.startswith(bucket_prefix):
    raise ValueError(f"The catalogue path {s3_path!r} does not match bucket {S3_BUCKET!r}.")
product_path = s3_path.removeprefix(bucket_prefix).rstrip("/") + "/"
print(product_name)

4. Downloading TIFF Files from S3 and Saving Locally

# Keep each selected product in its own directory to avoid mixing acquisitions.
output_dir = Path("sen_1") / product_name
output_dir.mkdir(parents=True, exist_ok=True)
tiff_files = list_tiff_files_in_s3(product_path + "measurement/")
if not tiff_files:
    raise RuntimeError("No measurement TIFF files were found for the selected product.")

for key in tiff_files:
    file_name = Path(key).name
    local_file_path = output_dir / file_name
    s3.download_file(S3_BUCKET, key, str(local_file_path))
    print(f"Downloaded: {file_name}")

Downloaded: s1a-iw-grd-vh-20240516t162754-20240516t162819-053899-068d12-002.tiff
Downloaded: s1a-iw-grd-vv-20240516t162754-20240516t162819-053899-068d12-001.tiff

Step 2 Urban areas extraction by thresholding

1. Reading and preprocessing data

Read the measurements with rasterio from the product directory created above. Keep the original profiles and any ground control points (GCPs) for exporting the results.

def find_files_with_names(folder_path):
    files = sorted(Path(folder_path).glob("*"))
    vv_files = [p for p in files if "-vv-" in p.name.lower() and p.suffix.lower() in (".tif", ".tiff")]
    vh_files = [p for p in files if "-vh-" in p.name.lower() and p.suffix.lower() in (".tif", ".tiff")]
    if len(vv_files) != 1 or len(vh_files) != 1:
        raise ValueError("The product directory must contain exactly one VV and one VH measurement TIFF.")
    return vv_files[0], vh_files[0]

input_file_vv, input_file_vh = find_files_with_names(output_dir)
print(input_file_vv)

with rasterio.open(input_file_vv) as src_vv:
    vv_band = src_vv.read(1)
    profile_vv = src_vv.profile
    gcps_vv, gcps_crs_vv = src_vv.gcps

with rasterio.open(input_file_vh) as src_vh:
    vh_band = src_vh.read(1)
    profile_vh = src_vh.profile
    gcps_vh, gcps_crs_vh = src_vh.gcps

if vv_band.shape != vh_band.shape:
    raise ValueError("VV and VH measurements must have matching dimensions.")
# Retain the source profiles and ground control points; do not invent a CRS or transform.

2. Pixel value conversion to decibels

After reading the data, we can convert the pixel values to decibels. Data from Sentinel-1 is converted to the decibel scale (dB) to better visualize differences in the intensity of radar wave reflection, which aids in the analysis of urban areas due to their specific reflection characteristics.

The highest peak of the VV polarization histogram is in the range of 15-20 dB, suggesting that most surfaces in the studied area have moderate radar backscatter. The threshold of 23.5 dB is chosen to highlight pixels with the highest radar backscatter values, which may correspond to areas with intense urban infrastructure or other features with high reflectivity.

The highest peak of the VH polarization histogram is in the range of 13-18 dB, suggesting that this channel is more sensitive to surfaces with lower radar backscatter. The threshold of 21.5 dB for the VH channel is appropriate to identify pixels with higher backscatter, which helps to better differentiate terrain structures with slightly lower but still significant backscatter.

def convert_to_db(image):
    # Create a copy of the image to avoid modifying the original array
    image = np.array(image, dtype=np.float32, copy=True)

    # Replace non-positive values with a small positive value to avoid taking log of zero or negative
    image[image <= 0] = np.finfo(float).eps

    # Calculate dB values
    with np.errstate(divide='ignore', invalid='ignore'):
        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)

3. Data analysis and thresholding

The next step involves plotting histograms of the decibel pixel values for both the VV and VH bands. These histograms provide insights into the distribution of radar reflection intensities, aiding in understanding the characteristics of the observed area. The chosen thresholds, 23.5 dB for VV and 21.5 dB for VH, are selected based on prior knowledge or analysis requirements specific to urban areas to segment regions of interest in the images. The thresholding method applied here converts pixel values above the threshold to 1 and those below to 0, aiding in the delineation of features or areas with specific radar reflection characteristics typical of urban environments.

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()

# 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 = 23.5
threshold_value_vh = 21.5

# 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)
Histograms of the logarithmic VV and VH pixel values

4. Images plotting

Next on the subplots we can visualize the differences in radio signal reflection for various polarization types. By displaying the original and decibel-transformed images for both VV and VH polarizations, we can observe distinct features and intensities in urban areas. The change in signal reflection in the original images is barely visible, which proves the usefulness of using a decibel scale. Furthermore, by combining information from both polarizations through thresholding, we enhance our ability to mask urbanized regions, providing a more comprehensive understanding of the observed area’s urban characteristics.

# 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()
Original, logarithmic and thresholded VV and VH measurements 
# Display the combined threshold image using a specific colormap
plt.figure(figsize=(18, 9))
im = plt.imshow(combined_threshold, cmap='twilight_shifted')
plt.title('Combined Threshold Image')
cbar = plt.colorbar(im)
cbar.set_label('Signal Strength')
plt.axis('off')
plt.show()
Combined VV and VH threshold mask for the Warsaw area

This method of data analysis facilitates straightforward masking of urbanized areas. In the processed image, the Warsaw agglomeration is clearly visible on the right side. A star-shaped accumulation of population is apparent in Warsaw along major transportation corridors such as railways and expressways. Smaller urban centers are visible in the image. Furthermore, it can be observed that the number of masked pixels accumulates as we approach the city center, indicating denser urban development in the downtown area compared to the suburbs

5. Quantitative analysis

After processing our satellite images, we can conduct quantitative analysis. The characteristics of the satellite data enable us to analyze individual pixel values, which is crucial for detailed and statistical studies. The area calculation assumes 10 × 10 metre pixels and treats the threshold mask as an urban-area proxy. It is an illustrative estimate, not a validated land-cover classification.

# Calculate the number of urban pixels by summing up the binary values in the combined threshold image
urban_pixels = np.sum(combined_threshold)

# Define the area covered by a single pixel (in square meters)
pixel_area = 10 * 10  # Assuming each pixel represents a 10m x 10m area

# Calculate the total urban area in square meters
urban_area = urban_pixels * pixel_area

# Calculate the total number of pixels in the combined threshold image
total_pixels = combined_threshold.size

# Calculate the percentage of urbanized area in the image
urban_percentage = (urban_pixels / total_pixels) * 100

# Convert urban area from square meters to hectares
urban_area_hectares = urban_area / 10_000

# Convert total area from square meters to hectares
total_area_hectares = (total_pixels * pixel_area) / 10_000

print(f"Percentage of Urbanized Area: {urban_percentage:.2f}%")
print(f"Urbanized Area: {urban_area_hectares:.2f} ha")
print(f"Total Image Area: {total_area_hectares:.2f} ha")

Percentage of Urbanized Area: 9.08%
Urbanized Area: 401728.18 ha
Total Image Area: 4426150.50 ha

6. Save the results

The output files retain the original profiles and ground control points. These steps do not perform terrain correction or reproject the measurements.

output_file_db_vv = "vv_db_urban_warsaw.tif"
output_file_db_vh = "vh_db_urban_warsaw.tif"
output_file_threshold_combined = "warsaw_urban_threshold_combined.tif"

def write_result(path, data, profile, gcps, gcps_crs):
    output_profile = profile.copy()
    output_profile.update(dtype=data.dtype.name, count=1, nodata=None)
    with rasterio.open(path, "w", **output_profile) as dst:
        dst.write(data, 1)
        if gcps:
            dst.gcps = (gcps, gcps_crs)

write_result(output_file_db_vv, vv_band_db.astype(np.float32), profile_vv, gcps_vv, gcps_crs_vv)
write_result(output_file_db_vh, vh_band_db.astype(np.float32), profile_vh, gcps_vh, gcps_crs_vh)
write_result(output_file_threshold_combined, combined_threshold, profile_vv, gcps_vv, gcps_crs_vv)

What To Do Next

The following article uses similar techniques for processing Sentinel-5P data on air pollution:

Processing Sentinel-5P data on air pollution using Jupyter Notebook on Creodias.

Also of interest:

Analyzing and monitoring floods using Python and Sentinel-2 satellite imagery on Creodias.