Help needed for OLCI sentinel 3

[HafezAhmad]

I hope you are doing well. I am currently working with OLCI in Google Earth Engine and have a few questions:
My study area near coast. does not cover srtm dem.
I have set up the following link https://py6s.readthedocs.io/en/latest/installation.html
I am trying to convert radiance to reflectance, here is the code, please review it.

import ee
from Py6S import *
import datetime
import math
import os
import sys
import pytz
from pysolar.solar import get_altitude

# Initialize Earth Engine
ee.Initialize()

# Define the dates to download
datesstart = ['2024-03-10', '2024-03-11', '2024-03-12', '2024-03-13', '2024-06-24', '2024-06-25', '2024-06-26']

# Define the study area
study_area = ee.FeatureCollection('projects/ee-wfahydrologymsu100/assets/WMS')
geometry = study_area.geometry()

# Function to calculate solar zenith angle
def calculate_solar_z(geometry, scene_date):
   latitude = geometry.centroid().coordinates().get(1).getInfo()
   longitude = geometry.centroid().coordinates().get(0).getInfo()
   solar_altitude = get_altitude(latitude, longitude, scene_date)
   solar_z = 90 - solar_altitude
   return solar_z

# Function to calculate surface reflectance
def surface_reflectance(bandname, image):
   # Instantiate the 6S model
   s = SixS(r"C:\Users\hafez\build\6SV\1.1\6SV1.1\sixsV1.1")
   
   # Extract image metadata
   info = image.getInfo()['properties']
   scene_date = datetime.datetime.utcfromtimestamp(info['system:time_start'] / 1000)
   scene_date = scene_date.replace(tzinfo=pytz.UTC)
   
   # Calculate solar zenith angle
   solar_z = calculate_solar_z(geometry, scene_date)
   
   # Get atmospheric data
   h2o = Atmospheric.water(geometry, ee.Date(scene_date)).getInfo()
   o3 = Atmospheric.ozone(geometry, ee.Date(scene_date)).getInfo()
   aot = Atmospheric.aerosol(geometry, ee.Date(scene_date)).getInfo()
   
   # Set the atmospheric constituents
   s.atmos_profile = AtmosProfile.UserWaterAndOzone(h2o, o3)
   s.aero_profile = AeroProfile.Continental
   s.aot550 = aot
   
   # Set the Earth-Sun-satellite geometry
   s.geometry = Geometry.User()
   s.geometry.view_z = 0               # always NADIR
   s.geometry.solar_z = solar_z        # solar zenith angle
   s.geometry.month = scene_date.month # month and day used for Earth-Sun distance
   s.geometry.day = scene_date.day     # month and day used for Earth-Sun distance
   s.altitudes.set_sensor_satellite_level()
   
   # Define the spectral response function for Sentinel-3 OLCI
   def spectralResponseFunction(bandname):
       bandSelect = {
           'Oa02_radiance': PredefinedWavelengths.S3A_OLCI_02,
           'Oa03_radiance': PredefinedWavelengths.S3A_OLCI_03,
           'Oa04_radiance': PredefinedWavelengths.S3A_OLCI_04,
           'Oa05_radiance': PredefinedWavelengths.S3A_OLCI_05,
           'Oa06_radiance': PredefinedWavelengths.S3A_OLCI_06,
           'Oa07_radiance': PredefinedWavelengths.S3A_OLCI_07,
           'Oa08_radiance': PredefinedWavelengths.S3A_OLCI_08,
           'Oa09_radiance': PredefinedWavelengths.S3A_OLCI_09,
           'Oa11_radiance': PredefinedWavelengths.S3A_OLCI_11,
           'Oa12_radiance': PredefinedWavelengths.S3A_OLCI_12,
           'Oa13_radiance': PredefinedWavelengths.S3A_OLCI_13,
       }
       return Wavelength(bandSelect[bandname])

   def toa_to_rad(bandname, image):
       rad = image.select(bandname)
       return rad

   # Run 6S for this waveband
   s.wavelength = spectralResponseFunction(bandname)
   s.run()
   
   # Extract  Edir 
   Edir = s.outputs.direct_solar_irradiance  # Direct solar irradiance
   Edif = s.outputs.diffuse_solar_irradiance  # Diffuse solar irradiance
   Lp = s.outputs.atmospheric_intrinsic_radiance  # Path radiance
   absorb = s.outputs.trans['global_gas'].upward  # Absorption transmissivity
   scatter = s.outputs.trans['total_scattering'].upward  # Scattering transmissivity
   tau2 = absorb * scatter  # Total transmissivity
   
   # Radiance to surface reflectance
   rad = toa_to_rad(bandname, image)
   ref = rad.subtract(Lp).multiply(math.pi).divide(tau2 * (Edir + Edif))
   
   return ref

# Process each date and download the image
for date in datesstart:
   date_obj = ee.Date(date)
   asset_name = 'S3OLCI_' + date_obj.format('YYYY_MM_dd').getInfo()
   
   # Get the Sentinel-3 OLCI image
   S3_image = ee.Image(
       ee.ImageCollection('COPERNICUS/S3/OLCI')
       .filterBounds(geometry)
       .filterDate(date_obj, date_obj.advance(1, 'day'))
       .sort('system:time_start')
       .first()
   )
   
   # Calculate surface reflectance for each band
   bands = ['Oa02_radiance', 'Oa03_radiance', 'Oa04_radiance', 'Oa05_radiance', 'Oa06_radiance', 'Oa07_radiance',
            'Oa08_radiance', 'Oa09_radiance', 'Oa11_radiance', 'Oa12_radiance', 'Oa13_radiance']
   
   surface_reflectance_images = [surface_reflectance(band, S3_image) for band in bands]
   
   # Combine all surface reflectance bands into a single image
   combined_reflectance = ee.Image.cat(surface_reflectance_images)
   
   # Clip the image to the study area geometry
   clipped_reflectance = combined_reflectance.clip(geometry)
   
   # Define the export parameters
   export_params = {
       'image': clipped_reflectance.toFloat(),
       'description': asset_name,
       'scale': 10,  # Adjust the scale as needed
       'region': geometry,
       'fileFormat': 'GeoTIFF'
   }
   
   # Export the image to Google Drive
   task = ee.batch.Export.image.toDrive(**export_params)
   task.start()
   
   # Print the export task status
   print(f'Exporting {asset_name}...')

print('All exports queued successfully.')

`


Best,
Hafez