Sentinel-2 ToA reflectance to radiance¶
Author: Gonzalo Mateo-García
In this notebook we show how to convert ToA reflectance from Sentinel-2 to radiance using the metadata fetched from the GEE or directly computing the conversion factors.
import ee
import matplotlib.pyplot as plt
from georeader import plot
from shapely.geometry import box
from georeader.readers import ee_image
from datetime import datetime, timezone
from rasterio import Affine
from georeader.readers import S2_SAFE_reader
from georeader import reflectance
import os
os.environ["GS_NO_SIGN_REQUEST"] = "YES"
# ee.Authenticate()
ee.Initialize()
Fetch a S2 image from GEE¶
collection_name = "COPERNICUS/S2_HARMONIZED"
tile = "S2A_MSIL1C_20240417T064631_N0510_R020_T40RCN_20240417T091941"
img_col = ee.ImageCollection(collection_name)
image = img_col.filter(ee.Filter.eq("PRODUCT_ID", tile)).first()
info_img = image.getInfo()
# projgee = image.select("B2").projection().getInfo()
aoi = box(55.325, 25.225, 55.415, 25.28)
bands_0bname = list(S2_SAFE_reader.BANDS_S2_L1C)
bands = [b.replace("B0","B") for b in bands_0bname]
crs = info_img["bands"][1]["crs"]
transform = Affine(*info_img["bands"][1]["crs_transform"])
img_local = ee_image.export_image(info_img['id'],
crs=crs, transform=transform,
bands_gee=bands,
geometry=aoi)
img_local
Fetch same image from Google bucket¶
s2_safe_folder = S2_SAFE_reader.s2_public_bucket_path(tile, check_exists=False)
print(f"File is located at: {s2_safe_folder}")
s2obj = S2_SAFE_reader.s2loader(s2_safe_folder, out_res=10)
s2obj
Convert DN to reflectances¶
toa_refl = img_local / 10_000
Fetch metadata to do toa reflectance to radiance conversion¶
solar_irradiance_gee = {k:v/1_000. for k, v in info_img['properties'].items() if k.startswith("SOLAR_IRRADIANCE")}
for b in bands:
print(f"{b} -> {solar_irradiance_gee[f'SOLAR_IRRADIANCE_{b}']}")
date_of_aquisition = datetime.fromtimestamp(info_img['properties']["system:time_start"]/ 1000).replace(tzinfo=timezone.utc)
date_of_aquisition
mean_solar_zenith_angle = info_img['properties']['MEAN_SOLAR_ZENITH_ANGLE']
U = info_img['properties']["REFLECTANCE_CONVERSION_CORRECTION"]
mean_solar_zenith_angle, U
import numpy as np
observation_date_corr_factor = np.pi / (np.cos(mean_solar_zenith_angle/180*np.pi) * U)
print(observation_date_corr_factor)
rad = reflectance.reflectance_to_radiance(toa_refl,
solar_irradiance = [solar_irradiance_gee[f"SOLAR_IRRADIANCE_{b}"] for b in bands],
observation_date_corr_factor=observation_date_corr_factor)
rad
Metadata from the SAFE file¶
s2obj.read_metadata_tl()
solar_irr_safe = s2obj.solar_irradiance()
U_save = s2obj.scale_factor_U()
print(f"SZA: {s2obj.mean_sza} U: {U_save}")
solar_irr_safe
for b in bands_0bname:
print(f"{b} -> {solar_irr_safe[b]}")
Calculate Irradiance and correction factor using georeader¶
satellite = tile.split("_")[0] # S2A, S2B or S2C
srf_satellite = S2_SAFE_reader.read_srf(satellite)
irradiances = reflectance.integrated_irradiance(srf_satellite) / 1_000 # Convert from mW/m2/sr/nm to W/m2/sr/nm
solar_irradiance_calc = dict(zip(srf_satellite.columns, irradiances))
solar_irradiance_calc
Show differences in solar irradiance¶
for b in bands_0bname:
bnozero = b.replace("B0", "B")
print(f"{b} -> GEE: {solar_irradiance_gee[f'SOLAR_IRRADIANCE_{bnozero}']} SAFE: {solar_irr_safe[b]:.4f} Calc: {solar_irradiance_calc[b]:.4f}")
observation_date_corr_factor_calc = reflectance.observation_date_correction_factor(img_local.footprint("EPSG:4326").centroid.coords[0],
date_of_acquisition=date_of_aquisition)
observation_date_corr_factor_calc
rad2 = reflectance.reflectance_to_radiance(toa_refl,
solar_irradiance = [solar_irradiance_calc[b] for b in bands_0bname],
observation_date_corr_factor=observation_date_corr_factor)
rad2
Show relative differences between the calculated and official factors¶
for i, b in enumerate(bands_0bname):
rad_b_or = rad.isel({"band":i})
rad_b_2 = rad2.isel({"band":i})
ratio = rad_b_or.values / rad_b_2.values
mean_ratio = np.mean(ratio)
std_ratio = np.std(ratio)
print(f"{b} {mean_ratio:.4f} {std_ratio:.6f}")
Licence¶
The georeader package is published under a GNU Lesser GPL v3 licence
georeader tutorials and notebooks are released under a Creative Commons non-commercial licence.
If you find this work useful please cite:
@article{portales-julia_global_2023,
title = {Global flood extent segmentation in optical satellite images},
volume = {13},
issn = {2045-2322},
doi = {10.1038/s41598-023-47595-7},
number = {1},
urldate = {2023-11-30},
journal = {Scientific Reports},
author = {Portalés-Julià, Enrique and Mateo-García, Gonzalo and Purcell, Cormac and Gómez-Chova, Luis},
month = nov,
year = {2023},
pages = {20316},
}