Using HDA to Find and Download Data for Urban Area Monitoring Using Sentinel-1 Data

import datetime
import requests
import numpy as np
import rasterio
import matplotlib.pyplot as plt
import requests
import zipfile
from getpass import getpass
import io
from rasterio.mask import mask
import os
import destinelab as deauth

# Timerange of data that we want to recieve
start_date = '2023-07-01'
end_date = '2023-07-07'
# Output directory of our desired data
output_dir = 'sen_1/'
# Geometry in form of a BBOX
bbox = [20.8510, 52.0976, 21.2712, 52.3345,]

COLLECTIONS_URL = 'https://hda.data.destination-earth.eu/stac/collections'
SEARCH_URL = 'https://hda.data.destination-earth.eu/stac/search'

def get_stac_collections(api_url):
    response = requests.get(api_url)
    if response.status_code == 200:
        stac_data = response.json()['collections']
        collections = [x['id'] for x in stac_data]
        return collections
    else:
        print(f"Error: {response.status_code} - {response.text}")
        return None
    
get_stac_collections(COLLECTIONS_URL)

collections = ['EO.ESA.DAT.SENTINEL-1.L1_GRD']

DESP_USERNAME = input("Please input your DESP username or email: ")
DESP_PASSWORD = getpass("Please input your DESP password: ")

auth = deauth.AuthHandler(DESP_USERNAME, DESP_PASSWORD)
access_token = auth.get_token()
if access_token is not None:
    print("DEDL/DESP Access Token Obtained Successfully")
else:
    print("Failed to Obtain DEDL/DESP Access Token")

auth_headers = {"Authorization": f"Bearer {access_token}"}

def search_items(access_token: str, search_url: str, collection: str, 
                 bbox: list[float | int], start_date: str, end_date: str):
    body = {
        'datetime': f'{start_date}T00:00:00Z/{end_date}T23:59:59Z',
        'collections': [collection],
        'bbox': bbox
    }
    response = requests.post(search_url, json=body, headers={'Authorization': 'Bearer {}'.format(access_token)})
    if response.status_code != 200:
        print(f'Error in search request: {response.status_code} - {response.text}')
        return None
    else:
        print("Request successful! Reading data...")
        products_list = [(feature.get('assets').get('downloadLink').get('href'), feature.get('links')[0].get('title')) for feature in response.json().get('features', [])]
        return products_list

collections_items = []
for c in collections:
    collections_items.append(search_items(access_token, SEARCH_URL, c, bbox, start_date, end_date))

def hda_download(access_token: str, url: str, output: str):
    response = requests.get(url,stream=True,headers={'Authorization': 'Bearer {}'.format(access_token), 'Accept-Encoding': None})
    if response.status_code == 200:
        print('Downloading dataset...')
        with zipfile.ZipFile(io.BytesIO(response.content)) as z:
            z.extractall(output)
        print('The dataset has been downloaded to: {}'.format(output))
    else:
        print('Request Unsuccessful! Error-Code: {}'.format(response.status_code))

for collection in collections_items:
    for item in collection:
        url = item[0]
        product_id = item[1]
        download_path = output_dir + product_id
        hda_download(access_token, url, download_path)
        break

# Function to find files containing specific substrings ("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):
        if dirs:  # Check if there are subdirectories
            first_subfolder = dirs[0]
            measurement_folder = os.path.join(root, first_subfolder, "measurement")
            
            if os.path.isdir(measurement_folder):
                # Look through the files in the "measurement" folder
                for file_name in os.listdir(measurement_folder):
                    if "vv" in file_name:
                        vv_path = os.path.join(measurement_folder, file_name)
                    elif "vh" in file_name:
                        vh_path = os.path.join(measurement_folder, file_name)
                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")

# 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

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)

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)

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

# Display the combined threshold image using a specific colormap
plt.figure(figsize=(18, 9))  
plt.imshow(combined_threshold, cmap='twilight_shifted')
plt.title('Combined Threshold Image')  
plt.axis('off')  
plt.show()

# 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")

# Define output file names
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'

# Update the data type in the profile to float32 for dB conversion
profile_vv.update(dtype=rasterio.float32)
profile_vh.update(dtype=rasterio.float32)

# Write the result data in dB to the output file
with rasterio.open(output_file_db_vv, 'w', **profile_vv) as dst:
    dst.write(vv_band_db.astype(np.float32), 1)
    
with rasterio.open(output_file_db_vh, 'w', **profile_vh) as dst:
    dst.write(vh_band_db.astype(np.float32), 1)

# Update the data type in the VV profile to uint8 for the combined thresholded image
profile_vv.update(dtype=rasterio.uint8)

# Write the combined thresholded image data to the output file
with rasterio.open(output_file_threshold_combined, 'w', **profile_vv) as dst:
    dst.write(combined_threshold, 1)