Creating and Leveraging a Synthetic Dataset of Cloud Optical Thickness Measures for Cloud Detection in MSI

Cloud formations often obscure optical satellite-based monitoring of the Earth’s surface, thus limiting Earth observation (EO) activities such as land cover mapping, ocean color analysis, and cropland monitoring. The integration of machine learning (ML) methods within the remote sensing domain has significantly improved performance on a wide range of EO tasks, including cloud detection and filtering, but there is still much room for improvement. A key bottleneck is that ML methods typically depend on large amounts of annotated data for training, which is often difficult to come by in EO contexts. This is especially true when it comes to cloud optical thickness (COT) estimation. A reliable estimation of COT enables more fine-grained and application-dependent control compared to using pre-specified cloud categories, as is commonly done in practice. To alleviate the COT data scarcity problem, in this work we propose a novel synthetic dataset for COT estimation, that we subsequently leverage for obtaining reliable and versatile cloud masks on real data. In our dataset, top-of-atmosphere radiances have been simulated for 12 of the spectral bands of the Multispectral Imagery (MSI) sensor onboard Sentinel-2 platforms. These data points have been simulated under consideration of different cloud types, COTs, and ground surface and atmospheric profiles. Extensive experimentation of training several ML models to predict COT from the measured reflectivity of the spectral bands demonstrates the usefulness of our proposed dataset. In particular, by thresholding COT estimates from our ML models, we show on two satellite image datasets (one that is publicly available, and one which we have collected and annotated) that reliable cloud masks can be obtained. The synthetic data, the newly collected real dataset, code and models have been made publicly available at https: //github.com/aleksispi/ml-cloud-opt-thick .


Introduction
Space-borne Earth observation (EO) has vastly improved our options for collecting information from the world we live in.Not only atmospheric parameter from remote places, but also layers underneath the atmosphere are subject to these remote sensing applications.For example, land use and land cover classification [2], damage assessment for natural disasters [21,1], biophysical parameter-retrieval [5], urban growth monitoring [8], and crop yield estimation [31] are globally derived on a regular basis.Most of these applications rely on optical sensors of satellites gathering images, and the cloud coverage is a barrier that hampers the exploitation of the measured signal [9].Accurate cloud cover estimation is essential in these applications, as clouds can critically compromise their performance.
The cloud detection and estimation solutions in multispectral images range from rule-based statistical methods to advanced deep learning approaches.On the one hand, rule-based thresholding methods exploit the physical properties of clouds reflected on 1 arXiv:2311.14024v3[cs.CV] 15 Mar 2024 different spectral bands of satellite images to generate cloud masks.FMask [33] and Sen2Cor [20] are examples of such approaches implemented for cloud masking for Landsat and Sentinel-2 multispectral imagery, respectively.On the other hand, recently a burst of work in the literature has been using ML approaches to solve the cloud detection and estimation problem [22,32,15,14].These works use convolutional neural networks with large manually annotated datasets of satellite images and perform better than statistical rule-based methods.
Clouds are inhomogeneous by nature and the spatial inhomogeneity of clouds, or their varying cloud optical thickness (COT), affects not only the remote sensing imagery but also atmospheric radiation.The COT is an excellent proxy for a versatile cloud mask.By choosing a threshold value for the cloud/clear distinction (or even finer-grained categories), both a cloud-conservative and a clearconservative cloud mask can be implemented.Most literature estimates COT using independent pixel analysis (IPA, also known as ICA, independent column approximation [27]).IPA considers the cloud properties homogeneous in the pixel and carries no information about the neighboring pixels.Statistical approaches [19,34,13] have explored the potential causes influencing COT and defined parameters to mitigate their effect.With the advancement of deep learning [18,11], some approaches for COT using neural networks have been proposed [24,28].Most of these methods require neighboring information of the pixels and large amounts of annotated spatial data.
Common for statistical and ML-based methods is that both require benchmark datasets to evaluate their performance and find areas of improvement.Each pixel in the satellite images must be labeled manually or by a superior instrument (e.g. an active LIDAR [10]) as cloudy or cloud-free (or using even finer-grained labels).The complex nature of clouds makes the task even more difficult.This type of labeling is generally time-consuming and requires expert domain knowledge, leading to limited publicly available datasets.

Contributions
In the below, we provide details of the contributions of this work, which can be summarized as (i) the creation and public release of a synthetic 1D dataset simulating top-of-atmosphere (TOA) radiances for 12 of the spectral bands of the Multispectral Imagery (MSI) onboard the Sentinel-2A platform; and (ii) an extensive set of experiments with different ML approaches, which show that COT estimation can be used as a flexible proxy for obtaining reliable and versatile cloud masks on real data.
A new synthetic 1D dataset of TOA radiances for MSI.We propose a new 1D dataset simulating TOA radiances for 12 of the spectral bands1 of the MSI onboard the Sentinel-2A platform.The data points have been simulated taking into consideration different cloud types, cloud optical thicknesses (COTs), cloud geometrical thickness, cloud heights, water vapour contents, gaseous optical thicknesses, as well as ground surface and atmospheric profiles.For this, we connected the resources natively available in RTTOV v13 [26] with external resources, such as a dataset of atmospheric profiles provided by ECMWF, and a dataset of spectral reflectances from the ECOSTRESS spectral library [3,23].More details about the dataset are provided in §2.
The public release of this data implies that the threshold of research participation is lowered, in particular for non-domain experts who wish to contribute to this field, and offers reproducible and controllable benchmarking of various methods and approaches.Potential users of our dataset should be aware of the possible limitations the IPA could introduce in their study.In that regard, we note that if the users are interested in developing high spacial resolution applications, IPA will introduce a systematic error in COT, due to 3D effects [35].Nevertheless, the derivation of COT allows for more flexible use of data (e.g.clear-conservative vs cloud-conservative cloud-mask).
Extensive ML model experimentation based on synthetic and real data.We present results on a range of ML models trained for cloud detection and COT estimation using our proposed dataset.The experimental results suggest that on unseen data points, multi-layer perceptrons (MLPs) significantly outperform alternative approaches, like linear regression.We have also investigated and demonstrated the performance of our trained models on two real satellite image datasets -the publicly available KappaZeta, [7] as well as a novel dataset provided by the Swedish Forest Agency (which has also been made publicly available).These results show that by thresholding COT estimates from our ML models trained on the proposed synthetic dataset, reliable cloud masks can be obtained on real images; they even outperform the scene classification layer by ESA on the dataset provided by the Swedish Forest Agency.The main ML models are described in §3 and the experimental results are included in §4.

Synthetic Cloud Optical Thickness Dataset
In this section we explain the COT synthetic dataset that we propose to use for subsequent thresholdbased cloud masking.The dataset consists of 200,000 simulated data points (akin to individual and independent pixels, as observed by a satellite instrument), which have been simulated taking into consideration different cloud types, cloud optical thicknesses (COTs), cloud geometrical thickness, cloud heights, as well as ground surface and atmospheric profiles.Variations in water vapour content, gaseous optical thickness, and angles (azimuth difference as well as satellite and solar zenith angle) have also been included in the dataset.
To generate this data, the connection had to be established between surface and atmospheric properties on the one hand, and on the other hand, top-of-atmosphere (ToA) reflectivities, as observed by the MSI at Sentinel-2A.We did so using a radiative transfer model, namely the fast Radiative Transfer for TOVs (RTTOV) v13 [26].The simulations are performed for the instrument on the first of a series of Sentinel-2 platforms.Due to variations in between different realizations of the MSI and because of degradation of the instrument over time, it is recommended to add some noise to the reflection data before training to get more robust models (see also §3.1 and §4).
The backbone of each simulation is an atmospheric profile, defined by the vertical distribution of pressure, temperature, water vapor content and ozone.These atmospheric profiles are picked at random from an ECMWF (European Centre for Medium-range Weather Forecasts) dataset of 10,000 profiles compiled with the aim of mapping the widest possible range of atmospheric conditions [6].The surface is assumed to work as a lambertian reflector [17].All spectral reflectance are provided by the ECOSTRESS spectral library [3,23].Sun and satellite angles are varied randomly.
Each main surface type is composed of several sub-groups.The main surface types are nonphotosynthetic vegetation (e.g.tree bark, branches, flowers), rock (e.g.gneis, marble, limestone), vegetation (e.g.grass, shrub, tree), soil (e.g.calciorthid, quartzipsamment, cryoboroll) and water (ice, snow in different granularities, and liquid).The total number of different reflectivity profiles (including all not listed sub-groups) comprises 743 variations.These 743 fine classes appear in the dataset as 139 distinguishable categories.The distribution of the main types within the final dataset is shown in Table 1.The instrument-specific reflection was derived by convolving the spectral surface reflection with the spectral response function of the corresponding MSI channel.When generating the data, surface types were sampled at random.
The dataset consists of four equal-sized parts (50,000 points each).One part features clear sky calculations only (no clouds are considered here, only atmospheric gases and surface properties), while the three remaining parts respectively depict water clouds, ice clouds, and mixed clouds.Cloud geometrical thickness is set by random.The clouds geometrical thickness varies from 1 layer to 6 with the mode at 2 and 3. Clouds of the same type populate adjacent layers but vertical gaps are possible between ice and water clouds.The cloud type (for water clouds) and liquid water/ice content is varied randomly.Water clouds are parameterized following the size distributions of the five OPAC cloud types [12] with updates in the liquid water refractive index dataset according to [29].Ice clouds are parameterized by temperature and ice water content (IWC) according to [4].Variations in size and shape distributions are taken into account implicitly by the parametrization of particle optical properties.Distributions of COT for water, ice and combined clouds are shown in Fig. 1, where it can be seen that our focus is set on optically thin clouds.Note that when a consistency-check fails (i.e. the cloud type and background atmospheric temperature do not match), the cloud will be removed and the COT is set to 0, but the data point remains in its place.Thus, within the dataset there are some data points that are labeled cloudy even though they depict clear sky cases.These clear sky cases are removed from the plots in Fig. 1.
The overall cap of COT ≤ 50 is a feature of the radiative transfer model.This is not seen as a limitation for model training, because COT is mainly used as a means for subsequent cloud classification.For example, see §4.2 where we successfully perform classification of real image pixels into semi-transparent and opaque (optically thin and thick) clouds.We do this even though we understand that the pure prediction of the optical density, due to 3D effects, has little significance.The most relevant 3D effect for MSI with a high spatial resolution is the neglect of the net horizontal radiation transport.The error depends on the horizontal component of the photon movement and thus largely on the solar zenith angle.However, even the incorrect COT is suitable for a cloudy/cloud-free decision.Note that according to the international satellite cloud climatology project (ISCCP, see e.g.[25]), clouds are considered optically thin for COT < 3.6, medium for 3.6 ≤ COT < 23 and optically thick for COT ≥ 23.Keeping in mind that optically thick clouds are in general easier to detect, it is not expected that this limitation affects our conclusions.
In summary, our dataset D = {(x 1 , y 1 ), . . ., x n , y n )} consists of 200,000 pairs (x i , y i ) with x i ∈ R 19 and y i ∈ R + the i:th input and ground truth COT, respectively, i.e. each element in the 19-dimensional vector x i is a real scalar, and y i is a non-negative real scalar (since COT cannot be negative).By ground truth COT y i , we mean the true COT against which COT predictions can be compared (e.g.measure how accurately a machine learning model predicts y i if it gets x i as input).As for x i , it contains simulated measured reflectivities for the 12 spectral bands, but also additional optional features: (i) satellite zenith angle In this work we have focused on experiments using only the 12 band reflectivities as model input, which are least effected by aerosols, and thus from this point onwards we will let x i ∈ R 12 denote such a data point.One of the main reasons for not using other types of input information besides the band reflectivities is compatibility (to ensure that the resulting models would be able to be evaluated on publicly available real image datasets, see e.g.§4.2).As for the ground truth COT y i , these lie in the range y i ∈ [0, 50]; again, the upper limit is set by the radiative transfer model.We randomly set aside  Note that we use an MLP, rather than e.g. a convolutional architecture, since the data points x i are independent from each other in our synthetic dataset.160,000 data points for model training, 20,000 for validation, and 20,000 for testing (see §3.1 and §4).

Machine Learning Models
In this section we briefly introduce the different machine learning (ML) models that we have implemented for the COT estimation task.Recall from §2 that the model inputs x i ∈ R 12 contain only 12 band reflectivities, i.e. no auxiliary inputs such as surface profile information is used by the models.
As our proposed dataset can be seen as consisting of individual 'pixels', with no spatial relationship in-between them (cf.§2), we have mainly considered multi-layer perceptron (MLP) models. 2Extensive model validations suggested that a five-layer MLP with ReLU-activations and hidden dimension 64 (same for all layers) yields the best results.See Fig. 2 for an overview of this architecture.
When performing inference on real imagery I ∈ R H×W ×C with an MLP model, we let it operate at a pixel-per-pixel basis, as the models are trained (see §3.1) on individual synthetic pixels due to the nature of the synthetic data.To improve spatial consistency of a resulting COT prediction map C ∈ R H×W , we apply a simple post-processing scheme where a sliding window of size M ×M and stride 1 runs on top of C ∈ R H×W , and produces an average (among M 2 values) at each location.Overlaps that occur due to the stride being smaller than the width of the kernel are handled by averaging.We found M = 2 to produce robust results.Thus, via this post-processing each value within C is essentially replaced by the mean value of the surrounding values, which increases the smoothness of the predictions (nearby predictions are less likely to differ drastically).This is desirable, because in real images adjacent pixels are typically quite similar to each other, so the predictions for adjacent pixels are expected to be quite similar as well.

Model Training
As is standard within machine learning, we normalize the model inputs x i ∈ R 12 such that each of the 12 entries of x i has zero mean and unit standard deviation.This is done by first computing the mean x mean ∈ R 12 and standard deviation x std ∈ R 12 on the training set, and then for each x i applying the following transformation (the subtraction and division below is done element-wise for each of the 12 entries) before feeding the result to the model: We also explored other data normalization techniques but did not see any improvements relative to (1).To improve model robustness, independent Gaussian noise is also added to each input.This noise is zero-mean and has a standard deviation corresponding to 3% of the average magnitude of each input feature.Thus, the transformation (1) is replaced by where ε i ∈ R 12 denotes such an independently sampled Gaussian.
The ML models are trained with the commonly used Adam optimizer [16] with a mean squared error (MSE) loss for 2,000,000 batch updates (batch size 32, learning rate 0.0003).Since the models are very lightweight, training can be done even without a graphics processing unit (GPU) within about an hour.

Fine-tuning with Weaker Labels
For some datasets (see §4.2) we have access to weaker labels in the form of pixel-level cloud masks.For instance, in the KappaZeta dataset [7], there are labels for 'clear', 'opaque cloud', and 'semitransparent cloud', respectively.In addition to evaluating models on such data, we can also set aside a subset of that data to perform model refinement based on the weaker labels.
Let τ semi and τ opaque (where 0 < τ semi < τ opaque ) denote the semi-transparent and opaque COT thresholds, respectively.Then, we refine a model using a loss L which satisfies the following criteria.If p denotes the prediction of a pixel which is labeled as • 'opaque cloud', then • 'semi-transparent cloud', then

Experimental Results
In this section we extensively evaluate various ML models discussed in §3.For COT estimation, we evaluate the following models: • MLP-k is a k-layer MLP (the main model has k = 5, cf.Fig. 2), where results are averaged among 10 identical k-layer MLP (only differing in the random seed for the network initialization); • MLP-k-ens-n is an ensemble of n k-layer MLPs, each identically trained but with a unique random network initialization; • MLP-k-no-noise is the same as MLP-k, but is trained without adding any noise on the training data (recall that we train with 3% additive noise by default); • lin-reg is a linear regression model.

Results on Synthetic Data
For our proposed synthetic dataset (cf.§2) we have access to ground truth COTs, and thus in Table 2 we report the mean absolute error (MAE) between the predicted outputs and corresponding ground truths on unseen test data.We see that training on data with artificially added noise improves model robustness significantly, whereas ensembling only marginally improves performance.Also, MLPs significantly outperform linear regression models.
As mentioned in §2, our synthetic dataset covers five main surface profiles: vegetation, rock, nonphotosynthetic vegetation, water and soil.Also recall that vegetation is by far the most common main profile in the dataset (70.5% of the data), whereas rock is the second most common one (10.7% of the data).To investigate how the prediction accuracy of the best model MLP-5-ens-10 performs on these five different main profiles, we also compute five main profile-conditioned test results (averaged across the different noise percentages, cf.bottom

Results on Real KappaZeta Data
This publicly available dataset of real satellite imagery [7] is used to assess the real-data performance of models that were trained on our synthetic dataset.Since ground truth information about COT is not available, we consider the weaker pixel-level categories 'clear', 'semi-transparent cloud' and 'opaque cloud'.We set aside a subset (April, May, June; 3543 images that we refer to as the training set) on which we search for COT thresholds which correspond to semi-transparent and opaque cloud, respectively.These thresholds were respectively found to be 0.75 and 1.25 (both values in the range [0, 50], cf.§2).
Given these thresholds, we evaluate a model's semantic segmentation performance (into the three aforementioned cloud categories) on July, August, September (2511 images in total, which we refer to as the test set; there is no spatial overlap between the training and evaluation sets).We also conduct similar experiments where we first refine the MLP models on the training set of KappaZeta, cf.§3.2.Such a model is denoted MLP-k-tune.This comparison allows us to quantify how much is gained in prediction accuracy if one trains on the target data domain (the training set of KappaZeta), as opposed to training only on off-domain data (our synthetic dataset).We also compare our MLP models with a U-net segmentation network3 for the same task.Different to an MLP, a U-net is a fully convolutional architecture that explicitly incorporates contextual information from surrounding pixels in an image when performing predictions.Thus this comparison with a U-net allows us to quantify how much is gained in prediction accuracy if one uses a model that takes into account spatial connectivity among pixels -we expect the U-net to be superior, since U-nets incorporate and take advantage of spatial patterns in images, while the MLPs inherently cannot.
The results are shown in Table 3, where two metrics are reported for each cloud category: • F1-score, i.e. the harmonic mean of precision and recall, where precision is the fraction of pixels that were correctly predicted as the given category, and recall is the fraction of pixels of the given category that were predicted as that category; • intersection-over-union (IoU), i.e. the ratio between the number of pixels in the intersection and union, respectively, of the predicted and ground truth cloud masks.
Both metrics are bounded in the [0, 1]-range (higher is better).For both metrics we also show the average result over the different cloud types (these are respectively denoted F1-avg and mIoU), where the average is computed uniformly across the different cloud types, i.e. each category contributes equally Table 3: Results on the KappaZeta test set.MLP approaches do not perform as well as U-nets, as is expected since U-nets integrate information from other pixels, which MLPs inherently cannot.See also qualitative results in Fig. 3.Among the MLP approaches, model fine-tuning on the KappaZeta training set only slightly improves results (column 4 vs 3), i.e. our synthetic dataset allows for model generalization to real data, while ensembling tuned models negligibly improves results (column 4 vs 2).   to the mean, independently of how common the category is in the dataset.
We see that the best MLP-approach (MLP-5ens-10-tune) obtains a lower mIoU and F1-score compared to the U-net (0.47 vs 0.54 and 0.52 vs 0.66, respectively).We re-emphasize however that the MLP approaches perform independent predictions per pixel, whereas the U-net takes into consideration spatial context.Furthermore, we see that model fine-tuning on the KappaZeta training set only slightly improves results on the test set (column 4 vs column 3), which shows that our In each example, columns 1 and 4 are the same as in Fig. 3, while column 2 shows pixel-level cloud type predictions based on thresholding the COTs of an MLP ensemble approach that was trained only on our synthetic data.Column 3 is the same as column 2, but the MLPs were refined on KappaZeta training data.Fine-tuning sometimes yields better results (e.g.top two rows on the left).In many cases, the results are however very similar before and after fine-tuning (e.g.third row on the left and right side), and sometimes results get worse after fine-tuning (e.g.fourth row on the left).synthetic dataset can be used to train models which generalize well to cloud detection on real data.It can also be seen that model ensembling only marginally improves results for KappaZeta-tuned models (column 4 vs column 2).

Qualitative results.
We further examine the results by inspecting qualitative examples in Fig. 3.The examples on the left side focus on our MLP approach (the one that was only trained on synthetic data) and includes the estimated COT values (column 2).On the right side we include comparisons to the U-net.From these examples we see that the ground truth segmentation masks are 'biased' towards spatial connectivity among pixels -in many cases, the ground truth does not seem entirely accurate, and rather places emphasis on spatial consistency at the cost of finer-grained contours.This bias is easily incorporated in an architecture which can leverage spatial connectivity (the U-net in our case), whereas it will cause confusion for per-pixel models (the MLPs in our case).In many of the examples it appears as if the MLP predictions are more accurate than the ground truth, when visually comparing with the associated input images (the main exception is shown at row #4 (right) where the MLP approach underestimates the COTs).
In Fig. 4 we show qualitative examples of the effect of model refinement on the KappaZeta training set.As can be seen, in many cases, results do not improve by model fine-tuning, i.e. our synthetic dataset can on its own provide models that perform well on many real examples.Recall that this can also be seen quantitatively by comparing the third and fourth columns of Table 3, where one sees that the KappaZeta-tuned variant obtains only slightly higher F1 and mIoU scores.
Model runtimes.The CPU runtime of an Table 4: Results on test data from the Swedish Forest Agency (SFA).The various MLP approaches were not trained on a single SFA data point.Our MLP-5-ens-10 model yields results comparable to ResNet18-cls (which was trained on the SFA training data), and it outperforms the SCL by ESA.Note that MLP-5 represents the average result among ten MLP-5 models that were trained using different network parameter initializations, whereas MLP-5-ens-10 is an ensemble of those ten different models.

Results on Real Data from the Swedish Forest Agency
The Swedish Forest Agency (SFA) is a national authority in charge of forest-related issues in Sweden.Their main function is to promote management of Sweden's forests that enable the objectives of forest policies to be attained.Among other things, the SFA runs change detection algorithms on satellite imagery, e.g. to detect if forest areas have been felled.For this, they rely on ESA's scene classification layer (SCL), which also includes a cloud probability product.The SFA's analyses require cloud-free images, but the SCL layer is not always accurate enough.Therefore, we applied models that were trained on our synthetic dataset on the SFA's data in order to classify a set of images as 'cloudy' or 'clear'.
To achieve this, the SFA provided 432 Sentinel-2 Level 2A images4 of size 20 × 20 (corresponds to 200 × 200 m 2 ) that they had labeled as cloudy or clear (120 cloudy, 312 clear), where an image was labeled as clear if no pixel was deemed to be cloudy.Fig. 5 shows the locations of these images within Sweden.The 432 images were randomly split into a training, validation and test split, such that the respective splits have the same ratio between cloudy / clear images (i.e.roughly 28% cloudy and 72% clear images per split).The training, validation and test sets contain 260, 72 and 100 images, respectively.
The results on the test set are shown in Table 4, which also includes the results of using the ESA-SCL.For our MLPs, we use the validation data in order to set a COT threshold above which a pixel is predicted as cloudy (we found 0.5 to be the best threshold).For the SCL, a pixel is predicted to be cloudy if the SCL label is 'cloud medium probability', 'cloud high probability', or 'thin cirrus'.For both our MLPs and the SCL approach, if a pixel is predicted as cloudy, the overall image is predicted as cloudy.We also compare with a ResNet-18 model [11] for binary image classification (cloudy and clear, respectively), trained on the union of the training and validation sets.Left-right and bottom-up flipping was used as data augmentation.
From Table 4 we see that the our main model (an ensemble of ten 5-layer MLPs) is on par with the dedicated classification model, despite not being trained on a single SFA data point (except for COT threshold tuning), and despite that the training data represents top-of-atmosphere data (i.e.Level 1C, not Level 2A as the SFA data).We also see that model ensembling is crucial (MLP-5-ens-10 vs MLP-5), and that our MLP-5-ens-10 model significantly outperforms the SCL (average F1-score 0.88 vs 0.68 for the SCL).Finally, note that even using only a single MLP-5 model yields a higher average F1-score than the SCL.

Conclusions
In this work we have introduced a novel synthetic IPA dataset that can be used to train models for predicting cloud types of pixels (e.g.clear, semi-transparent, opaque clouds).In order to enable broad application, in our dataset we chose to use cloud optical thickness (COT) as a proxy for cloud masking.Several ML approaches were explored on this data, and it was found that ensembles of MLPs perform best.Despite the fact that our proposed synthetic dataset (and thus associated models) is inherently pixel-independent, the models show promising results on real satellite imagery.In particular, our MLP approach trained on pixel-independent synthetic data can seamlessly transition to real datasets without requiring additional training.To showcase this, we directly applied (without further training) this MLP approach on the task of cloud classification on a novel real image dataset, and achieved an F1-score that is on par with a widely used deep learning model that was explicitly trained on this data.Furthermore, our approach is superior to the ESA scene classification layer at classifying satellite imagery as clear or cloudy, and can flexibly generate cloud type segmentation masks via COT thresholding.Code, models and our proposed datasets have be made publicly available at https: //github.com/aleksispi/ml-cloud-opt-thick.

Figure 1 :
Figure 1: Distribution of COTs in the range [0, 50] for different cloud types in our dataset: water clouds (left), ice clouds (center), and combined water and ice clouds (right).The focus of our data is on optically thin clouds.

Figure 2 :
Figure 2: Main MLP model architecture that we use for COT estimation.The input layer is shown on the left (4 nodes shown for the input x instead of all 12, to avoid visual clutter).Then follows four hidden layers (8 nodes shown for each, instead of 64).Finally, the output layer on the right produces the COT estimate ŷ.Each of the five layers is of the form z k = ReLU(W k z k−1 + b k ) = max(0, W k z k−1 + b k ) for k = 1, . . ., 5, with z 0 = x and z 5 = ŷ.Here, W k and b k respectively represent the weight matrix and bias vector for the k:th layer, and these parameters are calibrated through the training process (see §3.1).Note that we use an MLP, rather than e.g. a convolutional architecture, since the data points x i are independent from each other in our synthetic dataset.
Figure 2: Main MLP model architecture that we use for COT estimation.The input layer is shown on the left (4 nodes shown for the input x instead of all 12, to avoid visual clutter).Then follows four hidden layers (8 nodes shown for each, instead of 64).Finally, the output layer on the right produces the COT estimate ŷ.Each of the five layers is of the form z k = ReLU(W k z k−1 + b k ) = max(0, W k z k−1 + b k ) for k = 1, . . ., 5, with z 0 = x and z 5 = ŷ.Here, W k and b k respectively represent the weight matrix and bias vector for the k:th layer, and these parameters are calibrated through the training process (see §3.1).Note that we use an MLP, rather than e.g. a convolutional architecture, since the data points x i are independent from each other in our synthetic dataset.

Figure 3 :
Figure 3: Left: Examples of our main MLP approach (ensemble of ten 5-layer MLPs) on unseen KappaZeta test data.Column 1: Input image (only RGB is shown).Column 2: COT estimates (relative intensity scaling, to more clearly see variations).Column 3: Pixel-level cloud type predictions based on thresholding the COTs in column 2. Column 4: KappaZeta ground truth.Dark blue is clear sky, lighter blue is semi-transparent cloud, and turquoise is opaque cloud.Right: Similar to the left, but the 2nd column shows the thresholded model predictions (instead of the COT estimates), and the 3rd column is the U-net prediction.A failure case for both models is shown on the bottom row.

Figure 4 :
Figure 4: Eight additional qualitative examples on the KappaZeta test set (examples are shown in the same format on the left and right side of this figure).In each example, columns 1 and 4 are the same as in Fig.3, while column 2 shows pixel-level cloud type predictions based on thresholding the COTs of an MLP ensemble approach that was trained only on our synthetic data.Column 3 is the same as column 2, but the MLPs were refined on KappaZeta training data.Fine-tuning sometimes yields better results (e.g.top two rows on the left).In many cases, the results are however very similar before and after fine-tuning (e.g.third row on the left and right side), and sometimes results get worse after fine-tuning (e.g.fourth row on the left).

Figure 5 :
Figure 5: Locations in Sweden of the annotated imagery provided by the SFA.Larger dots indicate a higher density of images in the associated region.

Table 1 :
Distribution of the main surface types within our dataset.Details and examples of subgroups are given in the main text.

Table 2 :
Mean absolute errors (MAEs) on different variants of the test set of our synthetic dataset.Test-x% refers to the test set with x% added noise.Ensemble methods marginally improve over single-model ones.Models trained with additive input noise yield significantly better average results.Linear regression performs worst by far.Note that for the single-model variants MLP-5, MLP-5-no-noise and Lin-reg, we show the mean MAE over 10 different network parameter initializations, and standard deviations.