How to use HDA to find and download data for conducting monitoring of Śniadrwy lake

import requests
import zipfile
import io
import destinelab as deauth
from getpass import getpass

# 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 = 'output/'
# Geometry in form of a BBOX
bbox = [21.61868,53.66627,21.90926,53.82351]

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-2.MSI.L1C', 'EO.ESA.DAT.SENTINEL-2.MSI.L2A']

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

import rasterio
from osgeo import gdal, gdal_array, osr
import numpy as np
import os

def getFullPath(dir: str, resolution: int, band: str):
    if not os.path.isdir(dir):
        raise ValueError(f"Provided path does not exist: {dir}")
    elif resolution not in [10,20,60]:
        raise ValueError(f"Provided resolution does not exist: R{resolution}m")
    else:
        full_path = dir
        while True:
            content = os.listdir(full_path)
            if len(content) == 0:
                raise ValueError(f"Directory empty: {full_path}")
            elif len(content) == 1:
                if full_path[-1] != '/':
                    full_path = full_path + '/' + content[0]
                else:
                    full_path = full_path + content[0]
            else:
                if 'GRANULE' in content:
                    full_path = full_path + '/' + 'GRANULE'
                    break
                else:
                    raise ValueError(f"Unsupported dir architecture: {full_path}")
        full_path = full_path + '/' + os.listdir(full_path)[0]
        full_path = full_path + '/' + "IMG_DATA"
        if len(os.listdir(full_path)) == 3:
            full_path = full_path + '/' + f'R{resolution}m'
            images = os.listdir(full_path)
            for img in images:
                if band in img:
                    return full_path + '/' + img
            raise ValueError(f'No such band {band} in directory: {full_path}')
        else:
            images = os.listdir(full_path)
            for img in images:
                if band in img:
                    return full_path + '/' + img
            raise ValueError(f'No such band {band} in directory: {full_path}')

# Get transformation matrix from raster
def getTransform(pathToRaster):
    dataset = gdal.Open(pathToRaster)
    transformation = dataset.GetGeoTransform()
    return transformation

# Read raster and return pixels' values matrix as int16, new transformation matrix, crs
def readRaster(path, resolution, band):
    path = getFullPath(path, resolution, band)
    trans = getTransform(path) # trzeba zdefiniować który kanał
    raster, crs = rasterToMatrix(path)
    return raster.astype(np.int16), crs, trans

def rasterToMatrix(pathToRaster):
    with rasterio.open(pathToRaster) as src:
        matrix = src.read(1)
    return matrix, src.crs.to_epsg()

# Transform numpy's matrix to geotiff; pass new raster's filepath, matrix with pixels' values, gdal file type, transformation matrix, projection, nodata value
def npMatrixToGeotiff(filepath, matrix, gdalType, projection, transformMatrix, nodata = None):
    driver = gdal.GetDriverByName('Gtiff')
    if len(matrix.shape) > 2:
        (bandNr, yRes, xRes) = matrix.shape
        image = driver.Create(filepath, xRes, yRes, bandNr, gdalType)
        for b in range(bandNr):
            b = b + 1
            band = image.GetRasterBand(b)
            if nodata is not None:
                band.SetNoDataValue(nodata)
            band.WriteArray(matrix[b-1,:,:])
            band.FlushCache
    else:
        bandNr = 1
        (yRes, xRes) = matrix.shape
        image = driver.Create(filepath, xRes, yRes, bandNr, gdalType)
        print(type(image))
        band = image.GetRasterBand(bandNr)
        if nodata is not None:
            band.SetNoDataValue(nodata)
        band.WriteArray(matrix)
        band.FlushCache
    image.SetGeoTransform(transformMatrix)
    image.SetProjection(projection)
    del driver, image, band

# List of products' paths. If your output directory contains more than just downloaded products, please provide them in a list, just like in the commented lines below
#dataset = [output_dir+x for x in os.listdir(output_dir)]
dataset = [
     'output/S2B_MSIL1C_20230701T094549_N0509_R079_T34UEE_20230701T104205'
 ]
# Output directiry for new images
compution_output = 'output/ndwi_rasters'

# Iterating over single product
for item in dataset:
    # Reading name from path
    name = item.split('/')[-1]
    # Reading green band into matrix
    green = readRaster(item, 10, 'B03')
    # Reading NIR band into matrix
    nir = readRaster(item, 10, 'B08')
    # Calculating NDWI matrix
    ndwi = (green[0]-nir[0]) / (green[0]+nir[0])
    # Setting treshhold for water-containing pixels
    ndwi[ndwi >= 0] = 1
    ndwi[ndwi < 0] = 0
    # Creating SpatialReference object and setting it to match original's raster CRS
    projection = osr.SpatialReference()
    projection.ImportFromEPSG(green[1])
    # Creating raster from matrix in GeoTiff format
    npMatrixToGeotiff(f'{compution_output}/{name}.tif', ndwi, gdal_array.NumericTypeCodeToGDALTypeCode(np.float32), projection.ExportToWkt(), green[2])

img = rasterio.open('output/ndwi_rasters/S2B_MSIL1C_20230701T094549_N0509_R079_T34UEE_20230701T104205.tif')
from rasterio.plot import show
show(img)