July 01, 2026
Monitoring canopy height change is essential for understanding carbon sinks and forest dynamics. Remote sensing enables consistent, large-scale observations of such changes, increasingly integrated with deep learning architectures such as Geospatial Foundation Models (GFMs). However, existing methods and datasets frame the problem as binary change detection, which overlooks both the continuous nature of change, especially for vegetation, and the inherent uncertainty in labels. We present the Canopy Height Change (CHC) dataset, providing 3 m resolution continuous canopy height differences and associated spatially resolved uncertainties across 10598 km2 of northern and western Spain. The dataset is paired with a co-located time series of PlanetScope satellite imagery. Based on the dataset, we introduce the task of uncertainty-aware change regression, associated metrics and strategies for fine-tuning GFMs. Furthermore, we evaluate state-of-the-art GFMs and highlight promising directions and remaining challenges for advancing continuous canopy height change estimation.
Understanding how canopy height changes over time is essential for quantifying carbon sinks and sources, ecosystem dynamics, monitoring disturbances, and detecting tree growth or deforestation. Recent advances in remote sensing, particularly in high resolution satellite imagery and deep learning, now enable consistent characterization of canopy height as wall-to-wall maps at continental and global scales [1]–[4]. Typically, these maps are produced through computer vision models, trained in a supervised approach on satellite imagery and target data originating from Airborne Laser Scanning (ALS) or spaceborne LiDAR, such as GEDI, over a wide range of geographies and biomes.
Recently, increasing attention has been directed toward estimating the temporal dynamics of canopy height, leading to the development of geospatial models that incorporate a time dimension [3], [5], [6]. However, the validation data of the resulting products is not public, and either limited to single-time-step evaluation or change evaluation without quantification of uncertainty, which may pose challenges, as the signal of height increase is typically small. Other works usually simplify the problem by reducing change to a binary predicate (e.g., loss/no loss [7]), ignoring the continuous nature of height changes where both incremental growth and partial canopy loss are essential variables. Hence, robust and large-scale benchmark datasets capturing continuous height change remain scarce, constraining efforts to evaluate the accuracy of height models and remote sensing-based datasets. To our knowledge, no existing public benchmark dataset provides continuous canopy height change with change uncertainty over a large region and at high spatial resolution.
Here, we present the Canopy Height Change (CHC) dataset (fig. 1), a benchmark dataset for continuous canopy height change at 3 m resolution on an area of 10598 km2 in northern and western Spain. The changes, including continuous values on loss and gain, are derived from two national ALS campaigns from 2018 and 2023 [8], [9] and associated with a simulated spatially resolved change uncertainty, using ALS acquisition geometry and sampling density to model systematic and statistical error propagation. We paired the change data with co-located time series of PlanetScope optical imagery between 2018 and 2023. This enables models to harness spatial and temporal patterns in the imagery to predict canopy height change.
We apply the CHC dataset to benchmark Geospatial Foundation Models (GFMs). Large-scale, pretrained GFMs have increasingly been adopted in deep learning and Earth observation [10]. In contrast to traditional remote sensing models, typically optimized for a single product and sensor, GFMs are designed to encode broad, semantically rich feature representations of land surface processes, which can then be fine-tuned with comparatively small datasets representing the target domain, such as canopy height [11]. The models are typically trained in a self-supervised approach on a wide range of multi-sensor Earth observation data, such as from optical and radar at various resolutions. Because it requires temporal reasoning, cross sensor transfer, and fine-scale spatial sensitivity, the CHC dataset represents a relevant use case for evaluating GFMs. The dataset is compatible with the PANGAEA benchmark framework [12] and can be accessed here: https://sid.erda.dk/sharelink/eP4ENGhKTv.
The main contributions of this study are as follows:
A new CHC dataset, covering 10598 km2 with 3 m resolution labels of continuous height change with uncertainty quantification, and PlanetScope satellite imagery time series across 6 years.
Accurate change uncertainty assessment during training and evaluation through quantified systematic and statistical data uncertainties based on simulation at pixel level.
Performance analysis of state-of-the-art GFMs on the CHC dataset.
A large number of remote sensing change datasets exist for the benchmarking of change detection models in a variety of application domains. Typically, they combine manual change labels with co-located bi-temporal remote sensing imagery, between which the change is assessed. The change detection problem is usually classified into Binary Change Detection (BCD), Semantic Change Detection (SCD), and Remote Sensing Image Change Captioning (RSICC). Application domains span urban settings (construction and demolition) (, [13]–[16]), land cover/land use change (, [17]–[19]), and disaster assessments (, [20]–[22]).
3D change detection adds an additional dimension, for example height, depth, or volume changes, where classic planimetric approaches might not be sufficient. The data is often sourced from ALS, Digital Elevation Models (DEMs), 3D models, or stereo/multi-view images and either represented as 3D difference (, Euclidean), or as “2.5D”, when changes are projected onto a plane (, difference between DEMs) [23]. Change detection from repeat ALS is used in a wide range of applications from coastlines to canopy structure [24]. However, public datasets that pair remote sensing imagery with ALS-based height change targets are rare (tab. 1).
In the context of urban areas, 3D changes are often DSM or point cloud-based [25]. The 3DCD dataset [26], [27] comprises DSM differences of artificial objects, such as from construction or demolition in a semantic change detection dataset of 18.8 km2 area, combined with aerial imagery. Also in the urban domain, the SMARS dataset [28] provides synthetically generated DSM changes and optical imagery from simulated 3D scenes.
| Domain | ||||||||
| imgs. | ||||||||
| (km\(^2\)) | Sensor | |||||||
| (m) | Task | Uncertainty | ||||||
| 3DCD [26] | Urban | 2 | 18.8 | Aerial | 0.5 | Continuous | \(\times\) | |
| SMARS [28] | Urban | 2 | 12 | Synthetic | 0.3–0.5 | Ternary | \(\times\) | |
| OpenCanopy-\(\Delta\) [7] | Natural | 2 | 166 | SPOT 6–7 | 1.5 | Binary | \(\times\) | |
| CHC (ours) | Natural | 6 | 10.6k | PlanetScope | 3.7 | Continuous | ✔ |
Open-Canopy [7] is an ALS-based tree height dataset at 1.5 m spatial resolution, paired with SPOT 6-7 imagery. A subset of the data provides binary change between 2022 and 2023 on an area of 166 km2. However, only significant height reductions of more than 15 m and of areas larger than 200 m2 are included in the dataset. To our knowledge, no metric vegetation height dataset that continuously quantifies tree growth has been published so far.
Discrete airborne LiDAR (Aerial Laser Scanning, ALS) is widely used to generate 3D representations of entire landscapes. Mounted on an aircraft, the LiDAR scanner can generate dense wall-to-wall point clouds with high accuracy, which in turn are used to create Digital Surface Models (DSMs) and Digital Terrain Models (DTMs) through rasterization [29]. Uncertainties in LiDAR measurements can originate from several sources: positional errors include the horizontal and vertical uncertainty of points induced by uncertainty in the GPS and inertial measurement units and the angular and range accuracy of the sensor. Uncertainties in surface representation are created by the discrete sampling process of rough and complex surfaces, leading to occlusion and undersampling. Finally, the point cloud classification (, into vegetation, buildings, infrastructure, etc.) introduces uncertainty depending on the classification algorithms applied [30], [31].
In a multi-temporal setting, it is usually assumed that changes can be detected despite varying acquisition parameters, or that the data is consistent across collections [32]. However, repeated ALS acquisitions typically come from different sensors and flight configurations. Variation in pulse density, driven by flight and sensor parameters, affects both DTM and DSM accuracy, and vertical and horizontal uncertainties must be evaluated across swaths [33]. Multivariate analyses have shown that scan geometry is a key driver of LiDAR height uncertainty and bias, with measured height, local variability of neighboring cells, and pulse density among the strongest predictors [34]. Additional influences include beam divergence, peak pulse power, canopy structure, terrain slope, and canopy edges [32], [35]. ALS systematically underestimates canopy height because the highest crown points are often missed [34], [36]–[40], which is exacerbated by lower point densities [41], [42]. Furthermore, systematic biases have been attributed to DSM algorithm choices, pulse penetration characteristics, and understory structure [43]. Species‑related traits such as crown shape also influence under/overestimation, and reduced pulse density (, through point cloud thinning) exacerbates these species‑dependent biases [44].
Bias mitigation in LiDAR‑derived canopy heights typically relies on statistical correction models for varying acquisition conditions. Several approaches have been tested, although their explanatory power often remains modest. For example, [34] evaluated linear and non‑linear formulations to correct height underestimation and found limited improvements, reflecting the probabilistic nature of the discrete sampling. [41] proposed to correct density-induced biases through modeling the bias-density relationship by thinning high-density samples to lower densities to capture a parameterized relation, noting diminishing biases above roughly 7 ( points per square meter ) A probabilistic framework has also been introduced by [38] , who modeled the stand ‑ level vertical forest structure and simulated the sampling process as a function of pulse density and footprint size.
Geospatial foundation models (GFMs) have become ubiquitous in recent years, challenging the established paradigm of fully-supervised learning for remote sensing tasks [10]. Powered by advances in self-supervised learning, they serve the purpose of generalistic feature extractors that can be applied to any modality and for a variety of downstream tasks. The models are typically trained on large, diverse remote sensing datasets through contrastive learning [45], generative learning [46], [47], or supervised learning [48]. However, their limitations remain poorly understood due to their recency and bias towards classification tasks for which labels are more readily available [49]. On the PANGAEA benchmark [12], which only includes two regression tasks (BioMassters [50] and Open-Canopy [7]), it has been noted that regression remains challenging, with supervised baselines still outperforming complex self-supervised GFMs.
Heteroscedastic regression refers to supervised learning settings in which target noise varies across data points [51]. Although deep neural networks can tolerate a high degree of label noise when provided with a sufficiently large set of clean labels [52], [53], they also exhibit high capacity in fitting random noise, which can degrade performance under uncertain labels [54], [55]. Various approaches have been explored to mitigate the impact of noisy labels [56]. For regression, re-weighting the mean squared error (MSE) loss by target uncertainty can outperform a standard MSE by reducing the contribution of noisy training samples [51], [57].
The CHC dataset is designed to assess a model’s capability to regress the metric change in canopy height, both decrease and increase, and hence its ability to detect tree removal and growth. In this section, we describe the underlying data, the processing, and the characteristics of the final benchmark dataset.
We made use of the publicly available point cloud data from the Spanish National Plan of Air Orthophotography (PNOA), conducted by the Insituto Geográfico National (IGN) [8], [9]. The ALS data was collected in different nation-wide campaigns, out of which Cobertura 2a (2015-2021) and 3a (2022-2026) fall into our observation period, as defined by the availability of the PlanetScope data. The measurements are organized in regional lots with varying flight dates, sensors, and pulse densities (fig. 2), depending on region and time.
We selected the areas that had double LiDAR coverage from both campaigns, leaving us with measurements from 2018 (2a) and 2023 (3a) between which we calculated change. Furthermore, we only selected areas covered during leaf-on period between May and October, to avoid cross-seasonal effects originating from the absence of leaves in one of the campaigns. The filtering resulted in an area of 10598 km2, covering the region Cantabria in the north of the country, and parts of Extremadura at the western border to Portugal. We used processed point clouds in format LAS v1.4, as delivered by CNIG [8], [9] as input for the subsequent processing.
We paired the data with PlanetScope optical imagery from the same regions. The imagery is available with approximately 3 m resolution at daily frequency, and RGB and Near-Infrared bands. We selected cloud-free scenes within a seasonal window based on the MODIS phenology product [59] in late summertime, as done by [2]. If no coverage was reached at that time, we expanded the window until at least one image matched the criteria.
We rasterized the LiDAR point clouds to digital surface models (DSM) by calculating the 95th height percentile on a regular \(3\times 3 \si{\meter}\) grid for each of the two years. The points were filtered based on the encoded ASPRS classes [60], and only the vegetation-related classes (, classes 3-5) were included. Through vertical differencing of the DSMs, we computed the surface height change \(\Delta h^{(95)}\). Compared to calculating the difference between canopy height models (CHMs), which represent the height above ground, this approach has the advantage of avoiding additional errors introduced by DTM measurements and interpolation. We masked areas where vegetation was below 3 m in both years, to avoid impacts of non-woody vegetation changes, field harvests. We calculated the explanatory features on the same grid, using 2D binned statistics (bin size was equal to GSD, , 3 m).
DSM measurements contain both systematic and statistical uncertainties, which hinder the direct interpretation of observed height differences. Systematic uncertainty \(\sigma_\mathrm{syst}\) typically originates from technical constraints and flight parameters, whereas statistical uncertainty \(\sigma_\mathrm{stat}\) arises from the probabilistic process of measuring height at discrete locations through the LiDAR pulses and aggregating them into grid cells.
We estimated uncertainties and offsets between the campaigns through simulation and by analyzing the residuals between ALS campaigns at static locations. We estimated \(\sigma_\mathrm{syst}\) from the residuals between the two ALS campaigns on flat road surfaces that we assumed remained unchanged throughout the observation period, based on the Spanish Land Cover and Land Use Information System (SIOSE) (IGN, 2015). We obtained a distribution of residuals with standard deviations between 0.17 m and 0.23 m, and offsets between −0.09 m and −0.01 m, depending on the combination of ALS campaigns.
In contrast, \(\sigma_\mathrm{stat}\) is heteroscedastic and originates from the discrete LiDAR point measurements of complex geometric surfaces, and the aggregation of their 95th height percentile (\(h^{\left(95\right)}\)) into grid cells. We assumed that grid cells \(i\) with high pulse density (\(\rho_i\geq \SI{10}{\square\ppm}\), which equals 90 pulses in the \(3 \times 3\) grid cells) represent the surface geometry sufficiently well and used their heights (\(h_i^{\left(95\right)}\left(\rho_i^{\left(10\right)}\right)\)) as reference values. We then simulated the sampling process with \(j\) lower pulse densities (\(\SI{1}{\square\ppm}<\rho_i^{\left(j\right)}<\SI{10}{\square\ppm}\)) and calculated the corresponding height difference \(d_i^{\left(j\right)}=h_i^{\left(95\right)}\left(\rho_i^{\left(j\right)}\right)-h_i^{\left(95\right)}\left(\rho_i^{\left(10\right)}\right)\), along with other explanatory features \(\mathbf{v}_i\left(\rho_i^{\left(j\right)}\right)\) that represent surface geometry and ALS acquisition geometry (tab. S2). We trained a multi-layer perceptron (MLP) adjustment model \(f(\rho, \mathbf{v})\) on 7.935 × 105 training grid cells to predict height offset \(\hat{b}_i\) and variance \(\hat{\sigma}_i^2\) based on pulse density and explanatory features. Using the Gaussian negative log likelihood (NLL) loss (eq. 1 ), the model is incentivized to represent every prediction as a Gaussian distribution conditional on the input data, allowing the model to attenuate the effect of uncertain targets [61]. \[\begin{align} \mathcal{L}_{\mathrm{NLL}}=\frac{1}{N}\sum_{i=1}^{N}\left(\frac{1}{2}\exp{\left(-\log{\hat{\sigma}_i^2}\right)}\left(b_i-\hat{b}_i\right)^2+\frac{1}{2}\log{\hat{\sigma}_i^2}\right) \label{eq:nll} \end{align}\tag{1}\] Finally, we combined systematic and statistical uncertainties through variances summation for all pixels in both campaigns. To illustrate the effect of the adjustment model, we built a building-specific model using static buildings to represent unchanged complex geometries. The model attained an Expected Calibration Error (ECE) of 0.21 m for the predicted residual standard deviation, though it generally underestimated uncertainty relative to the observed residuals. Further details on uncertainty and offset estimation are provided in the supplementary materials.
The dataset is organized in 1123 tiles of \(1024\times 1024\) with 3 m resolution, aligned with the PlanetScope imagery. PlanetScope is available as a time series between 2018 and 2023 with one to nine unordered images per year, allowing for averaging or temporal shift augmentation. The height-change tiles represent changes between 2018 and 2023 and have three data layers: height change (in m), maximum canopy height (height above ground) in both years (in m), and change variance (in m2) (fig. 3). We apply a threshold of min. 3 m height to separate trees from other vegetation. 27 % of the area is covered with vegetation higher than 3 m in at least one of the years. The dataset covers a variety of land cover types, with Grassland, Croplands, and Forests being the most prominent ones (fig. 4d). The canopy height shows a long tail distribution and is on average slightly lower in the train set compared to validation and test set. The height changes at vegetated locations are centered around zero, steeply sloping on both sides of the distribution (fig. 4a). The relative changes show a similar pattern, while a local peak at −100 % represents complete tree removal. These changes represent raw changes and might be obscured by non-significant variations in canopy height, which can occur, , at canopy edges.
We observe different median canopy heights in the various land cover classes represented in the dataset. Forest and natural woodlands show the highest canopy (median 10.4 m) and orchards the lowest (median 3.5 m), just above the minimum height threshold (fig. 4c). Note that non-forest land cover can still include scattered trees that do not meet the criteria for forests. The median changes are positive in all land cover classes, ranging from median increases of 0.2 m in silviopastures to 0.7 m in urban green areas, which, however, might be confounded by misclassified buildings in the LiDAR point clouds. The largest median height change uncertainty is found for forests and natural woodlands (1.1 m), the lowest for orchards (0.8 m).
The dataset is split (fig. 2) into training (592 tiles), validation (262 tiles) and test set (269 tiles) by region in concentric squares with the validation set acting as a buffer between train and test set [58]. The test split only contains low-uncertainty change values, we applied a vegetation mask and retained only pixels with an absolute change \(z\)-score above 1.65, corresponding to cases where the observed deviation would occur with less than 10 % probability under a no-change scenario.
We used the CHC dataset to assess the performance of GFMs on the task of canopy height change regression. We employed pretrained GFM encoders and finetuned only the decoder on the training split. Model performance was assessed on the heldout test split, examining several approaches for integrating label uncertainty during training.
Scale-MAE [46] is a scale‑aware autoencoder that learns multiscale geospatial features using a ViT encoder with Ground Sample Distance positional encoding and a Laplacian‑pyramid decoder. It is pretrained via masked image modeling on multiscale remote‑sensing data. RemoteCLIP [45] adapts CLIP-style dual encoders for remote sensing by jointly training image and text encoders using largescale image-caption data, including UAV imagery. This enables learning semantically aligned visual-text representations suitable for zeroshot classification and other downstream tasks. DOFA [47] is designed to handle inputs from diverse remote sensing sensors. It uses a wavelength conditioned dynamic hypernetwork to generate modality adaptive parameters, allowing for arbitrary input channel counts. DOFA is pretrained on five heterogeneous modalities for cross sensor generalization and robust performance on unseen modalities. Prithvi-EO 2.0 [63] is a temporal ViT foundation model pretrained using a masked autoencoder objective on Harmonized Landsat–Sentinel2 (HLS) data. It uses inter-patch spatial attention and intra-patch temporal attention, enabling spatiotemporal modeling of satellite time series. CROMA [64] uses contrastive learning with masked autoencoding to exploit spatially aligned optical and radar imagery from Sentinel-1 and 2. It integrates separate unimodal encoders for optical and radar inputs (contrastive learning) and a joint multimodal encoder (reconstruction loss). SatlasNet [48] is trained on a large collection of remote‑sensing datasets spanning segmentation, regression, detection, and classification tasks. It processes time‑series inputs with Swin‑Transformer backbones and fuses features via temporal max‑pooling at multiple scales. SpectralGPT [65] uses masked autoencoding to learn visual representations in a progressive manner, using various Sentinel-2-based data sources and a ViT-based backbone architecture. The fine-tuning includes downstream tasks such as image classification, segmentation, and change detection. TerraMind [66] uses modality specific ViT tokenizers to embed data from various sources, such as Sentinel-1 and 2, land cover maps, or image captions. During pretraining, the model learns to reconstruct missing modalities through a masked modeling approach. DINOv3 [67] is a vision transformer pretrained through self-distillation, producing dense features that have been shown to generalize well across remote sensing tasks. ECHOSAT [3] is a global CHM time series dataset at 10 m resolution, generated with a spatiotemporal transformer model that performs pixel-wise regression of canopy height from optical and radar imagery, using a growth-constrained loss and GEDI-based canopy height to enforce physically realistic forest dynamics over time.
We evaluated the performance of the GFMs using the PANGAEA benchmark framework [12]. Each model received time series of annual PlanetScope imagery as input and predicted canopy height change in a dense prediction task. We used the frozen pre-trained GFM encoders to create embeddings at different depths of the encoder, which we then assembled into feature pyramids. These were passed on to a trainable UPerNet decoder [68] which performed multi-scale aggregation and produced dense predictions of canopy height change. We upscaled the output features with bilinear interpolation and added a residual refinement block at the end.
The encoders were natively pre-trained following different objectives and data regimes. Hence, not all encoders were compatible in terms of required input channels or time-series capacities. PANGAEA aligns input channels with those used to train the GFM, zero‑padding any unmatched channels [12]. For encoders without native support of time series input, PANGAEA introduces an intermediate Lightweight Temporal Attention Encoder (LTAE) [69] after the encoder. Each time step is embedded individually and then passed to the LTAE, which applies attention with temporal position cues to fuse these features across time steps. The fused feature pyramids are then passed to the UPerNet decoder. \[\begin{align} \mathrm{wMSE}=\frac{\sum_{i=1}^{N}\left[\frac{1}{\sigma_i^2}\left(\hat{h}_i-h_i\right)^2\right]}{\sum_{i=1}^{N}\frac{1}{\sigma_i^2}} \label{eq:wmse} \end{align}\tag{2}\] The training of the learnable UPerNet decoder parameters was performed for each encoder separately over 50 epochs on the training split of the CHC dataset, with AdamW optimizer, learning rate warmup, and cosine scheduler. We applied a vegetation mask to the dataset to only train and validate at locations with vegetation higher than 3 m. The mask was eroded by one pixel to remove confounding effects on vegetation edges.
We compared two loss‑weighting schemes to assess how GFMs behave with and without pixel‑level uncertainty. The weighted MSE (wMSE) [51] (eq. 2 ) scales errors by the inverse target variance, down‑weighting uncertain height‑change estimates so the model focuses on more reliable signals. The thresholded MSE (tMSE) instead computes the loss only for pixels with a z‑score above 1.65, excluding low‑signal‑to‑noise samples. Finally, the standard MSE treats all targets as equally reliable. We then evaluated the models on the test split of the CHC dataset, using root mean squared error (RMSE), mean absolute error (MAE), normalized MAE (nMAE), and coefficient of determination (\(R^2\)). In addition to overall metrics, we separately masked increases and decreases to assess performance on positive and negative height changes.
In addition to the GFMs, we included a UNet [70] and a ResNet-50 [71] with UPerNet decoder as baseline models that were trained fully supervised without any frozen parameters. We also included a randomly initialized ResNet-50 encoder that was not updated during training, as a null model. Training and evaluation followed the same logic as for the GFMs. Furthermore, the ECHOSAT dataset [5] acted as an independent baseline trained globally baselineon GEDI-derived canopy height.
The performance of the evaluated baselines and GFMs on the CHC test split showed that none of the GFMs could outperform a supervised UNet trained with the thresholded MSE (tMSE) loss, which achieved an RMSE of 5.65 m (tab. 2). Prithvi and DOFA achieved the best overall results among the GFMs, with an RMSE of 5.90 m and 5.99 m, respectively. On the other hand, Scale-MAE was inferior to the baseline of a randomly initialized and non-fine-tuned ResNet-50 encoder across all loss weighting schemes. The ECHOSAT dataset outperformed both GFMs and supervised baselines with an RMSE of 5.25 m. The results suggest that the fine-grained, subpixel structural cues driving continuous height change remain challenging for all of the frozen encoders when used for regression on PlanetScope. Moreover, task-specific models, such as the ECHOSAT Temporal-Swin-UNet trained on GEDI-derived canopy height, show superior performance, despite the coarser spatial resolution of 10 m. Current remote sensing GFMs often struggle to outperform supervised baselines and specialized models on pixel level tasks [12], [72], [73].
| Encoder | MSE | Increase | Decrease | All | |||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 5-8 (lr)9-12 (lr)13-16 | s | t | w | RMSE \(\downarrow\) | MAE \(\downarrow\) | nMAE \(\downarrow\) | \(R^2\) \(\uparrow\) | RMSE \(\downarrow\) | MAE \(\downarrow\) | nMAE \(\downarrow\) | \(R^2\) \(\uparrow\) | RMSE \(\downarrow\) | MAE \(\downarrow\) | nMAE \(\downarrow\) | \(R^2\) \(\uparrow\) |
| ResNet-50 | x | 3.51 | 2.57 | 0.68 | -1.05 | 11.47 | 9.05 | 0.94 | -1.09 | 6.98 | 4.55 | 0.82 | 0.21 | ||
| UNet | x | 4.11 | 2.95 | 0.78 | -1.81 | 8.13 | 6.03 | 0.63 | -0.05 | 5.65 | 3.89 | 0.70 | 0.48 | ||
| ECHOSAT | 3.81 | 2.83 | 0.75 | -1.43 | 7.55 | 5.68 | 0.59 | 0.09 | 5.25 | 3.70 | 0.67 | 0.55 | |||
| CROMA | x | 3.22 | 2.33 | 0.62 | -0.73 | 10.94 | 8.91 | 0.93 | -0.90 | 6.61 | 4.34 | 0.78 | 0.29 | ||
| DINOv3 | x | 4.16 | 3.19 | 0.85 | -1.89 | 8.76 | 6.49 | 0.68 | -0.22 | 5.95 | 4.20 | 0.76 | 0.42 | ||
| DOFA | x | 4.53 | 3.17 | 0.84 | -2.42 | 8.42 | 6.35 | 0.66 | -0.13 | 5.99 | 4.14 | 0.75 | 0.41 | ||
| Prithvi | x | 4.27 | 3.44 | 0.91 | -2.05 | 8.52 | 6.34 | 0.66 | -0.15 | 5.90 | 4.32 | 0.78 | 0.43 | ||
| RemoteCLIP | x | 5.39 | 4.02 | 1.06 | -3.84 | 10.00 | 7.57 | 0.79 | -0.59 | 7.12 | 5.10 | 0.92 | 0.17 | ||
| SatlasNet | x | 4.97 | 3.66 | 0.97 | -3.12 | 8.90 | 6.83 | 0.71 | -0.26 | 6.43 | 4.63 | 0.83 | 0.33 | ||
| Scale-MAE | x | 4.58 | 3.87 | 1.02 | -2.49 | 12.38 | 9.51 | 0.99 | -1.44 | 7.83 | 5.59 | 1.01 | 0.00 | ||
| SpectralGPT+ | x | 4.59 | 3.80 | 1.01 | -2.51 | 11.11 | 8.61 | 0.90 | -0.96 | 7.23 | 5.27 | 0.95 | 0.15 | ||
| TerraMind | x | 4.40 | 3.63 | 0.96 | -2.23 | 11.15 | 8.59 | 0.90 | -0.97 | 7.17 | 5.14 | 0.93 | 0.16 | ||
| Encoder | MSE | Increase | Decrease | ||||||
|---|---|---|---|---|---|---|---|---|---|
| 5-7 (lr)8-10 | s | t | w | Prec. | Recall | F1 | Prec. | Recall | F1 |
| ResNet-50 | x | 0.79 | 0.71 | 0.75 | 0.46 | 0.57 | 0.51 | ||
| ResNet-50 | x | 0.76 | 0.94 | 0.84 | 0.70 | 0.34 | 0.45 | ||
| UNet | x | 0.85 | 0.80 | 0.83 | 0.60 | 0.69 | 0.64 | ||
| ECHOSAT | 0.83 | 0.85 | 0.84 | 0.63 | 0.59 | 0.61 | |||
| CROMA | x | 0.77 | 0.96 | 0.85 | 0.79 | 0.33 | 0.47 | ||
| DINOv3 | x | 0.83 | 0.76 | 0.80 | 0.55 | 0.65 | 0.59 | ||
| DOFA | x | 0.83 | 0.77 | 0.79 | 0.54 | 0.63 | 0.58 | ||
| Prithvi | x | 0.81 | 0.83 | 0.82 | 0.60 | 0.57 | 0.58 | ||
| Prithvi | x | 0.84 | 0.79 | 0.81 | 0.57 | 0.65 | 0.61 | ||
| RemoteCLIP | x | 0.69 | 0.41 | 0.51 | 0.30 | 0.59 | 0.40 | ||
| RemoteCLIP | x | 0.77 | 0.63 | 0.69 | 0.40 | 0.57 | 0.47 | ||
| SatlasNet | x | 0.69 | 0.99 | 0.82 | 0.28 | 0.01 | 0.02 | ||
| SatlasNet | x | 0.82 | 0.74 | 0.78 | 0.51 | 0.63 | 0.57 | ||
| Scale-MAE | x | 0.69 | 0.00 | 0.00 | 0.31 | 1.00 | 0.47 | ||
| Scale-MAE | x | 0.62 | 0.07 | 0.12 | 0.30 | 0.91 | 0.45 | ||
| SpectralGPT+ | x | 0.79 | 0.76 | 0.77 | 0.49 | 0.53 | 0.51 | ||
| TerraMind | x | 0.78 | 0.88 | 0.82 | 0.61 | 0.42 | 0.50 | ||
| TerraMind | x | 0.76 | 0.94 | 0.84 | 0.69 | 0.33 | 0.44 | ||
Furthermore, we investigated the models’ performance in spatially detecting the direction of change on the test set (tab. 3). For detecting height increases, CROMA outperformed the supervised baselines with an F1 score of 0.85, while a supervised UNet achieved highest performance for decreasing heights (\(\mathrm{F1} = 0.64\)). Several GFMs could not outperform the randomly initialized and non-fine-tuned ResNet-50 encoder in detecting decrease, such as CROMA, RemoteCLIP, and Scale-MAE.
We investigated if providing target uncertainty to the models at training time by using wMSE or tMSE improves their performance on the test split. The results show that the best-performing models benefit from uncertainty-based weighting or thresholding. When trained with the tMSE loss, the RMSE of DOFA and the baseline UNet decreased by 4 % and 9 %, respectively, compared to the same models trained with the MSE loss. Training Prithvi using wMSE decreased the RMSE by 8 % on the test split. GFMs with lower overall performance show a less clear response to uncertainty weighting.
We presented the CHC canopy height change dataset, to our knowledge the first public benchmark dataset for continuous canopy height change regression with associated uncertainty. We benchmarked popular Geospatial Foundation Models (GFMs) on the task of height change regression and reported an overall low performance compared to a supervised UNet baseline and the task-specific ECHOSAT dataset. All models had difficulties regressing the change magnitude, while to some degree being able to estimate the direction of change. This indicates a challenging task, that we hope will foster research to address the notable gap in performance compared to classification tasks.
Notably, we showed that uncertainty estimates may play an important role, by (1) eliminating low-quality test samples in a fundamentally noisy task, and (2) providing an adjustment variable to tune the quality/quantity balance in the fine-tuning sets for GFMs. Our work provides essential data for measuring real canopy height change, and supports future research in environmental monitoring with remote sensing.
RF and MB acknowledge funding from the Danish National Research Foundation, Center for Remote Sensing and Deep Learning of Global Tree Resources (TreeSense), DNRF192. We thank Planet Labs PBC for the provision of the PlanetScope imagery.