TESTING A COMBINED MULTISPECTRAL-MULTITEMPORAL APPROACH FOR GETTING CLOUDLESS IMAGERY FOR SENTINEL-2

: Earth observation and land cover monitoring are among major applications for satellite data. However, the use of primary satellite information is often limited by clouds, cloud shadows, and haze, which generally contaminate optical imagery. For purposes of hazard assessment, for instance, such as flooding, drought, or seismic events, the availability of uncontaminated optical data is required. Different approaches exist for masking and replacing cloud/haze related contamination. However, most common algorithms take advantage by employing thermal data. Hence, we tested an algorithm suitable for optical imagery only. The approach combines a multispectral-multitemporal strategy to retrieve daytime cloudless and shadow-free imagery. While the approach has been explored for Landsat information, namely Landsat 5 TM and Landsat 8 OLI, here we aim at testing the suitability of the method for Sentinel-2 Multi-Spectral Instrument. A multitemporal stack, for the same image scene, is employed to retrieve a composite uncontaminated image over a temporal period of few months. Besides, in order to emphasize the effectiveness of optical imagery for monitoring post-disaster events, two temporal stages have been processed, before and after a critical seismic event occurred in Lombok Island, Indonesia, in summer 2018. The approach relies on a clouds and cloud shadows masking algorithm, based on spectral features, and a data reconstruction phase based on automatic selection of the most suitable pixels from a multitemporal stack. Results have been tested with uncontaminated image samples for the same scene. High accuracy is achieved.


INTRODUCTION
Today, the availability of remotely sensed data as provided by satellite missions is key information for Earth Observation (EO) and land use/land cover (LULC) monitoring. Moreover, the increasing concern about phenomena related to climate change and natural disasters such as flooding, drought, and seismic events, among others, is requiring more and more the availability and usability of such an information, either raw as well as processed. However, collecting primary remotely sensed information is not as immediate as desirable due to several atmospheric occurrences. Clouds, cloud shadows, and haze generally strongly contaminate optical remote-sensing images. Hence, for EO purposes, such contamination is intended as missing data and should be replaced or at least reduced (Colaninno, Marambio, & Roca, 2019). Different algorithms are available for automatically detect and mask clouds, and cloud shadows. As reported by Lin, Lai, Chen, & Chen (2014) reconstruction approaches are referred to three main categories. In-painting-based methods, which rely on data reconstruction by propagating the spectro-geometrical information of the same image to fill no-data from surrounding uncontaminated data (Lorenzi, Melgani, & Mercier, 2011). Multispectral-based methods, which rely on the analysis of spectral features of clouds across the visible, Near Infrared, and the thermal spectrum, for cloud detection and information restoration (Roy et al., 2008;Wang, Jin, Liang, Yan, & Peng, 2005). The automated cloud-cover assessment (ACCA) algorithm (Irish, Barker, Goward, & Arvidson, 2006), and the FMask (Zhu, Wang, & Woodcock, 2015) that also uses the view angle of the satellite sensor and the illuminating angle to assess possible cloud shadows, are among the main multispectral-based methods. While, multitemporal-based methods rely on temporal and spatial coherence to combine information from different time periods based either on threshold approach (Min Li, Soo Chin Liew, & Leong Keong Kwoh, 2003), regression tree (Helmer & Ruefenacht, 2005), or a contextual prediction approach to determine spectro-temporal relationships among image sequences (Benabdelkader & Melgani, 2008). This allows to get cloud-uncontaminated and non-shadowed pixels from images acquired at different times and reconstruct a cloud-free image scene. Besides, for reconstructing clouds contaminated imagery, geostatistical methods have been investigated that rely, for instance, on independent component analysis (Shen, Wang, Lv, & Qian, 2015), or interpolation methods based on ordinary cokriging and standardized ordinary co-kriging (C. Zhang, Li, & Travis, 2009). However, some limitations such as handling large cloud covered and shadowed areas in a heterogeneous landscape, or small-scale applications, as well as issues related to thermal response, or the need for effective clouds detection and masking algorithms, make the matter still challenging and necessary.
For certain purposes, the availability of cloudless optical data is fundamental. It has been widely demonstrated the effectiveness of such an information in post-disaster monitoring for instance. From 2000 on, satellite data has been widely investigated for disaster risk management (Bessis, Béquignon, & Mahmood, 2004;Kaku, 2019). Also, because over the last years, extreme earthquakes have struck extensive parts of the world, causing massive damages and destruction; monitoring earthquakes based on the use of satellite observations, both optical and radar, have gained interest globally ("Earthquakes from space: Earth observation for quantifying earthquake risks," 2017).
This work aims at testing the effectiveness of an approach that combines multispectral and multitemporal processes. First, masking clouds and cloud shadows is undertaken. Then, cloudless (uncontaminated) optical satellite images is obtained based on a multitemporal stack along one scene. We tested the suitability of the multispectral multitemporal cloud-free tool (M2CTool), as suggested by Colaninno, Marambio, & Roca (2019), and previously explored for Landsat data, for different mid-resolution satellite information. Here the test is assessed for Sentinel-2 MSI derived imagery. The accuracy of the masking step resulted quite high even at a higher spatial resolution with respect to Landsat data. Correlation with sample cloud-free image reveals high performance even for haze or thin clouds.

CASE STUDY AND DATA
The methodology is tested on multispectral imagery from the Multi-Spectral Instrument (MSI) on-board the satellite Sentinel-2, which is a constellation of two twin satellites, Sentinel-2A and Sentinel-2B, operating since 2015 within the European Copernicus program for EO missions, coordinated by the European Space Agency (ESA). On the other hand, in order to emphasize potentials of optical satellite images for monitoring post-disaster events; we worked on the Lombok Island, in Indonesia, as case study, because on 5 th of August 2018, a destructive and shallow earthquake heavily affected the island.

Study Area
As a case study, we worked on Lombok, an island in Indonesia, Southeast Asia, situated at 8°33'54"S and 116°21'04"E, and part of the chain of the Lesser Sunda Islands. The island is roughly circular, with 70 kilometers across, covering an area of around 4,514 square kilometers. Census 2014 has recorded about 3.35 million of inhabitants. Between July and August 2018, the Island was hit by different destructive earthquakes. On August 5, 2018, it was the most destructive. The epicenter was inland, near Loloan village, in the north of the island. The earthquake also caused tsunamis. Several facilities was severely damaged or destroyed. Hundreds of people died, injured, or displaced. Here, we emphasize the effectiveness of clouds-and shadowuncontaminated data for monitoring post-disasters events. Therefore, we have considered the earthquake of the 5 th of August in Lombok for testing an approach capable of getting cloudless imagery for pre-and post-monitoring event analysis.

Sentinel-2 Multi-Spectral Instrument
Within the EO Copernicus program, Sentinel-2 is a wide-swath, high-resolution, multi-spectral imaging satellite equipped with an optical instrument named Multi-Spectral Instrument (MSI). The mission addresses several operational applications including soil, vegetation, and water monitoring, and it is consistent with other EO missions, namely SPOT and Landsat (ESA, n.d.). Sentinel-2 provides a high revisiting period, i.e. every 5 days under the same viewing angles. The MSI instrument is capable of systematically collecting optical data at mid-to high-spatial resolution in 13 spectral bands along the visible, near, and short-wave infrared part of the spectrum. Different spatial resolutions are provided, namely 10, 20, and 60 meters per pixel, as reported in Figure 1.  Then, all images have been processed in SNAP, using the Sen2Cor tool, which performs atmospheric-, terrain and cirrus correction, of TOA Level 1C data, to get Level-2A Bottom-Of-Atmosphere (BOA) reflectance processed images.
For this study, 24 images, for the same scene, have been downloaded and processed, of which 12 refer to the pre-disaster (before August 5, 2018), and 12 refer to the post-disaster situation (after August 5, 2018). Sentinel-2 image scenes are defined by granules, or tiles, which are 100x100 km 2 orthoimages in UTM/WGS84 projection (ESA, n.d.). The tile used for this experiment is identified by the TILE_ID: 50LMR, and EPSG: 32750.

METHODOLOGY
The M2CTool, as experimented on Landsat imagery (Colaninno et al., 2019), relies on different steps. After data pre-processing, an algorithm is designed for masking clouds and cloud shadows. Hence, masked images are combined to reconstruct a composite cloud-and shadow-free image scene based on multitemporal uncontaminated pixels. A multitemporal correction is applied to reduce seasonal variation along the stack.
When working with multitemporal data, even if for the same scene, and despite orthorectification procedure, a slight shift can occur among some images alternately. So, pixels may be not perfectly aligned. Such a concern is more relevant as the spatial resolution increases. In this study, we overlooked this problem and mainly focused on the testing the M2CTool methodology. At any rate, when the objective is of analyzing more images in a time series, although images refer to the same sensor and same spatial resolutions, images co-registration, based on selected Ground Control Points (GCP), should be performed (Gao, Zhang, & Gu, 2017;Scaioni, Barazzetti, & Gianinetto, 2018).

The Masking Algorithm
Recognized cloud masking algorithms rely on the use of thermal images for discriminating clouds. Instead, we tested an automatic algorithm for masking both clouds and cloud-shadows by means of spectral properties, for optical imagery. The approach is designed on the evidence that clouds are highly reflective at wavelength in the blue slice of the electromagnetic spectrum. Hence, for optical data, blue band is the best source for clouds detection. Similarly, near infrared (NIR) band is reported as the best source for detecting cloud shadows (Song & Civco, 2002). Consequently, as indicated by Lymburner (2010), if clouds are present in a stack of multitemporal images, blue band pixel with values smaller than the mean of all pixels at the same point along the stack, is plausibly clouds-free. While, if shadows are present in a multitemporal stack, NIR band pixel greater than the mean less the standard deviation of all pixels at the same point along the stack, is reasonably shadows-free (Colaninno et al., 2019;Lymburner, 2010). Based on such an assumption, an algorithm for automatically masking clouds and cloud-shadows has been set up. As outlined in Figure 4, and according to the previous statements, the algorithm relies on a statistical approach built upon the use of the Blue and NIR bands. For both bands, perpixel mean, and standard deviation are calculated, at each point of the image scene, through the stack of multitemporal images. Therefore, the mean image for the Blue and NIR stack, and standard deviation image for the NIR stack are obtained. Under the assumption that, along a multitemporal stack, clouds do not keep the same geographical position from one image to another, if clouds are present at any point, per-pixel mean, and standard deviation values increase significantly. Such an assumption allows identifying those pixels that are covered by clouds and cloud-shadows by comparing the covered pixels with same pixels at different time period. In other words, significant difference in terms of mean, and standard deviation, allow effectively identifying and masking clouds-and cloud-shadows contaminated pixels (Colaninno et al., 2019). The masking algorithm has been implemented as a script model in QGIS, setup on Sentinel-2 data, thus providing an automatic tool for masking clouds and cloud-shadow for optical imagery.

Reconstruction phase of the cloudless image
In this study we used 10 meters spatial resolution data, i.e. Blue, Green, Red, and NIR bands. Once masked all spectral bands, for each image scene, at the different time periods, a reconstruction phase is undertaken to achieve a cloud-free image for pre-and post-disaster situation. The reconstruction relies on two assumption (Colaninno et al., 2019), i.e. along the visible and the near infrared spectrum, high range values reasonably signify clouds, snow, thick haze, cirrus, and cumulus (Lymburner, 2010); while, because low range values signify high light absorption, it can be assumed that they can signify shadows. Based on this, for both pre-and post-disaster images, a stack of each band has been made and ordered depending on reflectance values, i.e., from the highest to the lowest. By ordering the images, we get a shaded stack ranging from the image of maximum to the image of minimum reflectance.
Because the maximum image along the stack can capture clouds, while the minimum image capture shadows, we have analyzed, by photointerpretation, all the images between maximum and minimum. We divided the ordered stack into percentiles, and selected the 30 th percentile as the most suitable, for this study, for reducing the possibility of keeping cloud contaminated and shadowed pixels. This allows discarding both high-ranges and low-ranges reflectance pixels, which have been not completely masked through the previous masking step.
However, because at this latitude images are strongly affected by cloud cover, the effect of masking caused areas of no-data to arise. On the other hand, thin clouds persist at the 30 th percentile. So, in order to recover both empty areas, as well as for replacing remaining thin clouds, the minimum image retrieved from the original unmasked images, is used to fill-in, and replace undesired pixels. The minimum image is completely, or almost completely free of clouds, although strongly affected by shadows, so it is only partially practical. Therefore, the reconstruction process, as depicted in Figure 6, aims at combining the minimum of the original unmasked images, with the 30 th percentile of the masked image stack per each band. Figure 6. Workflow of the M2CTool reconstruction process The reconstruction algorithm ( Figure 6) is organized as follow: the minimum of the original image stack is first obtained, then the 30 th percentile of the stack of the masked images is provided.
Because the 30 th percentile still shows small and thin clouds, the standard deviation image is calculated among the minimum of the original images, and the 30 th percentile. Hence, on a threshold value, pixels with high standard deviation are further masked, together with the no-data caused by the masking process. An enhanced masked 30 th percentile is obtained and overlapped to the minimum image, which is used to fill-in gaps. Then, the overlay is made by mosaicking the two images in ENVI, where the minimum image is put behind the 30 th percentile. A color balancing for reducing seasonal variation, and a feathering buffer is applied to obtain the most homogeneous result. Figure 7 and Figure 8 show the obtained multispectral cloudless image scenes for both pre-and post-disaster situation, respectively.

ANALYSIS OF RESULTS
In order to assess the effectiveness of the method, we consider either the capability of the masking algorithm, as well as the goodness of both pre-and post-disaster cloudless image scenes, obtained by applying the whole process. A statistical approach is employed for assessing the results.

The effectiveness of the masking algorithm
In order to evaluate the effectiveness of the employed approach for clouds/shadows masking, we provide a 20 square kilometers sample, of image scene taken at July 17, 2018 (Figure 9), where the original image is compared, by photointerpretation, with the mask (red color) achieved by applying the automatic algorithm. Besides, the effectiveness of the algorithm is revealed by the analysis of the histograms of the frequency distribution (Shen et al., 2015) at each band ( Figure 10). Either for high reflectance values (clouds), as well as for low values (shadows), which reasonably contaminate the scene, the tool shows good performance in terms of capability for masking clouds and cloudshadows. Indeed, the solid line, i.e. the masked image, shows a significant decrease of the frequency at both sides (low and high values), and at each wavelength, with respect to the original unmasked image. However, the effect is most visible at the NIR band. Although here we roughly estimate the effectiveness of the algorithm in terms of clouds/shadows masking, it has been already widely explored. Indeed, a deep analysis is provided by Colaninno et al. (2019), where the M2CTool is systematically compared with the FMask algorithm (Zhu et al., 2015). The FMask is a widely recognized algorithm, which shows the best overall accuracy among many other algorithms they have tested (Foga et al., 2017).

Assessment of the cloudless image scenes
In order to assess the performance of the approach for obtaining a clouds-and shadow-free imagery, pre-and post-disaster case studies have been tested. The objective is to validate the results with respect to a cloud-free (original) image. However, because image scenes completely uncontaminated are not available, samples of 20 km square have been considered along the less contaminated portions of the images selected for the test. A statistical approach is provided for estimating the goodness of the results. In particular, we analyzed scatterplots of reflectance values among bands, and the profile of reflectance along the blue band for emphasizing the effect of thin clouds and haze (Shen et al., 2015;Sun, Latifovic, & Pouliot, 2017).

Pre-disaster image scene analysis:
The image selected for assessing the result obtained for pre-disaster scene, is a sample of 20 km from the image July 02, 2018 ( Figure 11). The performance of the M2CTool has been assessed first by analyzing scatterplots of reflectance values. In particular, the reconstructed cloudless image has been compared (band-byband) with an uncontaminated squared sample of the Sentinel-2 image. Scatterplots for the four bands are given in Figure 12.  It has also to be taken into account that, optical imagery (mostly along the visible spectrum) are often contaminated by thin clouds and haze, which are not easy to identify and will always be challenging to mask (Foga et al., 2017). In this sense, because the reconstruction approach relies on low percentile values, it reduces the probability of getting thin clouds and haze. In order to assess this, and based on the assumption that at clear sky conditions Blue and Red bands are highly correlated for most land cover types (Shen et al., 2015;Sun et al., 2017;Tupas, 2015;Y. Zhang, Guindon, & Cihlar, 2002), we provide a regression analysis between the two bands, both for the reference image, as well as for the reconstructed image ( Figure 13). With respect to the reference image, the reconstructed image shows a slightly higher coefficient of determination. This accounts for an improved clear sky condition. In fact, if we analyze the spectral profile of reflectance at the Blue band (thin clouds and haze particularly affect the Blue spectrum), along an area supposedly covered by those atmospheric concerns, the reconstructed image shows lower values with respect to the reference Sentinel-2 image ( Figure 14).  In this case, correlation among bands is considerably higher. However, as stated for the previous case, low percentiles of the stack emphasize more the vegetation. Hence, highest and lowest correlation are provided, for Red and NIR band, respectively.
Thin clouds and haze contamination are also assessed for this case study, by means of a regression analysis among Blue and Red bands, for both reconstructed image and reference Sentinel-2 ( Figure 17). In this sense, the reconstruction approach further reduced the probability of getting thin clouds and haze, thus increasing clear sky conditions. In fact, correlation here provides a slightly higher R2 for the reconstructed image, in comparison with the previous case. Also, the spectral profile of reflectance at the Blue band, reveals that, even for this case, at certain points of the image some thin clouds or haze is present. Again, the reconstructed image shows lower values with respect to the reference image, along a transect supposedly affected by atmospheric concerns (Figure 18). Figure 18. Spectral profile of the Blue band reflectance along a transect for both reconstructed image and reference Sentinel-2

CONCLUSIONS
In this work we tested an accurate and cost/time-effective algorithm for masking and obtaining cloud-free optical satellite imagery, particularly useful for large and heterogeneous landscapes. The experiment is undertaken upon high resolution Sentinel-2 imagery, at a spatial resolution of 10 meters. Actually, because the methodology does not rely on the use of thermal data to mask clouds, it has been demonstrated the suitability for different optical satellite derived data (Colaninno et al., 2019).
However, we report that, depending on the availability of data, the methodology effectively allows obtaining effective images free, or almost completely free, of clouds/shadows. Indeed, the more the available images, the more is the capability of getting uncontaminated images. In fact, the availability of a stack of images for the same image scene, increases the possibility of selecting the most useful data. Here, we have experimented the possibility of working with percentiles along the stack, instead of quartiles as previously investigated. This increases the probability of further reducing thin clouds and haze. In this sense, we have demonstrated that, the effectiveness of the approach is obtained even when the images, supposedly uncovered by clouds, are affected by thin clouds and/or haze.
According to the analysis of the results provided at chapter 3, the approach appears quite effective in most of the areas along the whole image scene. However, we report that other sample areas have been tested, besides those presented here, and in some samples the results do not provide same improvements, with respect to the original uncontaminated image. We assume that it is reasonable because if an uncontaminated image already exists, there is not the need for producing further information.
At any rate, main encountered limits are mostly due to a higher difficulty to handle with shadows more than clouds, and because the algorithm is based on a statistical approach. In fact, as suggested by Lymburner (2010), a moving window approach should be explored in order to provide site-specific statistics.