Description of class S2RutOp

Step 1: Initialize function

This function defines the product reading, the setup of the uncertainty product output and the global parameters extraction. Specific tasks:

The source product is read (source_product) and verified as a valid Sentinel-2 product.
Parsing Sentinel-2 metadata at three levels: product_meta, datastrip_meta and granule_meta. It also identifies the spacecraft for either S2A or S2B units.
Gets the selected bands (toa_band_names) and verifies that it is not empty
Extraction of global metadata parameters: distance to sun (u_sun),reflectance quantification (quant) and sun zenith angle (source_sza). The later contains the per-pixel sun zenith angle for 10, 20 and 60m spatial resolutions.
Extraction of user parameters: coverage factor (k) and uncertainty contributors flags (unc_select).
Setup of the uncertainty product including L1C scene_width, scene_height and geolocation: rut_product = snappy.Product(self.source_product.getName() + '_rut', 'S2_RUT', scene_width, scene_height)
Add each one of the selected bands (unc_toa_band) in the uncertainty product with the correct name and Nodata value.
Different metadata fields ar populated and added to the uncertainty product.

Step 2: computeTile function

This function selects the specific band in a tile, extracts the required parameters and sets the uncertainty samples after calling the S2RutAlgo operator. Specific tasks:

Selects the source_band and toa_band_id.
Selects source_sza, cloudmask and cirrusmask according to the spatial resolution of the band.
Extracts the metadata parameters for the specific band: absolute calibration coefficient (a), sun irradiance (e_sun), noise model coefficients (alpha and beta), and estimated diffuser degradation (u_diff_temp).
Extracts the L1C samples toa_samples
Call the class S2RutAlgo for the per-pixel uncertainty calculation unc = self.rut_algo.unc_calculation(np.array(toa_samples, dtype=np.uint16), toa_band_id)
Obtains the rest of the S2 masks and combines those values with unc.
Set the samples in the destination tile

Description of class S2RutAlgo

Step 1: Undo reflectance conversion

This step converts the L1C reflectance in L1B CN counts.

cn = (self.a * self.e_sun * self.u_sun *
math.cos(math.radians(self.tecta)) / ( math.pi * self.quant)) *
band_data.

Step 2: Orthorectification process

This step evaluates the impact of the uncertainty propagation through the resampling process.

RUTv1 does not include this step. This is to be implemented in next versions of the RUT

Step 3: L1B uncertainty contributors: raw and dark signal

NOTE: Each the following uncertainty contributors are pre-checked and included if unc_select == True

Noise in % units
u_noise = 100 * np.sqrt(self.alpha ** 2 + self.beta * cn) / cn
Out-of field straylight - systematic in radiance units
u_stray_sys = 0.3 * rad_conf.Lref[band_id] / 100
Out-of field straylight - random in % units
u_stray_rand = rad_conf.u_stray_rand_all[band_id]
Crosstalk in radiance units
u_xtalk = rad_conf.u_xtalk_all[band_id]
Instrument quantisation. Defined in _init_ function.
self.u_ADC = 0.5
Predefined DS stability in LSB units
u_DS = rad_conf.u_DS_all[band_id]

Step 4: L1B uncertainty contributors: gamma correction

Predefined relative gains residual in % units

u_gamma = 0.4

Step 5: L1C uncertainty contributors: absolute calibration coefficient

Pre-flight diffuser reflectance knowledge in % units
u_diff_abs = rad_conf.u_diff_absarray[band_id]
u_diff_temp calculated in computeTile function(see above)
Cosine diffuser knowledge. Defined in _init_ function.
self.u_diff_cos = 0.4
Diffuser straylight residual correction. Defined in _init_ function.
self.u_diff_k = 0.3

Step 6: L1C uncertainty contributors: reflectance conversion

Impact of the reflectance factor digitisation in % units

u_ref_quant = 100 * (0.5 / math.sqrt(3)) / band_data

Step 7: Combine uncertainty contributors

Calculation of the terms u'_DS and u'_ADC
u_adc = (100 * self.u_ADC / math.sqrt(3)) / cn ; u_ds = (100 * u_DS) / cn
Calculation of the straylight term and conversion of the crosstalk into % units
u_stray = np.sqrt(u_stray_rand ** 2 + ((100 * self.a * u_xtalk) / cn) ** 2)
Calculation of the diffuser uncertainty
u_diff = math.sqrt(u_diff_abs ** 2 + self.u_diff_cos ** 2 + self.u_diff_k ** 2)
Calculation of the standard uncertainty
u_1sigma = np.sqrt(u_ref_quant ** 2 + u_gamma ** 2 + u_stray ** 2 + u_diff ** 2 + u_noise ** 2 + u_adc ** 2 + u_ds ** 2)
Calculation of the expanded uncertainty in % units. Conversion of the Out-of-Field Straylight - systematic in % units. Resulting uncertainty values multiplied by 10
u_expand = 10 * (self.u_diff_temp + ((100 * self.a * u_stray_sys) / cn) + self.k * u_1sigma)
Clipping all the values to range [0,250] for uncertainty image UINT8 type
u_ref = np.uint8(np.clip(u_expand, 0, 250))