Top of Atmosphere Reflectance on Sentinel 3
Top of Atmosphere (TOA) Reflectance is a unitless measurement which provides the ratio of radiation reflected to the incident solar radiation on a given surface. It that can be computed from satellite measured spectral radiance using the mean solar spectral irradiance and the solar zenith angle.
When comparing EO images from different sensors, there are two advantages of using TOA reflectance instead of TOA spectral radiance.
First, it removes the cosine effect of different solar zenith angles due to the time difference between data acquisitions. Second, TOA reflectance compensates for different values of the solar irradiance arising from spectral band differences.
The TOA reflectance of the Earth is computed according to the equation:
where
Implementation
The code to do the TOA Reflectance conversion for sentinel 3 optical instruments is available at:
https://github.com/fer-marino/sentinel3_optical
It contains two scripts to convert the bands radiance to TOA reflectance and produce an RGB composition. It is a python script that uses numpy, scipy and xarray to do its job.
OLCI
For the OLCI instrument, each radiance band is located in a separated file named OaXX_radiance.nc, containing a variable with the same name.
red = open_dataset(product + "/Oa10_radiance.nc")["Oa10_radiance"].values[:]
green = open_dataset(product + "/Oa05_radiance.nc")["Oa05_radiance"].values[:]
blue = open_dataset(product + "/Oa03_radiance.nc")["Oa03_radiance"].values[:]
Here we use xarray library to open each file. Then we access the variable we need via the square bracket operator, then copy the whole content into a new variable. This basically copy in memory the entire band in a single step, which is way more efficient than accessing the single pixels throught this API. Another advantage is the automatic reformatting performed by the netCDF driver from uint16 to float32, applying a scale factor and an offset.
We do the same to ancillary data we need to apply the formula above:
instrument_data = open_dataset(product + "/instrument_data.nc")
tie_geometries = open_dataset(product + "/tie_geometries.nc")
detector_index = instrument_data["detector_index"].values[:]
sza = tie_geometries["SZA"].values[:]
solar_flux = instrument_data["solar_flux"].values[:]
We are extracting also a variable named detector_index. This is a 2D array that for each pixel gives us the detector number that was used for that pixel. We need this information to extract the solar flux at each detector.
Finally, we process each pixel applying the formula above. We skip invalid pixels doing a check on the detector_index variable, but a similar check could have been done on measurement data.
detector = int(detector_index[y][x])
angle = sza[y][int(x / 64)]
sol_flux_red = solar_flux[9][detector] # band 10 for red, has index 9
sol_flux_green = solar_flux[4][detector] # band 5 for green, has index 4
sol_flux_blue = solar_flux[2][detector] # band 3 for blue, has index 2
img[y][x][0] = max(0, min(1., math.pi * (red[y][x] / sol_flux_red) / math.cos(math.radians(angle))))*2 -1
img[y][x][1] = max(0, min(1., math.pi * (green[y][x] / sol_flux_green) / math.cos(math.radians(angle))))*2 -1
img[y][x][2] = max(0, min(1., math.pi * (blue[y][x] / sol_flux_blue) / math.cos(math.radians(angle))))*2 -1
The solar zenith angle variable does not fully cover the entire product, it has the same number of rows but a subsampled number of columns (1 out of 64). The reflectance is expected to go from zero to one, anyway this is not always true, especially over clouds. Anyway we do not care about clouds so we saturate anything above 1. We also rescale the value from the range [0, 1] to [-1,1] as scipy expect this range in float images.
Done? Almost. You can already see the results of this processing, anyway you may notice that for a better results it is necessary to apply some filter to make it looks nicer. The easier way is to apply a histogram equalization filter:
img_equalized = exposure.equalize_hist(img, nbins=512)
This filter will adjust the contrast of the image to make it looks nicer. The nbins parameter is the number of bins of the histogram the algorithm will build. Increasing this number will give you more color precision, but it will also make it slower and more memory time consuming during computation.
this sequence shows the reflectance, the radiance and the difference between the two
The results clearly show the correction applied by the sun angle, which make brighter the top right, and darker the bottom left. The differences also show that reflectance conversion has a limited effect on clouds which are usually almost saturated in the visible channels.
SLSTR
Also for SLSTR we have separate files for each band, here anyway we have also different views to pick the data from. The algorithm is almost the same as OLCI with few differences mostly because of the different data structure. The only remark is the Sun Zenith angle variable structure: it is a single 2D array which have the nadir view and the oblique view placed next to each other, so you need to do some math to place each SZA on the radiance grid. In this composition we decided to combine multiple bands to obtain more natural colors:
- Red: (S3 + S5) / 2
- Green: S3
- Blue: (S1 + S5) /2
comparison between radiance and TOA reflectance
In this RGB composition we see that the reflectance conversion corrects the excessive blue intensity on the radiance, due to the solar irradiance division.
Very interesting and complete “tutorial” . Looking at the SKIMAGE module of scipy you could also apply the rescale_intensity function that sets the contrast stretching depending on low and high saturation levels. It needs some more lines of code but the user can modulate the contrast with that.