Location-aware Adaptive Normalization: A Deep Learning Approach For Wildfire Danger Forecasting

Wildfire forecasting or wildfire-susceptibility mapping from remote sensing and Earth observations data is a very important topic for wildfire management[11]. We briefly review some prior related works in this direction. Iban et al. [19], Pham et al. [20], and Gholami et al. [21] relied on traditional machine learning approaches to generate susceptibility maps. Shang et al. [22] and Mitsopoulos and Mallinis [23] utilized Random Forests classifiers (RF). In their works, they studied the importance of biotic and abiotic predictors for wildfire forecasting. Jiang et al. [24] proposed a deep learning approach based on a Multi-Layer Perceptron (MLP) and included a comparison with traditional machine learning algorithms. In Le et al. [25], a similar MLP-based approach was presented to generate a forest fire danger map. Zhang et al. [26] used a convolutional neural network (CNN) and extended their work later to predict fire susceptibility at the global level [27]. Other works with CNN were conducted in Bjånes et al. [28] and Bergado et al. [29]. Furthermore, Huot et al. [30] approached the problem as a scene classification task using U-Net models to predict wildfire spreading. Their approach operates directly on the whole scene. A similar approach based on a U-Net++ model for global wildfire forecasting was proposed in Prapas et al. [31]. Yoon and Voulgaris [32] presented an approach that relies on a recurrent network with Gated Recurrent Units (GRU) to model past observations and on a CNN to predict wildfire probability maps for multiple time steps. More recently, Prapas et al. [16] and Kondylatos et al. [13] proposed to use LSTM-based (Long-Short Term Memory) approaches. They exploited both temporal and spatio-temporal context by applying recurrent LSTM and ConvLSTM models. They did not consider the whole scene at once but rather the classification of one pixel at a time (pixel-level).

Unlike these works, we do not treat all observation variables in the same way, but we propose a deep learning model that handles different types of variables in separated 2D and 3D CNN branches. In contrast to a ConvLSTM, which models spatial and temporal relations separately, a 3D CNN models spatio-temporal relations jointly. We assume that the dataset contains static and dynamic variables, which we argue is the case for most datasets. Similar to Prapas et al. [16] and Kondylatos et al. [13], we also formulate the problem as pixel-level classification taking into account the spatio-temporal local context around the target pixel.

When deep learning is applied to potentially multi-source remote sensing-based Earth observation data, multi-branch neural networks are a commonly used framework. This is mainly because such networks enrich representation learning and provide discriminative learning perspectives of the input variables [33]. In addition, an important aspect of the multi-branch design is the capability to adapt some parts of the model to a specific type of input. The general framework generates features from each branch and fuses these features in the network to obtain a unified feature vector. This fused representation is used as input to the subsequent layers. In Gaetano et al. [34], a two-branch 2D CNN network was proposed to handle panchromatic information along with a multi-spectral one for image classification. Tan et al. [35] reduced the depth of a semantic segmentation classifier by applying consecutive blocks, each containing three CNN branches. A similar objective can be found in Zhao et al. [33], where the network complexity was reduced via weight sharing and self-distillation (SD) embedding. In this way, only the main branch is used during inference, which inherits the knowledge of trained subbranches and has a close performance to an ensemble model. For hyperspectral image classification, Xu et al. [36] introduced a model called Spectral–Spatial Unified Network (SSUN). In their model, spectral features are learned by a grouping-based LSTM, and spatial features are learned by a 2D CNN. Shen et al. [37] used separated spectral and spatial convolutional branches for hyperspectral input ($\text{S}^{2}$CDELM). They based their framework on the extreme learning machine (ELM). Unlike common backpropagation algorithms, they used a single hidden layer feed-forward model. A multi-branch architecture was also explored for image fusion. Liu et al. [38] proposed a two-stream CNN called (StfNet). They investigated the task of spatio-temporal image fusion. Their network takes a coarse image input along with its neighboring images to predict the reconstructed fine image. Some works adapted a multi-branch architecture to construct a multi-scale feature vector. In Gan et al. [39], a dual-branch CNN with different filter kernel sizes was used as an autoencoder. Thus, the input image could be processed at different scales. Tang et al. [40] proposed a multi-scale Gaussian pyramid to handle hyperspectral input. They used the Gaussian pyramid to obtain multi-scale images which are then processed by ResNet modules [41]. In this way, spatial features can be learned at different scales. They further used a second branch, which performs a discrete wavelet transform on the spectral input followed by an LSTM module [42]. The spatial and spectral features are then fused and processed by an MLP to obtain the final classification result. For short-term multi-temporal image classification, Zheng et al. [43] addressed the task through a Multi-temporal Deep Fusion Network (MDFN). In their framework, the LSTM-based branch is used to learn temporal-spectral features. At the same time, temporal-spatial and spectral-spatial information is learned via a joint 3D-2D CNN with two branches. Furthermore, some works employed attention mechanisms with multi-branch architectures [44, 45, 46]. When attention is applied, it drives the model to focus more on regions of interest. The former described methods [34, 35, 36, 37, 38, 39, 40, 43, 44] except of Zhao et al. [33], Zhu et al. [45], and Deng et al. [46] did not consider linking information between branches during the feature learning stage. This limits the gradient flow and disentangles the correlations between the learned features.

This paper proposes an architecture composed of two CNN branches for a forecasting task. A 2D branch is used to learn spatial features from static variables. At the same time, a 3D branch is used to learn spatio-temporal features from dynamic variables, which vary along the input temporal dimension. The branches are further linked via adaptive modulation layers to model the causal effects of static variables on dynamic variables.

Since the introduction of normalization techniques in deep learning, they became a basic building block in many state-of-the-art models. Common normalization methods include batch normalization [47], group normalization [48], instance normalization [49], and layer normalization [50]. It has been shown empirically that normalization layers help with model optimization and regularization. Through normalization layers, the activation maps inside the model are normalized to follow a normal distribution with zero mean. After that, the normalized activation maps are modulated or denormalized by learnable affine transformation parameters. These parameters vary across channels and are learned based on the running training statistics together with the model parameters. Therefore, such a normalization method is called unconditional. Compared to popular unconditional normalizations, there exist conditional normalization techniques which aim to learn affine parameters conditionally on external input. In the field of computer vision, conditional normalization is often used for image synthesis and style transformation [51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64]. More recently, Marín and Escalera [65] adapted the conditional normalization from Wang et al. [52] and Park et al. [51] to generate high-resolution satellite images. We use the conditional normalization in a very different way than [65, 52, 51]. While these works focus on synthesizing an image using the segmentation map as conditional input, we aim to modulate dynamic features conditioned on static features.

A plethora of studies exist about temporal modeling in remote sensing [66, 67, 68, 69]. Recently, the self-attention model, also known as Transformer and first presented by Vaswani et al. [17] for natural language processing, has become a natural choice to handle sequential data, which includes a positional encoding. In the field of remote sensing, many works showed the benefits of adapting positional encoding for time-dependent image classification (Garnot et al. [70], Garnot and Landrieu [71]), panoptic segmentation (Garnot and Landrieu [72]), and image generation (Dress et al. [73]). In their work, each image was given an encoded time vector according to its position with respect to a reference point, i.e., the first acquisition time step. Nyborg et al. [74] used the calendar time (day of the year) to provide positional information within the sequence. They also proposed to learn or estimate time shifts between geographically distant regions to enhance the generalization further. In another work, Nyborg et al. [75] used the thermal time, which is obtained by accumulating daily average temperatures over the growing season, for crop classification. In general, transformer-based approaches require the positional encoding since the temporal information is otherwise lost within a transformer model. While 3D CNNs consider the temporal order of the input such that a positional encoding is not necessary, we show in this work that adding an absolute temporal encoding is also useful for time-dependent forecasting.

As shown in Fig. 1, our network consists of two branches that process dynamic and static variables, respectively. We denote the data cube with dynamic variables by $\mathcal{C}(\{V_{d},T,W,H\})$ and the data cube with static variables by $\mathcal{C}(\{V_{s},W,H\})$. As in previous works, we normalize the input channel-wise to the range $[0,1]$. Since the static variables do not have a time component, we use 2D convolutions for the static branch and 3D convolutions for the dynamic branch. More in detail:

Dynamic branch. The dynamic branch takes the variables $V_{d}$ which vary over time as input. It consists of $3$ blocks; each block has a 3D convolution with a $3\times 3$ kernel size followed by a ReLU activation function and a 3D max pooling layer. To reduce overfitting, we use global average pooling (GAP) [88] at the end of the last block. We denote the feature vector learned from this branch as $\textbf{X}_{d}\in\mathbb{R}^{256}$.

Static branch. In parallel to the dynamic branch, the static branch has a similar architecture. However, 2D convolutions are used instead of 3D ones. We denote the feature vector learned from this branch as $\textbf{X}_{s}\in\mathbb{R}^{128}$. Note that the dimensionality of the static feature vector is lower than the dimensionality of the dynamic feature vector since the input data cube is smaller.

In a nutshell, the dynamic- and static branch are two functions $f_{d}$ and $f_{s}$, respectively:

For the dynamic feature vector we add an absolute temporal encoding $\textbf{X}_{\tau}$, which will be described in Section IV-C. The two feature vectors are then concatenated and fed into 4 classification layers with 1D convolutions of kernel size 1. The layers reduce the dimensionality from 384 to 256, 128, 32, and 2. To reduce overfitting, we use dropout with a dropout probability $p=0.5$ for the 1D convolutional layers except the last two layers. Finally, a softmax activation is used after the last classification layer to predict the probability of a wildfire. In addition, we use a batch normalization layer [47] for the $1^{st}$ block of each branch. More implementation details are given in Section V.

For training, we use the binary cross entropy as loss function:

where $N$ denotes the batch size and $\hat{Y}_{T{+}1}^{(n)}\in\{0,1\}$ is the true label for sample $n$.

In the following, we discuss the Location-Aware Adaptive Normalization (LOAN) that modulates the dynamic features based on the static features and the already mentioned absolute temporal encoding.

In general, dynamic variables are correlated with the geographic location, i.e., temperature and pressure change with elevation, soil moisture and NDVI vary with land cover, and humidity is correlated with some static variables like the waterway distance. Since the dynamic variables depend on the static variables and not vice versa, we aim to exploit this knowledge in our approach. This is done by learning to normalize the dynamic features based on the location-dependent static features. To this end, we propose a conditional normalization technique for remote sensing data called Location-aware Adaptive Normalization (LOAN).

We first describe a batch-normalization [47] where the activation map is normalized before it is modulated by scale $\gamma$ and bias $\beta$. Let ${z_{d}}^{i}\in\ \mathbb{R}^{N\times K^{i}\times D^{i}\times W^{i}\times H^{i}}$ be an activation map in the $i\text{-th}$ block of the dynamic branch and ${z_{s}}^{i}\in\mathbb{R}^{N\times K^{i}\times W^{i}\times H^{i}}$ be an activation map of the corresponding $i\text{-th}$ block in the static branch, where $K^{i}$ denotes the number of channels and $D^{i}$, $W^{i}$ and $H^{i}$ denote the depth, width, and height of the activation map $z^{i}$, respectively. Using the indices $n\in\{1,\ldots,N\}$, $k\in\{1,\ldots,K^{i}\}$, $t\in\{1,\ldots,D^{i}\}$, $w\in\{1,\ldots,W^{i}\}$, and $h\in\{1,\ldots,H^{i}\}$, the normalization of the dynamic branch is performed by the following equation:

where ${z_{d}}_{(n,k,t,w,h)}^{i}$ is the activation before the normalization, $\mu_{k}$ and $\sigma_{k}$ are the computed mean and standard deviation of channel $k$, i.e., computed over the tensor $D^{i}\times W^{i}\times H^{i}$ and all samples $n$ in the batch, and $\gamma_{k}$ and $\sigma_{k}$ are the learnable modulation parameters.

In our case, we aim to learn a modulation of the dynamic features ${z_{d}}_{(n,k,t,w,h)}^{i}$ at the $i\text{-th}$ block where the modulation parameters ${\gamma_{s}}_{(n,k,w,h)}^{i}$ and ${\beta_{s}}_{(n,k,w,h)}^{i}$ depend on the corresponding static features ${z_{s}}_{(n,k,w,h)}^{i}$:

In contrast to (5), the modulation parameters ${\gamma_{s}}_{(n,k,w,h)}^{i}$ and ${\beta_{s}}_{(n,k,w,h)}^{i}$ vary with respect to sample $n$ in the batch, the location $(w,h)$, and across channel $k$, but they are constant over time $t$. Furthermore, they are conditioned on the static features ${z_{s}}^{i}$. We thus call ${z_{s}}^{i}$ the conditional map for the modulation.

The way how ${\gamma_{s}}_{(n,k,w,h)}^{i}$ and ${\beta_{s}}_{(n,k,w,h)}^{i}$ are computed is illustrated in Fig. 2. First, the conditional map ${z_{s}}^{i}$ is normalized channel-wise as following:

where

Afterward, ${\hat{z_{s}}}_{(n,k,w,h)}^{i}$ is projected by two convolutional layers with $K^{i}$ filters to compute ${\gamma_{s}}_{(n,k,w,h)}^{i}$ and ${\beta_{s}}_{(n,k,w,h)}^{i}$. In our implementation, these modulation parameters are then duplicated along the temporal dimension to match the depth $D^{i}$ of ${z_{d}}^{i}$ such that (6) can be computed.

We add the Location-aware Adaptive Normalization layer (LOAN) in the first two blocks as shown in Fig. 1. The activation maps of the dynamic branch are normalized only in the $1^{\text{st}}$ block and modulated in both the $1^{\text{st}}$ and $2^{\text{nd}}$ blocks. The impact of the blocks where the LOAN layer is added is evaluated in Table IV.

In the experimental section, we also evaluate a variant of LOAN that is not conditioned on the intermediate features of the static branch as shown in Figs. 1 and 2, but on the static variables directly as shown in Fig. 3. In this case, the conditional map has different spatial dimensions and number of channels compared to the features in the dynamic branch, i.e., the conditional maps consist of $C$ variables and have $W^{s}\times H^{s}$ spatial dimensions. In this respect, the conditional map ($W^{s}\times H^{s}$) is first resized to match the spatial dimensions ($W^{i}\times H^{i}$) of the activation map from the dynamic branch. We use the nearest-neighbor method for the down-sampling. The resized conditional map is then fed into a convolutional layer with $3\times 3$ kernel size to double the number of channels, i.e., $2\times C$. Finally, as in the previous version of LOAN, the conditional map is normalized, projected by two convolutional layers, and duplicated along the temporal dimension to compute ${\gamma_{s}}_{(n,k,w,h)}^{i}$ and ${\beta_{s}}_{(n,k,w,h)}^{i}$. The impact of different conditional maps using the variants of LOAN is evaluated in Table II.

Some extreme events in the climate model have a dependent relation on time [89, 14]. This is also the case for the FireCube dataset [18] where wildfire events vary from month to month and occur more frequently in the summer time as shown in Fig. 5. So far, the network does not consider an absolute time like the day of the year. Instead, for any forecast day $T{+}1$, the last 10 days are used as observations but the network does not have the information what day during the year $T$ is.

As shown in Fig. 1, we add this information to the dynamic branch before we concatenate the static and dynamic features. To encode the day of the year, we use the standard fixed sinusoidal-based encoding by Vaswani et al. [17]. We pre-compute for each day of the year $\tau{\in}\left[0,365\right]$¹¹1We consider February $29$ for the encoding., which is extracted from $T$, the absolute temporal encoding vector $\textbf{X}_{\tau}\in\mathbb{R}^{256}$:

where $j$ is the embedding dimension. Each even dimension results from a sine function, while each odd dimension results from a cosine function. This allows $\tau$ to have a smooth and yet unique encoding for every time step, i.e., each day of a year. Note that the vector has the same size as $\textbf{X}_{d}$.

In order to add the absolute time encoding vector $\textbf{X}_{\tau}$ to the dynamic feature vector $\textbf{X}_{d}$, we weight each element of the vector by a learnable weight vector $\textbf{W}\in\mathbb{R}^{256}$:

where $\circ$ denotes the Hadamard product, i.e., element-wise multiplication. Fig. 4 illustrates how the temporal embedding is added to the model.

We compare our approach to the approaches that have been evaluated on the described dataset in Kondylatos et al. [13]. This includes two deep learning models, namely LSTM [42] and ConvLSTM [66], and two classical machine learning classifiers, namely Random Forests (RF) [94] and XGBoost [95]. For further details regarding the architectures and hyper-parameters of the models, we refer to the work of Kondylatos et al. [13]. In order to demonstrate the benefit of treating static and dynamic variables differently, we also compare with a one-branch 3D CNN where we duplicate the static variables along the temporal dimension and concatenate them together with the dynamic variables to form a single data cube. In addition, we compare to the recent transformer models TimeSformer [92] and Video SwinTransformer [93], which use space-time attention. As mentioned in Subsection II-D, vision transformers are based on a self-attention mechanism to model spatio-temporal dependencies. Compared to CNNs, transformers have less inductive bias and need much more computational resources for training. To the best of our knowledge, no prior work has done a systematic study about the performance of video vision transformers for wildfire forecasting.

To ensure a fair comparison, all baselines were re-implemented and trained on the same samples with a fixed random seed. We do not use any augmentation technique. The quantitative results of our experiments are provided in Table LABEL:tab:table1. The results of the proposed 2D/3D CNN are shown with and without absolute temporal encoding (TE).

We can observe that the proposed 2D/3D CNN outperforms the other methods for most metrics, particularly FP, TN, Precision, F1-score, and OA, on the validation and testing sets. SwinTransformer and ConvLSTM achieve a slightly higher AUROC for 2019 and 2020, respectively. LSTM and SwinTransformer achieve a higher recall, but at the cost of a very low precision. In comparison with other deep learning methods, LSTM and SwinTransformer have even the highest number of false positives for all years. The main weakness of LSTM lies in the fact that it does not consider the spatial context, while large models like SwinTransformer are prone to overfitting, which results in a relatively poor performance for 2020. RF and XGBoost have the same disadvantage as LSTM, but even a weaker temporal model and thus perform worse than LSTM. Most interesting is the comparison to 3D CNN since it uses the same 3D CNN structure but only one branch, i.e., it treats static variables like dynamic variables. The results show that the proposed approach with two branches outperforms the single-branch architecture for all metrics and all years. This demonstrates the importance of treating static variables differently than dynamic variables. Adding the absolute temporal encoding (TE) to the model substantially reduces FP at the cost of decreasing TP. This is also reflected in the precision and recall.

In Table LABEL:tab:table2, we present additional experimental results alongside the memory footprint as the number of parameters (# Params), the estimated multiply-accumulate operations (multiply-adds) (MMACs), and the expected inference time, which is estimated as samples per millisecond (# SPmS). The performance metrics are calculated on both testing years $2020$-$2021$ as one set. Since the LOAN layers increase the amount of parameters of the 2D/3D CNN, we report the results of the one-branch 3D CNN using $323$k and $499$k parameters. The smaller 3D CNN has about the same amount of parameters as the 2D/3D CNN without LOAN, whereas the larger 3D CNN has more parameters than the proposed model. Due to the lack of a spatial modeling, LSTM has the fewest parameters and is the fastest, but the precision is very low. ConvLSTM, the small one-branch 3D CNN, and 2D/3D CNN without LOAN and TE, which consider the spatial context, perform similar but 2D/3D CNN is the fastest approach and ConvLSTM is the slowest approach among them. Compared to the other approaches, transformer models have considerably more parameters and require more operations, which makes them computationally expensive. Nevertheless, our approach outperforms both transformer models while requiring by far less parameters and computational operations.

Normalizing the dynamic features conditioned on the static features (LOAN) increases all metrics. It also outperforms the large 3D CNN in all metrics and inference time. Adding TE increases all metrics when LOAN is not used, while increasing the computational cost only very little. When LOAN is used, adding TE decreases the recall but increases all other metrics. Since LOAN and TE change the dynamic features, we observe a different trade-off between recall and precision if both are used. This change is consistent over the years as shown in Table LABEL:tab:table1. Nevertheless, the F1-Score, AUROC, and OA are highest if both are used.

To assess the importance of different static variables, we present the results obtained with different combinations of static variables in Table I. For this experiment, we use 2D/3D CNN with LOAN but without TE. All dynamic variables are used in this experiment and the results are reported for the year $2019$ and for the years $2020$-$2021$ as one set. The static variables are grouped into 3 main categories: topographic variables consisting of digital elevation model (DEM) and slope, anthropogenic-related variables consisting of distance to roads or waterway and population density, and land cover variables.

From the results in Table I, we can conclude that among the 3 categories topographic variables give the best results for the years $2020$-$2021$ when they are used without other variables. While land cover and anthropogenic-related variables provide the best results for the year $2019$. Overall, all static variables are relevant and the best results are obtained when all static variables are used (last row). This is also better than using the static variables of 2 out of the 3 categories, which is reported in rows 4-6 of Table I.

While Table I shows the importance of different static variables as input to the 2D/3D CNN, we also analyze the impact of different ways to modulate the dynamic features in Table II. For the experiments, we use the 2D/3D CNN without TE and all dynamic and static variables as input. While we use all variables, the different settings differ in the input that is fed to the LOAN layer, i.e., the static features that the modulation of the dynamic features is conditioned on.

The results of the proposed conditioning, where we use the features from the corresponding block of the static branch, are shown in the last row. In the first row, we show the results if we do not use LOAN at all, i.e., we do not use any modulation of the dynamic features. For the other rows of Table II, we modify LOAN such that it is not conditioned on the intermediate features of the static branch as shown in Figs. 1 and 2. Instead, we condition LOAN directly on static variables. Note that we need to slightly adapt LOAN as shown in Fig. 3 since the number of static input variables $C$ differs from the number of feature channels $K$ at the block where LOAN is added.

As seen from Table II, the best result is obtained when we use the activation maps, i.e., the intermediate features, from the static branch for conditioning the modulation. However, comparable results are obtained when all static variables are directly used by the variant of LOAN that is shown in Fig. 4. While this variant achieves slightly higher AUROC, the variant shown in Fig. 3 achieves higher F1-score and OA. Using the static variables directly, we can analyze how the three categories of static variables impact the modulation of the dynamic features and thus the results. We can conclude that the modulation of the dynamic features is very sensitive to its condition. For the year $2019$, all three categories (rows $2$-$4$) improve the results compared to the setup without feature modulation (first row). For the years $2020$-$2021$, this is not the case and only the combination of all static variables (row $5$) leads to an improvement. The reason is the mismatch between $C$ and $K^{i}$. Compared to the other categories, land cover has the highest number of variables per pixels ($C{=}10$) and shows the best performance. This indicates that conditioning the modulation of the dynamic features on the intermediate features of the static branch is a more practical approach than conditioning the modulation on the static variables directly, which seems to be sensitive to the number of variables.

Predicted wildfire danger-susceptibility maps are depicted in Figs. 6, 7, and 8. We take input from the $1^{st}$ and $2^{nd}$ days of three consecutive months in summer (June, July, and August) and predict for the $2^{\text{nd}}$ and $3^{rd}$ days of each month. We end up with around $500$k pixels (samples) per day. The output from the deep learning models (LSTM, ConvLSTM, Ours) is a probability $Y\in[0,1]$. In addition, we visualize the predicted maps produced by FWI with the provided spatial resolution $8\text{ km}\times 8\text{ km}$. The output of FWI is clipped to the range $[0,50]$ [13]. The ignition points of large wildfires at those days are represented as black circles on the map. The first observation is that regardless of the coarse resolution of FWI, the predicted maps produced by the deep learning models are more reliable. While FWI relies on meteorological observations and models functional relationships, the results show that the modeled functional relationships are insufficient and do not reflect the complexity of the problem of forecasting wildfire. We also find that our proposed model with TE discards spots that result in a high false positive error rate while it keeps extreme ones (cf. the results for July in $2020$ and $2021$). Another important observation is that the LSTM model, which does not account for the spatial context, tends to produce heterogeneous predictions where neighboring pixels often have very different wildfire danger probabilities. Consequently, it generates many false positives. In contrast, our proposed model and ConvLSTM produce more homogeneous and clustered predictions. The respective performance metrics for Figs. 6 and 7 are provided in Table III.

As shown in Fig. 1, the proposed network has three blocks and we add LOAN to the first and second block. We evaluate in Table IV different configurations where we add LOAN only to the first or to all three blocks. The results are reported without TE. If we add LOAN only to the first block, the performance increases for the year $2019$ but not for the years $2020$-$2021$ compared to our model without LOAN (first row). When adding LOAN to the first two blocks, we observe a consistent improvement for all years. For the year $2019$, Precision, Recall, F1-score, and AUROC are improved by $+3.96\%$, $+12.39\%$, $+8.75\%$, and $+1.14\%$, respectively, and by $+0.36\%$, $+3.62\%$, $+2.02\%$, and $+0.30\%$ for the years $2020$-$2021$, respectively. Adding LOAN to all three blocks performs worse than adding LOAN only to the first two blocks. This is due to the decrease of spatial resolution after each block by the pooling layers. In the third block the spatial resolution is too coarse for a location-specific modulation of the dynamic features.

As we have seen in Tables LABEL:tab:table1 and LABEL:tab:table2, absolute temporal encoding (TE) increases precision at the cost of lower recall. Depending on the use of the wildfire forecasting, recall or precision are more important. The precision also depends on the amount of negative samples. In order to show that TE gives consistently a higher precision, we varied the number of negative samples. In this experiment, we test on all positive samples in the test set (years $2020$-$2021$) and gradually increase the number of negative ones. As shown in Fig. 9, we start with a setup where the number of negative samples is equal to the number of positive samples, i.e., 5635 positive and 5635 negative samples. We then increase the number of negative samples. Since the number of negative samples increases, the precision decreases for all methods. Note that the recall does not change since the number of positive samples remains the same. As already observed in Table LABEL:tab:table2, LSTM has a very low precision. 2D/3D CNN with LOAN has in all settings a higher precision than ConvLSTM and a much higher recall as shown in Table LABEL:tab:table2. Transformer models on the other hand have less precision than ConvLSTM but provide an overall better recall as shown in Table LABEL:tab:table2. While adding TE decreases recall, it improves the precision substantially and the improvement increases when the number of negative samples increases. While other metrics like F1-score or AUROC combine precision and recall in a single measure, depending on the application a higher recall or a higher precision might be more important. If precision is more important, TE is very useful. If recall is more important, TE should not be used. We also point out that TE encodes only the day of the year since the dataset has a relatively small spatial extension ($10.2^{\circ}\text{Lon}\times 8^{\circ}\text{Lat}$). In case of larger datasets at continental scale, a consideration of the spatial location for the encoding would also become relevant as biogeographical regions occur [74], which are characterized by different climate variabilities and anthropogenic drivers over time. Finally, we plot in Fig. 10 the loss (4) curve during training.