OLCI Level 1B Reduced Resolution - Sentinel-3

import destinelab as deauth

import requests
import json
import os
from getpass import getpass

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}"}

HDA_STAC_ENDPOINT="https://hda.data.destination-earth.eu/stac/v2"

COLLECTION_ID = "EO.EUM.DAT.SENTINEL-3.OL_1_ERR___"

response = requests.post(HDA_STAC_ENDPOINT+"/search", headers=auth_headers, json={
    "collections": [COLLECTION_ID],
    "datetime": "2024-06-25T00:00:00Z/2024-06-30T00:00:00Z",
    "bbox":  [10,53,30,66]
})
if(response.status_code!= 200):
    (print(response.text))
response.raise_for_status()

from IPython.display import JSON

product = response.json()["features"][0]
JSON(product)

from tqdm import tqdm
import time

assets = ["downloadLink"]

for asset in assets:
    download_url = product["assets"][asset]["href"]
    print(download_url)
    filename = product["id"]
    print(filename)
    response = requests.get(download_url, headers=auth_headers)
    total_size = int(response.headers.get("content-length", 0))

    print(f"downloading {filename}")

    with tqdm(total=total_size, unit="B", unit_scale=True) as progress_bar:
        with open(filename, 'wb') as f:
            for data in response.iter_content(1024):
                progress_bar.update(len(data))
                f.write(data)

import zipfile
import os
import shutil
import uuid

# 1) Create a temporary extraction directory
tmp_dir = f"tmp_extract_{uuid.uuid4().hex}"
os.makedirs(tmp_dir, exist_ok=True)

# 2) Extract zip content safely inside the temporary folder
with zipfile.ZipFile(filename, 'r') as z:
    z.extractall(tmp_dir)

# 3) Identify the Sentinel‑3 product folder inside the temp folder
entries = os.listdir(tmp_dir)
product_dirs = [e for e in entries if os.path.isdir(os.path.join(tmp_dir, e))]

if not product_dirs:
    raise RuntimeError("No directory found inside the ZIP. Unexpected structure!")

product_name = product_dirs[0]  # usually only one
product_path_tmp = os.path.join(tmp_dir, product_name)

# 4) rename the extracted folder
final_name=product_path_tmp+".SEN3"
shutil.move(product_path_tmp, final_name)

print("Extraction complete!")
print("Product available at:", final_name)

pip install --quiet satpy pyspectral

from packaging.version import Version
import satpy
from datetime import datetime
from satpy import find_files_and_readers
from satpy.scene import Scene
from satpy import available_readers

print(satpy.__version__)
if Version(satpy.__version__) < Version("0.57"):
    from satpy.composites import GenericCompositor
    from satpy.writers import to_image
    from satpy.resample import get_area_def
elif Version(satpy.__version__) == Version("0.57"):
    from satpy.composites import GenericCompositor
    from satpy.area import get_area_def
else:
    from satpy.composites.core import GenericCompositor
    from satpy.area import get_area_def
from satpy import available_readers

import pyresample
import pyspectral

import warnings
warnings.filterwarnings("ignore")
warnings.simplefilter(action = "ignore", category = RuntimeWarning)

from datetime import datetime
from satpy import find_files_and_readers


files = find_files_and_readers(
    sensor="olci",
    start_time=datetime(2024, 6, 25, 0, 0),
    end_time=datetime(2024, 6, 30, 0, 0),
    base_dir=tmp_dir+"/",
    reader="olci_l1b",
)
scn = Scene(filenames=files, reader="olci_l1b")

scn = Scene(filenames=files)

scn.available_dataset_names()

# load bands 
scn.load(['humidity','total_ozone'])
scn['humidity']

import matplotlib.pyplot as plt
from pyresample.kd_tree import resample_nearest
from pyresample.geometry import AreaDefinition

#area definition
area_id = 'worldeqc30km'
description = 'World in 3km, platecarree'
proj_id = 'eqc'
projection = {'proj': 'eqc', 'ellps': 'WGS84'}
width = 820
height = 410
area_extent = (-20037508.3428, -10018754.1714, 20037508.3428, 10018754.1714)
area_def = AreaDefinition(area_id, description, proj_id, projection,
                          width, height, area_extent)

#scene 
lons, lats = scn["total_ozone"].area.get_lonlats()
swath_def = pyresample.geometry.SwathDefinition(lons, lats)
total_ozone = scn["total_ozone"].data.compute()
result = resample_nearest(swath_def, total_ozone, area_def, radius_of_influence=20000, fill_value=None)

#cartopy
plt.rcParams['figure.figsize'] = [15, 15]
crs = area_def.to_cartopy_crs()
fig, ax = plt.subplots(subplot_kw=dict(projection=crs))
coastlines = ax.coastlines()  
ax.set_global()

#plot
img = plt.imshow(result, transform=crs, extent=crs.bounds, origin='upper')
# Calculate (height_of_image / width_of_image)
im_ratio = result.shape[0]/result.shape[1]
 
# Plot vertical colorbar
plt.colorbar(fraction=0.047*im_ratio)
plt.show()