Space-scale exploration of the poor reliability of deep learning models: the case of the remote sensing of rooftop photovoltaic systems

Abstract

Photovoltaic (PV) energy grows rapidly and is crucial for the decarbonization of electric systems. However, centralized registries recording the technical characteristics of rooftop PV systems are often missing, making it difficult to monitor this growth accurately. The lack of monitoring could threaten the integration of PV energy into the grid. To avoid this situation, remote sensing of rooftop PV systems using deep learning has emerged as a promising solution. However, existing techniques are not reliable enough to be used by public authorities or transmission system operators (TSOs) to construct up-to-date statistics on the rooftop PV fleet. The lack of reliability comes from deep learning models being sensitive to distribution shifts. This work comprehensively evaluates distribution shifts’ effects on the classification accuracy of deep learning models trained to detect rooftop PV panels on overhead imagery. We construct a benchmark to isolate the sources of distribution shifts and introduce a novel methodology that leverages explainable artificial intelligence (XAI) and decomposition of the input image and model’s decision regarding scales to understand how distribution shifts affect deep learning models. Finally, based on our analysis, we introduce a data augmentation technique designed to improve the robustness of deep learning classifiers under varying acquisition conditions. Our proposed approach outperforms competing methods and can close the gap with more demanding unsupervised domain adaptation methods. We discuss practical recommendations for mapping PV systems using overhead imagery and deep learning models.

Impact Statement

This paper analyzes the effects of distribution shifts on deep learning models trained to detect rooftop photovoltaic (PV) systems on aerial imagery by combining explainable artificial intelligence methods. It then proposes practical solutions based on this analysis to enhance the robustness of these models, thereby improving their reliability and facilitating the use of remote sensing techniques to support the insertion of rooftop PV systems into the grid. The methodology laid out in this work can be replicated for other case studies.

1. Introduction

Photovoltaic (PV) energy grows rapidly and is crucial for the decarbonization of electric systems (Haegel et al., 2017). The rapid growth of rooftop PV systems makes it challenging to of estimate the global PV installed capacity, as centralized data is often lacking (Hu et al., 2022; Kasmi et al., 2022). Remote sensing of rooftop PV systems using orthoimagery and deep learning models is a blooming solution for mapping rooftop PV installations. Deep learning-based pipelines have become the standard method for remote sensing PV systems, with works like DeepSolar (Yu et al., 2018) paving the way for country-wide mapping of PV systems using deep learning and airborne or spaceborne orthoimagery. Recently, methods for mapping rooftop PV systems in many regions, especially in Europe, have been proposed (Frimane et al., 2023; Kasmi et al., 2022; Kausika et al., 2021; Lindahl et al., 2023; Mayer et al., 2020; Rausch et al., 2020; Zech and Ranalli, 2020). Some of these works (Kasmi et al., 2022; Mayer et al., 2022) introduced methods to estimate the technical characteristics of the PV systems (individual localization, orientation, PV installed capacity). The identification of rooftop PV systems improves their integration into the electric grid by improving the ability of transmission system operators (TSOs) to more accurately estimate their power production in real-time (Kasmi et al., 2024) but can also promote their future expansion, as this data helps understanding the drivers behind rooftop PV adoption (Alipour et al., 2020; Colas and Saulnier, 2024; Graziano and Gillingham, 2015; Wang et al., 2022).

However, deep learning-based detection methods are sensitive to so-called distribution shifts, i.e., differences between the training and testing data (Koh et al., 2021). This sensitivity manifests by unpredictable and sharp accuracy drops when the model is deployed on unseen images. It limits their practical usability as a trained model cannot be deployed without retraining to carry out registry updates. Besides, the unpredictability of the model’s behavior limits its reliability as it casts doubt on what it perceives as a PV panel (De Jong et al., 2020; Hu et al., 2022). In this work, we define the reliability of a model as its ability to rely on relevant features to identify PV systems, i.e., to be “right for the right reasons” and to be simultaneously able to rely on robust features (Kasmi, 2024; Ross et al., 2017). Steps towards improving the quality of registries (i.e., tables recording the location and some technical information on PV systems) of rooftop PV systems constructed using deep learning algorithms have been taken, with Hu et al., 2022 and Kasmi et al., 2022 discussing the practical evaluation of the mapping algorithms or Li et al., 2021 identifying the minimum resolution to detect rooftop PV systems from orthoimagery (whether this minimum resolution is the native image resolution or the resolution obtained after increasing the input image resolution with methods such as those proposed by Ho et al., 2022). To date, Wang et al. (2017) is the only work that studied the poor generalizability of PV mapping algorithms, though it was limited to two cities and one image dataset. More recently, Pena Pereira et al. (2024) analyzed how PV system typologies and backgrounds affect performance, recommending input patch size adjustments and data augmentation to improve detection accuracy. Despite these advances, further work is needed to understand what the model identifies as a PV panel during training and how distribution shifts—arising from variations in PV systems, backgrounds, or acquisition conditions—impact performance.

This work aims to improve the reliability of deep learning models deployed in real-world settings prone to distribution shifts, taking the remote sensing of rooftop PV systems as a case study. We introduce a novel methodology to understand and address the sensitivity to distribution shifts based on empirical experiments and throughout analysis of the model’s decision using explainable AI (XAI) methods. Empirical evaluation and XAI methods enable us to identify the most important sources of distribution shifts and grasp why they occur. We evaluate a wide range of popular domain adaptation techniques (i.e., methods that aim at reducing the sensitivity to distribution shifts of deep learning algorithms) and introduce a novel data augmentation method. This method, based on our empirical findings, aims at effectively and reliably reducing the sensitivity to distribution shifts of deep learning models trained to detect PV systems from orthoimagery. We discuss practical takeaways regarding the choice of training data and domain adaptation methods for remote sensing PV systems. Since the sensitivity to distribution shifts is a recurring issue with the real-world deployment of deep learning systems (Koh et al., 2021), we discuss the requirements for applying our methodology to alternative use cases.

The code for replicating the results of this paper can be found at https://github.com/gabrielkasmi/robust_pv_mapping, and model weights can be found at https://zenodo.org/records/14673918.

2. Related works

2.1. Remote sensing of rooftop photovoltaic installations

The remote sensing of rooftop PV systems is now a well-established field with early works dating back to Golovko et al., 2018; Malof et al., 2015, Malof et al., 2016; Yuan et al., 2016. The DeepSolar project (Yu et al., 2018) marked a significant milestone by mapping distributed and utility-scale installations over the continental United States using state-of-the-art deep learning models. Many works built on DeepSolar to map regions or countries, especially in Europe, covering areas such as North-Rhine Westphalia (Mayer et al., 2020), Switzerland (Casanova et al., 2021), Oldenburg in Germany (Zech and Ranalli, 2020), parts of Sweden (Frimane et al., 2023; Lindahl et al., 2023), Northern Italy (Arnaudo et al., 2023), the Netherlands (Kausika et al., 2021) or the surroundings of Berkeley in California (Parhar et al., 2021), Connecticut (Malof et al., 2019) or the surroundings of Sfax, in Tunisia (Bouaziz et al., 2024). Several works even included GIS data to construct registries of PV installations (Kasmi et al., 2022; Kausika et al., 2021; Mayer et al., 2022; Rausch et al., 2020). In the current context of rapid rooftop PV growth (Haegel et al., 2017; RTE France, 2022), remote sensing of rooftop PV installations using deep learning and orthoimagery offers the potential to address the lack of systematic registration of small-scale PV installations (Kasmi et al., 2022; Kausika, 2022).

However, current methods cannot be transposed from one region to another without incurring accuracy drops, thus limiting their practical usability (Hu et al., 2022), as the aim of these models is to be regularly deployed on new images to construct and maintain official registries of PV systems (De Jong et al., 2020). The unpredictability of the accuracy drops also casts doubt regarding the reliability of these methods in such applied settings. To address this issue, Kasmi et al., 2022 recently introduced a method aiming at indirectly assessing the accuracy of the detections by automatically comparing the registry generated by deep learning algorithms to reference data, which is often aggregated at the city scale. While this work enabled the quantification of the drop in accuracy encountered during deployment, no cues as to why the accuracy varied during deployment were discussed. Kasmi et al. (2023a) introduced a benchmark to disentangle the sources of distribution shifts occurring with orthoimages of PV systems and outlined some promising directions to improve the reliability of deep learning algorithms. This work builds on and deepens the analysis of Kasmi et al., 2023a to propose a methodology for identifying the main sources of distribution shifts when dealing with the remote sensing of rooftop PV systems, understanding how these shifts affect deep learning models and extensively discussing how explainable AI techniques and domain adaptation methodologies can help mitigate the sensitivity to distribution shifts while improving the end user’s trust towards deep learning black-boxes.

2.2. Distribution shifts and domain adaptation

Definition. Distribution shifts, i.e., the sensitivity to the fact that “the training distribution differs from the test distribution” (Koh et al., 2021) are ubiquitous in machine learning (Torralba and Efros, 2011). The sensitivity to distribution shifts causes unpredictable performance drops, which can have dire consequences as models are deployed in safety-critical settings such as autonomous driving (Sun et al., 2022b), medical diagnoses (Pooch et al., 2020) or finance (Thimonier et al., 2024). Distribution shifts formally consist of a break in the assumption that the training and testing (or deployment) data are independently and identically distributed (i.i.d., Zhou et al., 2023). This assumption is central when training models to minimize the empirical risk, as the underlying assumption is that the empirical risk is a good approximation of the true risk of the estimator. This assumption is true only if the data is i.i.d.; otherwise, the empirical risk no longer represents the true risk. It corresponds to epistemic uncertainty (Gal, 2016) as the model is exposed to data outside its prior training experience.

Distribution shifts in remote sensing. Due to its nature, remote sensing data often breaks the i.i.d. assumption (Tuia et al., 2024). For instance, the raw imagery consists of large image tiles cut into smaller thumbnails before being passed to the model. Therefore, the thumbnails exhibit a strong spatial correlation. Tuia et al. (2016) identified two primary sources of shifts in the input data to which models are sensitive: variations in the geographical scenery and varying acquisition conditions. Following Murray et al., 2019, we can add a third one: the ground sampling distance (GSD).

The acquisition conditions encompass the conversion of a scene into a digital image and include all sources of variability in the input images caused by different sensors, exposure, attitude and altitude during acquisition, and atmospheric conditions. The ground sampling distance (GSD) is the upper bound to the image’s effective resolution. The effective resolution considers the distortions induced by the angle of incidence of the sensor (e.g., RGB camera). The lower the ground sampling distance, the more detailed the image. In practice, the effective resolution is limited by the GSD and the image quality (noise, optical transfer function, and intrinsic geometric consistency). In this article, with a slight abuse of wording, we will use the terms “GSD” and “resolution” interchangeably. The GSD corresponds to the distance between two consecutive pixels measured on the ground and is expressed in meters per pixel. The resolution measures the number of pixels per unit of length (e.g., inches or centimeters), and the image size describes its dimension. For consistency with the related literature, we will use the term “resolution” when broadly referring to the GSD. However, we will explicitly use “GSD” in specific contexts, such as when expressing it with its unit, as it makes no sense to refer to a “resolution of 0.2 m/pixel”.

So far, the only work that investigated the poor reliability of deep learning models applied to the remote sensing of PV panels is Wang et al., 2017. The authors argued that the generalization ability from one city to another depends on how “hard” to recognize the PV panels are. However, no proper definition of the “hardness” to recognize PV panels or a proper disentanglement of the effect of each source of variability was carried out, and there was no prescription regarding model training or data preprocessing. More recently, Li et al., 2021 and Pena Pereira et al., 2024 studied the practical implications of having different resolutions or PV panel instances on the model’s performance. These studies focus on the observable impacts of factors such as system heterogeneity, ground sampling distance, or image resolution on performance but overlook the underlying mechanisms driving the performance degradation. Following Kasmi et al., 2023a, we consider that improving the reliability of PV mapping models requires a deeper understanding of the underlying reasons for their sensitivity to these distribution shifts. Finally, to the best of our knowledge, Kasmi et al., 2023a is the only work to have implemented some domain adaptation techniques in the context of mapping rooftop PV systems.

Distribution shifts and domain adaptation. Domain adaptation is the go-for approach to address the sensitivity of machine learning models to distribution shifts (Ben-David et al., 2006). The different distributions are referred to as the source domain S, on which the model is initially trained, and the target domain T, on which the model is deployed. Different approaches can be distinguished depending on the number of source and target domains or the availability of labeled data. In its most constrained setting, we assume that we have access to labeled data from the source domain and, at most, unlabeled data from the target domain. This setting is sometimes referred to as unsupervised domain adaptation. We refer the reader to surveys such as Csurka, 2017; Csurka et al., 2021; Guan and Liu, 2022; Tuia et al., 2016; Zhou et al., 2023 for an extensive discussion of the domain adaptation settings and techniques. The general idea of domain adaptation is to learn a representation of the data invariant across domains or, equivalently, insensitive to distribution shifts.

We can distinguish two broad approaches to domain adaptation: implicit and explicit regularization. Implicit regularization encourages the model to generalize across domains without imposing specific constraints on the loss function during the initial training. Data augmentations form a first class of implicit regularization techniques. By viewing multiple copies of the same image that have been altered, the model learns to be invariant to these alterations. The aim is that the model is no longer sensitive to a given set of perturbations of the input images. Popular data augmentation methods consist of defining a method to generate as many perturbed samples as possible while preserving the semantic content of the image. To this end, AugMix (Hendrycks et al., 2020) applies a random sequence with random weights of perturbations to the input image. Similarly, Hendrycks et al., 2022 augmented an input image with fractal patterns, and Sun et al., 2022a perturbed the Fourier spectrum of the input image. Cubuk et al., 2019 used a reinforcement learning framework to search for an optimal augmentation policy, selecting the type, magnitude, and probability of transformations based on a target validation set, and Cubuk et al., 2020 simplified this framework to make it computationally less demanding. Another approach for implicit regularization is to modify the model’s architecture by enforcing additional invariances, such as the invariance to various groups of translations, reflections, and rotations as done by Cohen and Welling, 2016.

On the other hand, explicit regularization techniques require access to several source domains or unlabeled samples from the target domain, making these approaches more demanding than the implicit ones. The most popular approach is CORrelation ALignment (CORAL), and its counterpart for deep models DeepCORAL (Sun and Saenko, 2016), which aligns the distributions or the representations across domains by aligning their second-order statistics. On the other hand, Ganin et al., 2016; Shen et al., 2018 or Tzeng et al., 2017 leveraged adversarial training to align the feature representations across domains. More recently, invariant risk minimization (Arjovsky et al., 2019) ensured that the model’s representation is invariant across environments by ensuring the model’s predictions remained the same across domains. This approach, however, required at least two source environments to compute invariant representations and struggled to scale to complex model architectures such as ResNets (Zhou et al., 2022).

Fundamentally, improving the robustness against distribution shifts is a long-tailed problem, meaning that unseen situations eventually arise, and not all situations can be accounted for (Recht et al., 2019; Torralba and Efros, 2011). Therefore, to improve the reliability of deep learning systems and not only their robustness, we need to be able to characterize the representation learned by the model and understand how it is affected by the distribution shifts. To this end, we propose to use explainable artificial intelligence (XAI) methods.

2.3. Explainable artificial intelligence (XAI)

Modern deep learning algorithms are often qualified as black boxes, meaning it is hard to grasp their inner workings fully. This black-box nature limits the applicability of machine learning in safety-critical settings (Achtibat et al., 2022). We can distinguish two main approaches for machine learning explainability: by-design interpretable models and post-hoc explainability (Parekh, 2023). Flora et al. (2022) note that there is no consensus yet in the literature regarding the use of the terms explainability and interpretability. Following Flora et al., 2022, we say that a model is interpretable if it is inherently or by design interpretable, and a model is explainable if we can compute a post-hoc explanation of its decision. By-design interpretability aims at constructing models that are transparent and self-explanatory (Sudjianto and Zhang, 2021), e.g., the decision boundaries of a decision tree. On the other hand, post-hoc explainability seeks to explain a model’s decision by highlighting important features contributing to this decision without explicitly stating how these features affected the model. Methods such as class activation maps (CAMs, Zhang et al., 2017), which plot a heatmap of the important image regions for the classification of this image, fall into this category.

XAI methods for model debugging. One of the main motivations for XAI is to inspect the decision of models to assess whether they rely on relevant factors to make predictions. Several works highlighted biases in the decision process, such as the reliance on spurious features. Lapuschkin et al., 2019 leveraged the GradCAM (Selvaraju et al., 2017) to show how classifiers could rely on watermarks rather than relevant areas of the input image for horses classification, thus highlighting a so-called “Clever Hans” (Pfungst, 1911) effect. CAMs (Zhang et al., 2017) have also been used to understand the behavior of convolutional neural networks (CNNs) in medical imagery classification by Zhang et al., 2021. Another example of the usage of XAI tools to understand and debug a model was proposed by Dardouillet et al., 2023, who leveraged SHapley Additive exPlanations (SHAP, Lundberg and Lee, 2017) to understand a model deployed for oil slick pollution detection on the sea surface. Going one step further, Andeol et al., 2023 recently used conformal predictions to improve the trustworthiness of railway signal detections, a case study where one needs to be sure that the model makes predictions for adequate reasons. In this work, we exploit the complementarities between post-hoc and by-design interpretable XAI methods to provide a thorough understanding of the sensitivity to distribution shifts of CNNs deployed for mapping PV systems from orthoimagery.

3. Data

To analyze the effect of distribution shifts on deep learning models in the context of the remote sensing of rooftop PV systems, we rely on the training dataset Base de données d’apprentissage profond pour les installations photovoltaiques (Database for deep learning applied to PV systems, BDAPPV, Kasmi et al., 2023c). BDAPPV contains nearly 50,000 annotated images of PV systems. A very interesting feature of our case study is that the database contains annotations for 28,000 unique PV systems in France and neighboring countries. The training images were also retrieved from two different sources: satellite images coming from the Google Earth Engine (hereafter referred to as “Google,” Gorelick et al., 2017) and aerial images coming from the IGN (IGN, 2024), the French public operator for geographic information. We have annotations for about 28,000 Google images and 17,000 IGN images. Both providers overlap, meaning we have two annotations for about 7,000 individual PV systems. The dataset is nearly balanced. We refer the reader to Kasmi et al., 2023c for more details regarding the dataset’s characteristics. Figure 1 presents some samples coming from BDAPPV. We refer the reader to Section 4.1 to understand how we used BDAPPV to disentangle the different sources of distribution shifts.

Figure 1. Examples of images of the same PV panels but with different providers and acquisition dates (Up Google, down: IGN).

4. Methods

We aim to explain why convolutional neural networks (CNNs) applied to detect PV panels on orthoimages are sensitive to distribution shifts. We first construct a benchmark to isolate the effect of the three main instances of distribution shifts on orthoimagery highlighted by Tuia et al., 2016 and Murray et al., 2019 using the BDAPPV dataset (see Section 4.1 for more details). These instances include the variability in the geographic location, varying acquisition conditions, and the varying ground sampling distance (GSD).

After quantifying the respective impact on prediction accuracy—measured by the F1 score—we leverage two XAI approaches to understand why these shifts affect the performance. Our working hypothesis is that analyzing the model’s prediction in terms of scales can help us understand why the model is sensitive to distribution shifts. Indeed, scales are located in space and, for each location, correspond to a dyadic partition of the frequency space. Therefore, given each location, we can identify which frequency ranges the model relies on. On the other hand, frequencies are unevenly affected by distribution shifts (for instance, high frequencies are more fragile, Chen et al., 2022). So, scales enable us to assess whether the model focuses on the PV panel to make a prediction and which frequencies it focuses on at this location. Using decomposition in terms of scale is particularly well suited in the case of remote sensing images since the scales, expressed in pixels on images, are indexed in meters and can thus point towards actual elements depicted in the images.

We combine two complementary approaches to explain the model’s decision. Both of these approaches are grounded in the wavelet theory. On the one hand, we leverage a by-design interpretable model, the Scattering transform (Bruna and Mallat, 2013, introduced in Section 4.2.2). We compare the predictions of this model—which are intrinsically interpretable—with those of CNN to see when the predictions are the same and when they differ. On the other hand, we decompose the decision of the model using a post-hoc explainability method, the wavelet scale attribution method (WCAM, Kasmi et al., 2023b), which is a post-hoc explainability method to isolate the important scales in the predictions of our black-box CNN model.

Finally, based on our findings, we propose a data augmentation method to improve the robustness of CNNs, compare our approach with popular domain adaptation methods, and draw some conclusions regarding the choice of image data.

4.1. Disentangling the sources of distribution shifts on orthoimagery

BDAPPV features images of the same installations from two providers and records the approximate location of the PV installations. Using this information, we can define three test cases to disentangle the distribution shifts that occur with remote sensing data: the GSD, the acquisition conditions, and the geographical variability. Natively, our dataset disentangles the effect of the spatial shift, thanks to Google images being roughly geolocalized. Both the resolution and acquisition conditions vary when shifting from Google to IGN images. To disentangle the two sources of shifts, we downsampled the Google images to a GSD of 0.2 m/pixel to match the GSD of the IGN images. We chose not to upsample the IGN images to a GSD of 0.1 m/pixel as it would require adding information to the images and making additional assumptions regarding the method used to carry out the super-resolution task.

We train a ResNet-50 model (He et al., 2016) on Google images downsampled at 0.2 m/pixel of resolution and evaluate it on three datasets: a dataset with Google images at their native 0.1 m/pixel GSD (“Google 0.1 m/pixel”), the IGN images with a native 0.2 m/pixel GSD (“IGN”) and Google images downsampled at 0.2 m/pixel located outside of France (“Google Spatial Shift”). We add the test set to record the test accuracy without distribution shift (“Google baseline”). We only do random crops, rotations, and ImageNet normalizations (i.e., with a mean of [0.485, 0.456, and 0.406] and a standard deviation of [0.229, 0.224, and 0.225]). Figure 2 plots examples of the different test images to disentangle the effects of distribution shifts. The baseline and IGN images represent the same panel at the same spatial resolution. The Google 0.1 m/pixel depicts the same scene but with the native resolution of Google images. Finally, the Spatial shift test set contains images from outside of France.

Figure 2. Test images on which a model trained on Google images (downsampled to 0.2 m/pixel of GSD, “Google baseline”) is evaluated. “Google 0.1 m/pixel” corresponds to the source Google image before downsampling and evaluates the effect of the varying image resolutions. “Google Spatial Shift” corresponds to Google images taken outside of France. “IGN” corresponds to images depicting the same installations as Google baseline but with a different provider.

4.2. Space-scale decomposition of a model’s decision process

4.2.1. Background: the wavelet transform of an image

Motivation and definition. We propose to analyze the decision process of an off-the-shelf CNN model through the lenses of space-scale or wavelet decomposition. Wavelets are a natural tool to decompose an image into scales while maintaining local analysis in space: they provide a single space-scale decomposition. As scales are indexed in terms of actual distances on the ground, we can directly identify the important objects contributing to a model’s decision by studying the important scales. In appendix A, we provide further evidence of the limitation of “traditional” feature attribution methods for explaining the false detection of deep learning models in our use case. Figure 3 illustrates the objects that can be found at different scales of an orthoimage.

Figure 3. Decomposition of a PV panel into scales.

A wavelet is an integrable function $ \psi \in {L}^2\left(\mathrm{\mathbb{R}}\right) $ with zero average, normalized, and centered around 0. Unlike a sinewave, a wavelet is localized in space and in the Fourier domain. This implies that dilatations of this wavelet enable to scrutinize different scales, while translations enable to scrutinize spatial location. In other words, scales correspond to different spatial frequency ranges or spectral domains.

To compute an image’s (continuous) wavelet transform (CWT), one first defines a filter bank $ \mathcal{D} $ from the original wavelet $ \psi $ with the scale factor $ s $ and the 2D translation in space u. We have

(1) $$ \mathcal{D}={\left\{{\psi}_{s,u}(x)=\frac{1}{\sqrt{s}}\psi \left(\frac{x-u}{s}\right)\right\}}_{u\in {\mathrm{\mathbb{R}}}^2,\hskip0.4em s\ge 0}, $$

The computation of the wavelet transform of a function $ f\in {L}^2\left(\mathrm{\mathbb{R}}\right) $ at location $ x $ and scale $ s $ is given by

(2) $$ \mathcal{W}(f)\left(x,s\right)={\int}_{-\infty}^{+\infty }f(u)\frac{1}{\sqrt{s}}{\psi}^{\ast}\left(\frac{x-u}{s}\right)\mathrm{d}u, $$

which can be rewritten as a convolution (Mallat, 1999). Computing the multi-level decomposition of $ f $ requires applying Equation 2 $ J $ times, with $ 1\le s\le J $ . $ J $ denotes the number of levels of decomposition. For each scale, the translation in space $ u $ corresponds to the orientations at a given level.

Mallat (1989) showed that one could implement the multi-level dyadic decomposition of the discrete wavelet transform (DWT) by applying a high-pass filter $ H $ to the original signal $ f $ and subsampling by a factor of two to obtain the detail coefficients and applying a low-pass filter $ G $ and subsampling by a factor of two to obtain the approximation coefficients. Iterating on the approximation coefficients yields a multi-level transform where the $ {j}^{th} $ level extracts information at resolutions between $ {2}^j $ and $ {2}^{j-1} $ pixels. The detail coefficients can be decomposed into various rotations (usually horizontal, vertical, and diagonal) when dealing with 2D signals (e.g., images).

Interpreting the wavelet transform of an image. Figure 4 illustrates how to interpret the (two-level) wavelet transform of an image. Reading is the same for any multi-level decomposition. The right image plots the two-level dyadic decomposition of the original image on the left. Following this transform, the localization on the image highlighted by the red polygon can be decomposed into six detail components (marked yellow and blue) and one approximation component (marked pink). Each detail component has three directions: horizontal, vertical, and diagonal. The yellow components correspond to details at the 1–2 pixel scale, and the blue components to the details at the 2–4 pixel scale. For each location, the wavelet transform summarizes the information in the image at this scale and location.

Figure 4. Image and associated two-level dyadic wavelet transform with indications to interpret the wavelet transform of the image. “Horizontal,” “diagonal,” and “vertical” indicate the direction of the detail coefficients. The direction is the same at all levels.

4.2.2. By design interpretable XAI method: the Scattering transform

The Scattering transform (Bruna and Mallat, 2013) is a deterministic feature extractor. CNNs and the Scattering transform share the same multi-level architecture, where the previous layer’s output is passed onto the next after a nonlinearity is applied. The nonlinearities in a CNN are generally rectified linear units (ReLU), whereas in the Scattering transform, it is a modulus operation. Unlike CNNs, whose kernel coefficients are learned during training, the coefficients of the Scattering transform are fixed. Bruna and Mallat (2013) showed that the Scattering transform computes representations from an input image that share the same properties of translational invariance as the representations computed with a CNN. The advantage of the Scattering transform is that as filters are fixed, we can know precisely what information they extract from the input image. Figure 5 summarizes the feature extraction process of the Scattering transform.

Figure 5. A scattering propagator $ {U}_J $ applied to $ x $ computes each $ U\left[{\lambda}_1\right]x=\mid x\star {\psi}_{\lambda_1}\mid $ and outputs $ {S}_J\left[0/\right]x=x\star {\phi}_{2^J} $ (black arrow). Applying $ {U}_J $ to each $ U\left[{\lambda}_1\right]x $ computes all $ U\left[{\lambda}_1,{\lambda}_2\right]x $ and outputs $ {S}_J\left[{\lambda}_1\right]=U\left[{\lambda}_1\right]\star {\phi}_{2^J} $ (black arrows). Applying $ {U}_J $ iteratively to each $ U\left[p\right]x $ outputs $ {S}_J\left[p\right]x=U\left[p\right]x\star {\phi}_{2^J} $ (black arrows) and computes the next path layer. Figure borrowed from Bruna and Mallat, 2013. Note: In the image, the input $ x $ corresponds to $ f $ and $ \lambda ={2}^jr $ is a frequency variable corresponding to the $ {j}^{th} $ scale with $ r $ rotations.

The input image $ x $ is downsampled, and a wavelet filter $ \phi $ is applied in $ J $ directions. The wavelet coefficients at that scale are retrieved (black arrows), and the image is passed onto the next layer (blue arrows). As the depth increases, the spatial extent covered by the filters decreases. At each spatial location, one takes the modulus of the wavelet transform to compute a scale-invariant representation that indicates the amount of “energy” in the image at this scale and localization.

The Scattering transform is parameterized by the number $ m $ of layers and the number $ J $ of orientations. We have a total of $ mJ+{m}^2J\left(J-1\right)/2 $ coefficients. At the end of the decomposition, the features, i.e., the scattering coefficients, are flattened into a single vector of size $ mJ+{m}^2J\left(J-1\right)/2 $ . We can identify to which scale, location, and orientation on the input image this feature corresponds.

We implement three variants of the Scattering transform with depths $ m $ varying from one to three levels. Bruna and Mallat (2013) stated that first-order coefficients were insufficient to discriminate between two very different images but that coefficients of order $ m=2 $ could. We consider $ J=8 $ orientations. We stack the scattering coefficients into a vector of dimension $ mJ+{m}^2J\left(J-1\right)/2 $ , akin to the penultimate layer of a CNN. We train a linear classifier on this feature vector. Our implementation of the Scattering transform is based on the Python library Kymatio (Andreux et al., 2020).

4.2.3. Post-hoc XAI method: the wavelet scale attribution method (WCAM)

Traditional feature attribution methods (Petsiuk et al., 2018; Selvaraju et al., 2020; Simonyan and Zisserman, 2015) highlight the important areas for the prediction of a classifier in the pixel (spatial) domain. The WCAM (Kasmi et al., 2023b) generalizes attribution to the wavelet (space-scale domain). The WCAM provides us with two pieces of information: where the model sees and what scale it sees at this location. The decomposition of the prediction in terms of scales points towards actual elements on the input image since orthoimagery scales are indexed in meters. For example, on Google images, details at the 1–2 pixel scale correspond to physical objects with a size between 0.1 and 0.2 m on the ground. Thus, we know what the model sees as a panel; we can interpret it and assess whether it is sensitive to varying acquisition conditions. We refer the reader to appendix B or to Kasmi et al., 2023b for more details on the computation of the WCAM.

Reading a WCAM. Figure 6 presents an example of an explanation computed using the WCAM. On the right panel, we can see the important areas in the model prediction highlighted in the wavelet domain. On the left panel, we can see the spatial localization of the important components. We can see two main spatial locations: the center of the image, which depicts the PV panel, and the bottom left, which depicts a pool. Disentangling the scales, we can see that the PV panel’s importance spreads across three scales (orange arrows), while the pool is only important at the 4–8 pixel scale. This underlines that the model focuses on the PV panel because it sees details ranging from small details in the PV modules to the cluster of modules.

Figure 6. Decomposition in the wavelet domain of the important regions for a model’s prediction with the WCAM.

4.3. Improving the robustness through implicit regularization

Improving the robustness to noise and scale perturbations. Since we know that varying acquisition conditions induce perturbations which primarily affect high-frequency components (i.e., the finest scales, Lone and Siddiqui, 2018, we primarily focus on implicit regularization and more precisely data augmentations. Indeed, data augmentation is sufficient to enforce invariance to alterations in the frequency domain. Besides, they are easier to implement for deep learning practitioners and do not require having access to samples from the target domain. For the sake of completeness, we compare our results with explicit regularization techniques. We evaluate popular data augmentation methods to improve the robustness of classification models to image corruptions (Cubuk et al., 2019; Cubuk et al., 2020; Geirhos et al., 2019; Hendrycks et al., 2020, 2022). We consider the AugMix method (Hendrycks et al., 2020) and the recently-proposed RandAugment (Cubuk et al., 2020) and AutoAugment (Cubuk et al., 2019) methods. We refer the reader to appendix D.1 for a detailed presentation of these methods.

Proposed data augmentation methods. As a baseline, we propose blurring the input image and refer to this method as Blurring. We apply a nonrandom Gaussian blur to the image. The value is set by comparing visually Google and IGN images and trying to remove details from Google images that are not visible on IGN images. After a manual inspection, we set the blur level to discard the details at 0.1–0.2 m scale from the image. It corresponds to a blurring value $ \sigma =2. $ in the ImageFilter.GaussianBlur method of the Python Imaging Library (PIL). Our proposed method consists of combining blurring, which removes the small-scale details of the image with a random perturbation of the wavelet transform of the image. We randomly set some wavelet coefficients to 0 sand reconstructed the image from its perturbed coefficients. The perturbation is done across all scales, and the set of coefficients set to 0 is determined using uniform sampling. This results in a random perturbation that removes information for some precise scales and locations. We then reconstruct the image from its perturbed wavelet coefficients. For each call, 20% of the coefficients are canceled. This value balances between the loss of information and the input perturbation. We perturb each color channel independently. The wavelet perturbation aims to disrupt information at specific scales, as it can happen with varying acquisition conditions. The resulting data augmentation method is referred to as Blurring + Wavelet perturbation (WP). Figure 7 presents examples of perturbed images using our method.

Figure 7. Illustration of the effect of our data augmentation method on a sample of images.

Domain adaptation. We complement our analyses by comparing our approach with popular domain adaptation techniques. These techniques are more demanding as unlabeled data from the target domain is required. We refer the reader to appendix E for a discussion of the results obtained with the domain adaptation techniques.

5. Results

5.1. Deep models are mostly sensitive to varying acquisition conditions, leading to an increase in the number of false negatives

Table 1 shows the results of the decomposition of the effect of distribution shifts into three components: resolution, acquisition conditions, and spatial shift. We can see that the F1 score drops the most when the model faces new acquisition conditions. The second most significant impact comes from the change in the resolution. However, the performance drop remains relatively small compared to the effect of the acquisition conditions (which can also be assimilated to variations in the image quality). In our framework, there is no evidence of an effect of geographical variability once we isolate the effects of acquisition conditions and resolution. This effect is probably underestimated, as images of our dataset that are not in France are near France. However, the effect of the acquisition conditions is sizeable enough to seek methods for addressing it.

Table 1. F1 Score and decomposition in true positives, true negatives, false positives, and false negatives rates of the classification accuracy of a CNN model trained on Google images (Google baseline) and tested on the three instances of distributions shifts: GSD (Google 0.1 m/pixel), the geographical variability (Google Spatial Shift) and the acquisition conditions (IGN).

5.2. The Scattering transform shows that clean, fine-scale features are transferable but poorly discriminative

Discriminative and transferable features. In the following, we distinguish between two kinds of features: the discriminative and the transferable features. Discriminative features enable the model to discriminate well between PV and non-PV images. Relying on discriminative features ensures a low number of false positives. On the other hand, transferable features correspond to features that generalize well across domains. If a model relies on transferable features, its performance should remain even across domains. Ideally, we would like a model to rely on discriminative and transferable features.

Accuracy of the Scattering transform. Table 2 presents the accuracy results of the Scattering transform and compares it with a random classifier and the ERM (which is the same model as the one evaluated in Table 1). We can see that the performance on the source domain lags behind the performance of the CNN, but the Scattering transform generalizes better to IGN than the CNN. However, this comes at the cost of a high false positive rate. Table F2 in appendix F.2 presents similar accuracy results for variants of the Scattering transform model in the depth and number of features.

Table 2. F1 Score and decomposition in true positives, true negatives, false positives, and false negative rate of the classification accuracy of the Scattering Transform model trained on Google images and deployed on IGN images. The best results are bolded.

Implications for the CNN. We know which features the Scattering transform relies on. It leverages information at the two-pixel scale after downsampling the input image. In other words, the Scattering transform makes predictions based on clean features at the two-pixel scale. Therefore, we can deduce that these features are transferable, as the performance remains even across datasets, but not very discriminative as the false positives rate is high (across both datasets). Therefore, the analysis of the errors of the Scattering transform and the CNN highlights a potential trade-off between transferable and discriminative features

On the other hand, the CNN should rely on discriminative features located at coarser scales than 8 pixels, and on noisy features. In Section 5.3, we investigate how the distortion of the input image’s coarse scales impacts the CNN’s decision process and the shift in its predicted probability. In Section 5.4.1, we discuss how noise in input images affects the generalization ability of the CNN.

5.3. CNNs are sensitive to the distortion of coarse-scale discriminative features

Predicted probability shifts. The CNN outputs a predicted probability of a PV panel on the input image. When evaluating the CNN on the same scene from two providers, we compute predicted probability shift $ \Delta p=\mid {p}_{ign}-{p}_{google}\mid $ when the model trained on Google is evaluated on IGN images. $ {p}_{google} $ denotes the predicted probability on Google images and $ {p}_{ign} $ on IGN images. By construction, $ \Delta p\in \left[0,1\right] $ . If $ \Delta p=0 $ , the predicted probability did not change when changing the provider. On the other hand, if $ \Delta p\to 1 $ , then it means that the model made a different prediction solely because of the new acquisition condition.

Correlations between the probability shift and low-scale similarity of the images. For all images in our test set ( $ n=4321 $ ), we compute the similarity between the low-scale components of the input image across the two domains. This enables us to assess how similar images depicting the same scene on Google and IGN are, with respect to their low-scale components, i.e., components larger than 8 pixels, which correspond to the approximation coefficients of a 3-level dyadic decomposition of the image.

On the other hand, we compute the predicted probability shift for each image across two domains. The predicted probability shift indicates how much the model’s prediction changed when facing the IGN image.

Suppose the CNN is indeed sensitive to low-scale perturbations of the input image. In that case, we expect a correlation between the dissimilarity between the approximation coefficients (which only contain the low-scale components of the image) and the predicted probability shift (which indicates whether the model changed its prediction once faced with a new image).

We evaluate the similarity between the approximation coefficients using two metrics: the Structural similarity index measure (SSIM, Wang et al., 2004) and the Euclidean distance between the approximation coefficients. The SSIM takes values between −1 and 1, where 1 indicates perfect similarity, 0 indicates no similarity, and −1 indicates perfect anti-correlation. On the other hand, the Euclidean distance takes positive values; the greater the distance, the greater the dissimilarity between the images.

We evaluate the correlation between the similarity of the approximation coefficients and the magnitude of the probability shift using the Pearson correlation coefficient (PCC, Pearson and Galton, 1895). The PCC is a statistical measure that quantifies the strength and direction of a linear relationship between two continuous variables. The PCC value ranges from $ -1 $ to $ 1 $ , where −1 indicates a perfect negative linear relationship, 1 is a perfect positive linear relationship, and 0 is no linear relationship (the variables are uncorrelated). Given two random variables X and Y, the PCC is given by

$$ r=\frac{\mathrm{Cov}\left(X,Y\right)}{\sigma_X{\sigma}_Y}, $$

where $ \mathrm{Cov}\left(X,Y\right) $ denotes the covariance between $ X $ and $ Y $ and $ \sigma $ is the standard deviation. In addition to computing the PCC, we report its p-value to assess whether the reported value significantly differs from 0, thus rejecting the hypothesis that the variables are uncorrelated.

As expected, we obtain a negative Pearson correlation coefficient equal to −0.41 (with a p-value $ <{10}^{-5} $ ) between the input images’ SSIMs and the predicted probability shift. Using the Euclidean distance, we obtain a correlation coefficient of 0.250 ( $ p<{10}^{-5} $ ). These results back the idea that the CNN is sensitive to low-scale perturbations of the input image, which results in a shift in the predicted probability.

Visualization of the model’s response with the WCAM. The WCAM disentangles the important scales in a model’s prediction. It enables us to see which scales were disrupted. On Figure 8, we present an example of an image that was initially identified as a PV panel but turned out to be no longer recognized on IGN images.

Figure 8. Analysis with the WCAM of the CNNs prediction on an image no longer recognized as a PV panel.

We can see that in both cases, the approximation details are important in the model’s prediction. The model responds to distortions at this scale by no longer focusing on a single area. Indeed, the model weights more components located at the 2–4 and 4–8 pixel scale (orange circles), which were not as important initially. At the level of the perturbed scales, we can also witness that the model is disrupted by factors lying next to the PV panel (green circle). We supply more examples of such cases in appendix G and discuss the quantitative analysis of this result in appendix C.

5.4. Pathways towards improving the robustness to acquisition conditions

5.4.1. Blurring and wavelet perturbation improve accuracy

Table 3 reports the results of our data augmentation techniques and compares them with existing methods. We can see that augmentations that explicitly discard small scales (high frequencies, i.e. Blurring and Blurring + WP) information perform the best. However, the blurring method sacrifices the recall (which drops to 0.6) to improve the F1 score. In Table 3, this can be seen by the increase in false positives rate. Therefore, this method is unreliable for improving the robustness to acquisition conditions. We recall that the true positive rate and the false negative rate divide the number of true positives (resp. false negatives) by the number of positive samples. Similarly, the true negative and false positive rates divide the number of true negatives (resp. false positive) by the number of negative samples in the dataset. The true positive rate corresponds to the recall, and the true negative rate to the specificity.

Table 3. F1 Score and decomposition in true positives, true negatives, false positives, and false negatives rate for models trained on Google with different mitigation strategies. Evaluation on IGN images. The Oracle corresponds to a model trained on IGN images with standard augmentations. Best results are bolded, second-best results are underlined, values highlighted in red indicate the worst performance, and values in orange indicate the second-to-last worst performance

On the other hand, adding wavelet perturbation (WP) contributes to restoring the accuracy of the classification model without sacrificing the precision or the recall. While the drop in accuracy is still sizeable compared to the Oracle, the gain is consistent compared to other data augmentation techniques. Compared to RandAugment, the best-benchmarked method, our Blurring + WP is closer to the targets regarding true positives and true negatives and makes lower false negatives. This experiment shows that it is possible to consistently and reliably improve the robustness of acquisition conditions using a data augmentation technique, which does not leverage any information on the IGN dataset.

5.4.2. On the role of the input data: towards some practical recommendations regarding the training data.

Generalizability of the feature representation. Our results show that lowering the reliance on high-frequency content in the image improves generalization. This content is located on the 0.1-0.2 m scale and only appears on Google images. In Table 4, we flip our experiment to study how a model trained on IGN images generalizes to Google images. Results show that the model trained on IGN generalizes better to the downscaled Google images than the opposite. This result further supports the idea that higher resolution is not necessarily better for good robustness to acquisition conditions.

Table 4. F1 Score and true positives, true negatives, false positives, and false negatives rates. Evaluation computed on the Google dataset. ERM was trained on Google and Oracle on IGN images

Reliability trade-offs. The training data is often considered as given in many practical settings, especially given the high cost of annotating samples. This motivates the use of domain adaptation techniques, such as those described in this work. From Tables 3 and E1 in appendix E, it should be noted that the F1 score can be misleading regarding how the different methods attenuate the effects of the distribution shifts. In particular, it should be noted that Blurring achieves a very high F1 score at the expense of the number of false positives. On the other hand, the domain adaptation method Wasserstein Distance Guided Representation Learning (WDGRL) exhibits a relatively false negative rate. Depending on the task at hand, different methods can be preferred. In the case of the remote sensing of rooftop PV systems, false negatives can be an issue, thus leading to favor solutions such as Blurring or Adversarial Discriminative Domain Adaptation (ADDA), even though we know that these methods generate a lot of false detections.

6. Discussion

6.1. Conclusion

This work aims to explain why convolutional neural networks (CNNs) applied to detect PV panels on orthoimages are sensitive to distribution shifts. We first set up an experiment to disentangle the effects of the three main distribution shifts occurring in remote sensing (Murray et al., 2019; Tuia et al., 2016), namely geographical variability, varying acquisition conditions, and varying resolution. We showed that the varying acquisition conditions contribute significantly to the observed performance drop. To explain why this drop occurs, we leverage space-scale analysis to disentangle the different scales from the input images. We combine two types of explainable AI methods grounded in the wavelet decomposition of the input images to show that the CNN relies on noisy features (at the finest scales) and features that are not very well transferable across domains (at the coarsest scales).

We then introduced a data augmentation technique to improve the model’s robustness to distortions of the coarse-scale features and remove noise from the fine-scale features. We compare this method against popular data augmentation techniques and show that our approach outperforms these baselines. We also compared our approach with more demanding domain adaptation techniques and showed that our approach remains competitive. We then discussed several practical takeaways of this study for the training or the choice of the training data for the initial training of the deep learning model.

Broader impact. Mapping rooftop PV systems is a recurring issue in many countries, and a lot of actors interested in such rooftop PV registries require reliable data (De Jong et al., 2020; Kasmi, 2024). While offering the possibility to quickly and cheaply map PV systems over vast areas, current methods for mapping rooftop PV installations lack reliability owing to their poor generalization abilities beyond their training dataset (De Jong et al., 2020). This work addresses this gap and thus demonstrates that remote sensing of PV installations is a reliable way to construct registries of rooftop PV systems.

The methodology introduced in this work, which consists of first isolating the main source of performance drop among possible types of distribution shifts, then leveraging XAI methods to grasp better the impact of these shifts on the model’s predictions to finally highlight how invariance to these shifts can be mitigated can be replicated to other case studies.

6.2. Limitations and future works

Further discussion of the geographical variability. Our training data was limited to a narrow area around France. Therefore, we suspect the effect of the geographical variability to be underestimated. For instance, Freitas et al., 2023 showed that fine-tuning a model with data that is not far from the target area (e.g., France when the goal is to map PV systems in Portugal) enables accuracy gains compared to directly transferring a model trained over the United States. It could be interesting to study how the performance varies with the distance between the training data and the target mapping area once all other factors (acquisition conditions, resolution) are accounted for.

Extensions to other models. Over the last couple of years, foundation models (Bommasani et al., 2022) have been redefining the standards in deep learning. These very large models, trained on large data corpora, have shown remarkable performance for many challenging tasks, especially for text (Brown et al., 2020) and image (Rombach et al., 2022) generation. These models are used for more conventional and specialized tasks such as image segmentation (Kirillov et al., 2023) and achieve superior performance to conventional approaches while only requiring a few samples to learn their new task. Extending this benchmark and evaluating the performance of foundation models fine-tuned for segmenting PV panels, such as Yang et al., 2024, under distribution shifts could be interesting.

Application to other case studies. The key ingredients to replicate our methodology are the disentanglement of the effect of the different types of shifts occurring in the case study at hand and the combination of various XAI methods to build a relevant intuition regarding the characterization of the feature representation of the model, in terms of semantic relevance and in terms of robustness to the types of shifts that are recurring in the case study. Introducing new domain adaptation methods should not be prioritized over thorough decomposition of the distribution shifts and analysis of their effects.