UAV Hyperspectral Mapping of Soil Organic Carbon: A Quality-Controlled Workflow from Field Sampling to Spatial Prediction

1. Introduction

Soil organic carbon (SOC) is a key indicator of soil fertility, agricultural productivity and ecological sustainability [1,2]. It influences nutrient cycling, soil aggregation, water-holding capacity and cation exchange, while also representing an important terrestrial carbon pool [3,4,5]. Spatially explicit information on SOC is therefore relevant for precision agriculture, soil management and carbon-sequestration strategies. Conventional laboratory analysis provides reliable SOC measurements, but collecting samples at the density required for field-scale mapping is costly and time-consuming. Proximal and remote sensing approaches have consequently attracted increasing interest as complementary tools for estimating soil properties over larger areas.

Hyperspectral remote sensing is particularly promising because it records reflectance in many narrow and contiguous spectral bands, allowing subtle variations related to soil composition, mineralogy, moisture, texture and surface condition to be observed [6,7,8,9]. While modern hyperspectral sensors capture vast amounts of data, there is often significant redundancy among contiguous bands, and a carefully selected subset of wavelengths can be sufficient for most soil-property retrievals. In the visible and near-infrared (VNIR) range, SOC-related spectral responses are commonly expressed through broad and overlapping absorption patterns associated with organic matter, soil colour, mineral composition and other soil characteristics that co-vary with carbon content. Compared with multispectral sensors, hyperspectral instruments can preserve more of this information and may therefore support more detailed soil-property retrieval [10]. UAV-borne hyperspectral systems are especially attractive for field-scale applications because they can acquire imagery at centimetre-scale spatial resolution under flexible acquisition conditions. This makes them suitable for heterogeneous and fragmented agricultural fields, where satellite observations may contain mixed pixels (i.e., individual pixels comprising a combination of multiple land cover or surface types) or be unavailable during suitable bare-soil conditions [11].

UAV hyperspectral imagery has been used to estimate a range of soil properties, including salinity, moisture, pH, texture, nitrogen, organic matter and trace elements. However direct UAV hyperspectral studies of soil organic matter (SOM) and SOC remain comparatively limited and reported performance varies substantially among sites. Yan et al. [12] demonstrated field-scale SOM mapping using a small calibration dataset and averaged centimetre-scale spectra within a $1\times 1$ m support around each sampling location. This approach recognised that a laboratory soil sample is more appropriately represented by a local image neighbourhood than by an individual pixel. Song and Zhang [13] further demonstrated the potential of UAV hyperspectral imagery for SOC mapping across multiple farmland areas, while showing that retrieval accuracy depends strongly on field conditions, spectral preprocessing and modelling choices. Similar sensitivities have been reported for UAV-based estimation of soil texture [14], salinity [15], available copper [16] and heavy metals [17]. Collectively, these studies show that UAV hyperspectral data can support detailed soil-property mapping, but also indicate that image quality, the way spectra are extracted around soil-sampling locations, and validation design are decisive for the credibility of the resulting models and maps.

The reliability of a predictive soil mapping model depends largely on the quality of the spectra extracted at the sampling locations, which imposes challenges with UAV based data acquisition and processing. Due to variable environmental conditions and imperfections in imaging hardware, the UAV hyperspectral data may contain radiometric inconsistencies, geometric misalignment, striping and mosaicking artefacts that introduce spectral variation unrelated to soil properties. If such effects are not identified and mitigated before modelling, learning algorithms may interpret them as meaningful soil variability. Radiometric calibration, geometric processing and systematic image-quality assessment are therefore prerequisites for physically meaningful soil-property prediction rather than preprocessing steps [18].

The confidence to a model outcome depends on how the model performance is evaluated. When multiple pixels are extracted around the same soil-sampling location, these observations are spatially related and cannot be treated as independent samples. Cross-validation must therefore respect the hierarchical structure of the data by keeping all pixel observations associated with a given soil sample within the same validation fold [19]. Evaluation across independent fields is additionally needed when the objective is to assess transferability beyond the calibration setting.

Against this background, this study evaluates a quality-controlled workflow for field-scale prediction of SOC from UAV hyperspectral imagery. It treats image preprocessing, the extraction of spectra around soil-sampling locations, and model validation as interdependent parts of a single workflow rather than isolated technical steps. First, reflectance is radiometrically calibrated, harmonised against colocated proximal spectroradiometer measurements and screened at pixel level before spectra are extracted around the soil-sampling locations. Second, three strategies for relating laboratory observations to high-resolution imagery are compared: single-pixel spectra, neighbourhood-aggregated features and local pixel-level augmentation. Third, all image observations associated with the same soil sample are assigned to the same cross-validation fold to prevent optimistic estimates arising from spatially related pixels. The aim is therefore not only to compare predictive accuracy, but to establish a reproducible and statistically sound procedure for translating sparse soil observations into trustworthy UAV-based field-scale predictions.

2. Materials and Methods

2.1. Overall Workflow

Figure 1 summarises the proposed workflow for quality-controlled UAV hyperspectral soil-property mapping.

2.2. Study Areas, Soil Sampling and Laboratory Analysis

Two agricultural regions with contrasting soil-forming conditions and climate were selected to test the workflow across independent settings.

2.2.1. Flanders Site (Belgium)

The first study area was located in East Flanders, Belgium, in an arable landscape in the municipality Merelbeke-Melle (approximately 10km south-east of Ghent, site center near $50.{98}^{\circ }$ N, $3.{78}^{\circ }$ E; UTM zone 31N, EPSG:32631). The region has a temperate maritime climate. Eight sandy-loam bare-soil agricultural fields were surveyed (B2, Caritas, L1, L3, M5B, R6, RvK1A–B and S24), covering a range of field sizes and organic-carbon levels (Table 1).

A total of 85 georeferenced topsoil samples were collected across the eight fields, with between 3 and 27 samples per field depending on field area and expected within-field variability. At each location, a soil sample was taken from the 0–10 cm layer; each sample was a composite of 10 cores within a diameter of approximately 2m. Geo-coordinates of the samples were recorded with a RTK GNSS receiver in WGS84 and subsequently projected to UTM zone 31N to align with the hyperspectral orthomosaics.

2.2.2. Central Macedonia Site (Greece)

The second study area was located in the Region of Central Macedonia, northern Greece, with agricultural fields distributed across the Regional Units of Imathia and Thessaloniki. The area is characterized by a Mediterranean climate with warm, dry summers and cool, wetter winters and comprises a heterogeneous agricultural landscape with a wide range of crop types. The surveyed fields included both annual crops and perennial orchards, predominantly olive and pomegranate orchards, as well as wheat and cotton fields. This diversity provided contrasting field conditions, vegetation structures, and soil-management practices across the study area.

A total of 42 georeferenced topsoil samples were collected across the selected agricultural fields. At each sampling location, soil was collected from the 0–10 cm layer by combining multiple soil subsamples taken within an approximately 2 m radius around the sampling point to obtain a representative composite sample.

2.2.3. Laboratory Analysis

Soil analyses were performed at separate laboratories for the two campaigns, each following its routine protocol. The SOC values of the Flemish samples were analysed at the accredited ILVO soil analysis laboratory by dry combustion following ISO 10694, and they were reported as a mass fraction (% of oven-dry soil). The Central Macedonia samples were air-dried prior to analysis and their SOC values were determined by Walkley–Black wet oxidation. These samples were additionally characterized for particle-size distribution (sand, silt and clay) by the Bouyoucos hydrometer method, for coarse fraction by sieving, for pH by potentiometric method, for electrical conductivity (EC) by conductometric method, for calcium carbonate (CaCO3) by the Bernard calcimeter method, and for total nitrogen (N) by the Kjeldahl method.

Across the Flanders site, measured SOC ranged from 0.86% to 2.14% (mean $1.41\%$, median $1.40\%$, standard deviation $0.34\%$; coefficient of variation $23.8\%$). Per-field summary statistics are provided in Table 1. Field-level mean SOC ranged from $0.93\%$ in S24 to $1.75\%$ in Caritas and M5B.

At the Central Macedonia site, SOC ranged from 0.60% to 1.30%, with a mean of $0.87\%$ and a standard deviation of $0.19\%$. Per-field SOC statistics are reported in Table 2, with field-level means ranging from $0.76\%$ in PM to $1.03\%$ in RD. These values indicate a generally lower and narrower SOC range than that observed at the Flanders site.

For the additional soil properties measured in Central Macedonia, sand content ranged from 47 to 59% (mean $52.93\%$, standard deviation $2.98\%$), silt content from 23 to 36% (mean $28.69\%$, standard deviation $3.22\%$), and clay content from 16 to 21% (mean $18.38\%$, standard deviation $1.36\%$). The coarse fraction ranged from 2.0 to 2.7% (mean $2.34\%$, standard deviation $0.21\%$), while CaCO₃ ranged from 0.8 to 7.0% (mean $3.89\%$, standard deviation $1.85\%$). Soil pH varied from 6.98 to 7.61 (mean $7.32$, standard deviation $0.19$), and EC from 408 to 604 $\mu$S cm⁻¹ (mean $512\mu$S cm⁻¹, standard deviation $50\mu$S cm⁻¹). Total N ranged from 0.30 to 1.10 g kg⁻¹, with a mean of $0.66$ g kg⁻¹ and a standard deviation of $0.23$ g kg⁻¹.

These laboratory values constitute the reference target for model calibration and validation.

2.3. UAV Hyperspectral Acquisition

The spectral imagery was acquired using a tuneable Fabry–Pérot interferometer (FPI) hyperspectral camera developed by VTT. The imager is able record wavelengths from 650 nm to 1200 nm and allows the user to programmatically select the spectral bands to be captured. The imager scans the wavelengths of interest and produces a 2D frame for each spectral band. The images produced by the 1.3 megapixel InGaAs sensor were binned by the camera once horizontally and once vertically to reduce the storage size of the recorded data. The instrument was mounted on a multi-rotor UAV and included a global navigation satellite system (GNSS) receiver for georeferencing the image data. Table 3 shows the key specifications of the imager.

Six bands with center wavelengths of 670, 730, 800, 930, 970 and 1100 nm in the red to near-infrared (VNIR) region were selected for this study. The objective was to obtain the spectral features most relevant to soil organic carbon while limiting acquisition time and data volume. The UAV was operated at 80–100 m altitude above ground level under fair-sky conditions, giving a ground sample distance of approximately 1.1–1.4 cm. Three reference panels of known reflectance (50%, 80% and 95%) were placed within each scene and imaged together with the field; their use for radiometric calibration is described in Section 2.4.

Flights were conducted over the Flanders fields between 19 September and 21 October 2025 and over the Central Macedonia fields between 16 and 27 February 2026. The two campaigns consequently differ in season: the Flemish fields were surveyed post-harvest under bare-soil conditions, whereas the Central Macedonia survey fell within the wet winter period, when autumn-sown cereals were actively growing and orchards carried a permanent canopy. The implications for spectral quality and retrieval performance are considered in Section 4. Per-flight dates and the resulting orthomosaics are available in the archived dataset [20,21]. Acquisitions were generally scheduled as close to local solar noon as practical (13:00-14:00 local time) to minimise variation in illumination geometry, although flight times varied across some fields due to cloud conditions and the limited availability of suitable bare-soil conditions. Because the FPI camera acquires the wavelength bands sequentially rather than simultaneously, a spatial offset arises between bands when the platform is moving; this was addressed by low flight speed. Individual frames were mosaicked into per-field orthomosaics using OpenDroneMap [22].

2.4. Quality Control and Spectral Extraction

The reliability of a predictive model depends on the quality of the spectra acquired, so the spectral extraction step was treated as a modelling decision rather than as preprocessing. This section describes, in the order applied, how reflectance was calibrated and harmonised across sensors, which pixels were considered valid, how vegetation was handled without sacrificing sample size, and how the pixels were related to each laboratory measurement.

2.4.1. Radiometric calibration

The reflectances were calibrated using the signal from the three reference panels of nominal reflectances 0.50, 0.80 and 0.95 imaged within each scene. The mean signal over a manually delineated region on each panel was fit with with zero-offset constrained regression to the nominal reflectance, so that zero signal maps to zero reflectance and a single multiplicative gain per band was obtained. Constraining the intercept avoids the negative reflectances that an unconstrained fit can produce at low signal, at the cost of absorbing any additive path-radiance term into the gain.

For the Central Macedonia campaign, field reference spectra were collected using a PSR+ 3500 spectroradiometer (Spectral Evolution Inc., Lawrence, MA, USA), covering the 350–2500 nm range via a 512-element silicon photodiode array (350–1000 nm; 2.8 nm spectral resolution) and two 256-element InGaAs arrays (970–1910 nm, 8 nm resolution; 1900–2500 nm, 6 nm resolution). Measurements were acquired with a bare fibre under ambient illumination, in nadir orientation matching the UAV viewing geometry, and calibrated against a Spectralon reference panel; PSR acquisitions were performed concurrently with the UAV flights to minimise illumination differences between the two datasets. PSR spectra were resampled to the six UAV band centres and bandwidths by averaging the high-resolution PSR reflectance over each band`s nominal bandwidth prior to comparison with UAV reflectance.

An additional harmonisation step was applied for the Central Macedonia campaign because residual differences remained between panel-calibrated UAV reflectance and colocated PSR measurements. UAV spectra were matched to PSR observations by searching a limited neighbourhood around each PSR coordinate to account for geolocation and orthomosaic misalignment, retaining only spectrally consistent and locally homogeneous candidates. A multiplicative gain was then estimated independently for each band from the ratio of PSR to UAV reflectance and applied uniformly to the corresponding orthomosaic. This step provides an independent, field-based check on the panel calibration specifically at the Central Macedonia site, where residual bias was detected; no comparable PSR deployment was available at the Flanders site, where the panel-based calibration alone was used throughout.

2.4.2. Pixel-Level Validity Criteria

The pixels within each orthomosaic were checked against the following validity criteria. Pixels flagged as no-data were excluded. Any band returning non-positive reflectance was treated as invalid. Due to factors such as sensor noise, vignetting, or data transmission errors, individual spectral bands may occasionally be missing for a given pixel. Because a partial spectrum cannot be used consistently across the derived features, a pixel was discarded entirely if any of the six bands was missing.

2.4.3. Neighbourhood-Based Extraction of UAV Pixels

To link each laboratory measurement with the high-resolution UAV imagery, a circular neighbourhood of radius r was defined around each recorded soil-sampling coordinate This approach avoids relying exclusively on a single nearest pixel, which may be affected by local surface heterogeneity, geolocation uncertainty or image noise. Figure 2 shows representative neighbourhoods over the UAV orthomosaic. Three radii were evaluated, $r\in \{0.25,0.5,1.0\}$ m, corresponding to circular supports of 0.5, 1.0 and 2.0 m diameter.

2.4.4. Vegetation Filtering

The Central Macedonia campaign included perennial orchards and standing annual crops, and therefore fully bare spectral support could not be assumed. NDVI was calculated for each valid pixel as

$$\mathrm{NDVI}={\displaystyle \frac{R_{800}-R_{670}}{R_{800}+R_{670}}},$$

where $R_{800}$ and $R_{670}$ denote the reflectance measured in the bands centred at 800 and 670 nm, respectively.

Rather than applying a fixed NDVI threshold, which can remove different numbers of pixels or even entire sampling locations depending on local vegetation conditions, vegetation influence was reduced using relative NDVI ranking within each sampling neighbourhood. Candidate pixels were ranked in ascending order of NDVI, with distance from the sampling coordinate used to resolve ties. A fixed fraction q of the pixels with the lowest NDVI values was then retained (Figure 3). In this study, $q=0.6$ was used.

This procedure preferentially retains the least vegetated spectral observations available around each soil sample while avoiding the loss of sampling locations that may occur with an absolute NDVI threshold. No additional vegetation selection was required for the Flanders campaign because imagery was acquired under predominantly bare-soil conditions.

2.4.5. Spectral Extraction Strategies

From this common neighbourhood definition, three strategies for relating the laboratory measurement to the imagery were compared.

• Single-pixel spectra (baseline): The reflectance of the single pixel containing the sampling coordinate was used directly. This represented the conventional point-based extraction and served as the baseline against which the neighbourhood strategies are assessed.
• Neighbourhood-aggregated features: For each sampling location, the per-band mean and standard deviation of reflectance were computed over all selected pixels in the circular neighbourhood, summarising both the average spectral response and its local variability. This yielded a single feature vector, and hence one training observation, per soil sample.
• Local pixel-level augmentation: Every selected pixel within the neighbourhood was treated as an individual training observation, each inheriting the laboratory SOC value of its parent sampling location. This increased the number of pixel-level training instances available to the models while retaining the fine-scale spectral variation within each support. To prevent information leakage, every pixel was tagged with the identifier of its parent sampling location so that all pixels from the same location could be constrained to the same cross-validation fold.
• Derived spectral features: In all strategies, the same set of derived spectral features was appended to the reflectance values: first-order spectral derivatives between consecutive bands, and normalised differences between every pair of bands (${\mathrm{ND}}_{ij}=(R_j-R_i)/(R_j+R_i)$), analogous to normalised vegetation indices but applied across the full six-band set. This yielded 26 features per observation for the single-pixel and pixel-augmentation strategies (six reflectances, five derivatives and fifteen normalised differences), and 32 features for the aggregated strategy, which additionally included the six per-band standard deviations. For the aggregated case, the derived indices were computed from the per-band mean reflectances. The feature definition was identical at both study sites.

2.5. Model Development and Validation

Because multiple pixels can originate from the same soil sample, all observations associated with a given sampling location were kept together within the same cross-validation fold. This prevented spatially related pixels from appearing in both training and validation data. Model performance was evaluated using grouped cross-validation.

Four regression algorithms were evaluated: Random Forest (RF) [23], gradient-boosted trees (XGBoost) [24], a one-dimensional convolutional neural network (CNN) [25], and TabPFN [26]. RF and XGBoost were included as established tree-based ensemble methods, while the CNN provided a conventional deep-learning baseline. TabPFN was additionally evaluated as a pretrained transformer-based foundation model for tabular prediction.

The RF model used 200 trees with unlimited depth. The XGBoost model used 200 trees with a learning rate of 0.1, a maximum depth of 4, and subsampling of both observations and features at 0.8. The CNN comprised two one-dimensional convolutional layers (16 and 8 filters, kernel size 3, Leaky-ReLU activation, $L_2$ weight regularisation of $4\times {10}^{-4}$), followed by two fully connected layers (32 and 16 units) and a single linear output. It was trained by minimising the mean squared error with the Adam optimiser (initial learning rate $2\times {10}^{-3}$) for up to 200 epochs with a batch size of 64. Hyperparameters were fixed at these values for all models and were not tuned per field or per strategy, to keep the comparison between extraction strategies unbiased.

Each model was finally retrained using all available observations for each strategy and algorithm to generate wall-to-wall SOC maps for the surveyed fields. All models used fixed random seeds to ensure reproducibility.

2.5.1. Validation Design

Model evaluation was organised at the soil-sampling-location level, matching the independent laboratory reference measurements. Performance was assessed using five-fold cross-validation with a common fold structure across extraction strategies, so that all strategies were evaluated on identical held-out sampling units. For pixel-level augmentation, all pixels from a given sampling location were assigned to the same fold, and their held-out predictions were averaged to the sampling-location level before any performance metric was computed. Predictive performance was quantified using RMSE and $R^2$ on the pooled out-of-fold predictions.

2.5.2. Statistical Significance Testing

The purpose of the significance testing was to determine, for each soil property and regression model, whether replacing single-pixel extraction with a spatially-informed extraction strategy produced a reduction in prediction error that reflects a genuine effect of incorporating spatial context, rather than variability attributable to the specific cross-validation split.

Each extraction configuration (a strategy–radius combination) was therefore compared against the single-pixel baseline of the same regression model on the same paired sample-level out-of-fold predictions, using an identical procedure applied uniformly across the full design: this comparison is the basic unit of inference, repeated once per property, model, strategy and radius. The null hypothesis, common to every such comparison, is that the baseline and method predictions are exchangeable with respect to the true value at each sampling location. This was tested with a paired permutation test on the RMSE difference, which avoids assuming a parametric error distribution and accommodates RMSE as a non-linear aggregate statistic; 95% confidence intervals for the RMSE difference were obtained by paired bootstrap resampling of the sampling locations.

Every configuration for a given property and model was tested against the same baseline. These comparisons therefore form a single family of hypotheses per property–model combination. Raw p-values were adjusted within each such family using the Holm step-down procedure [27] at $\alpha =0.05$, so that the reported significance controls the family-wise probability of at least one false positive among the configurations compared for that model.

3. Results

3.1. Comparison of Spectral Extraction Strategies

Table 4 summarises SOC prediction performance at the Flanders site across the three spectral extraction strategies, neighbourhood radii and regression models. RMSE reduction is expressed relative to the corresponding single-pixel baseline for each model, while the reported Holm-adjusted p-values quantify the significance of the paired improvement over that baseline.

The single-pixel strategy provided moderate predictive performance for RF, XGBoost and TabPFN, with $R^2$ values between 0.467 and 0.498 and RMSE values between 0.236 and 0.244 % SOC. In contrast, the CNN performed worse than other models when trained from a single pixel per sampling location ($R^2=-0.157$, RMSE $=0.359$ % SOC), indicating greater sensitivity to the limited point-based spectral representation.

Neighbourhood aggregation improved SOC prediction for most combinations of radius and regression model. For RF, XGBoost and TabPFN, performance generally increased as the neighbourhood radius increased. The best overall configuration was obtained with TabPFN using aggregated features at a radius of 1.0 m, reaching $R^2=0.784$ and RMSE $=0.155$ % SOC. This corresponded to a 36.5% reduction in RMSE relative to the TabPFN single-pixel baseline ($p_{\mathrm{Holm}}=0.001$). RF and XGBoost showed similar improvements at 1.0 m, reaching $R^2=0.708$ and $0.754$, respectively. The effect of aggregation was less consistent for the CNN, with performance improving at 0.25 and 0.50 m but deteriorating at 1.0 m.

Pixel-level augmentation produced a more consistent improvement across models. All augmentation configurations reduced RMSE relative to their corresponding single-pixel baselines and remained significant after Holm correction. The strongest augmented results were obtained at the larger neighbourhood sizes: at 1.0 m, RF, XGBoost and TabPFN reached $R^2$ values of 0.754, 0.749 and 0.763, respectively, with RMSE reductions of approximately 29–33%. The CNN benefited particularly from augmentation, reaching its best performance at 0.50 m ($R^2=0.540$, RMSE $=0.226$ % SOC), corresponding to a 37.0% reduction in RMSE relative to its single-pixel baseline.

Overall, 23 of the 24 neighbourhood-based configurations reduced RMSE relative to the corresponding single-pixel model, and 20 remained significant after Holm correction. Although no single extraction strategy was uniformly superior across all models and radii, the results consistently demonstrate the benefit of incorporating spatial context around the soil-sampling location. Larger neighbourhoods generally provided the strongest performance, although this relationship was model dependent and not strictly monotonic.

In Central Macedonia, SOC prediction was substantially more difficult than in Flanders. Because any combination of the used extraction strategy, used radius, and used algorithm did not produce a reliable model, Table 5 reports the best-performing configuration per strategy, selected by lowest pooled out-of-fold RMSE. Neighbourhood augmentation produced a numerical improvement, with the best configuration reaching $R^2=0.106$ and RMSE $=0.177$ % SOC, close to the site standard deviation of 0.19 %, but the improvement was not significant after Holm correction ($p_{\mathrm{Holm}}=0.501$). As even these best cases show no usable predictive skill, the remaining configurations do not affect any conclusion drawn here.

3.2. Generalization to Additional Soil Properties in Central Macedonia

The same neighbourhood-based framework was subsequently evaluated for the additional soil properties available in Central Macedonia. The clearest improvements were observed for EC and CaCO₃, for which pixel-level augmentation at a 1.0 m radius produced substantial reductions in prediction error. For EC, the best configuration was obtained with RF ($R^2=0.686$, RMSE $=28.13$), corresponding to a 38.5% reduction in RMSE relative to the closest-pixel baseline ($p_{\mathrm{Holm}}=0.001$). For CaCO₃, TabPFN with augmentation at 1.0 m achieved $R^2=0.625$ and RMSE $=1.119$, representing a 29.2% reduction in RMSE ($p_{\mathrm{Holm}}=0.0036$).

For pH, the lowest prediction error was obtained with CNN augmentation at 1.0 m ($R^2=0.609$, RMSE $=0.115$). Although this corresponded to a large 91.3% reduction relative to the same-model closest-pixel baseline, the magnitude of the relative improvement should be interpreted cautiously because the CNN baseline performed poorly. Prediction also remained weak for clay and total N, with the best configurations yielding $R^2$ values close to or below zero.

Table 6. Best-performing neighbourhood-based configuration for each soil property. RMSE improvement is calculated relative to the corresponding closest-pixel baseline using the same regression model.

Property	Model	Strategy	Radius (m)	${\mathit{R}}^2$	RMSE	RMSE improvement (%)	${\mathit{p}}_{\mathbf{Holm}}$
CaCO3	TabPFN	Augmentation	1.00	0.625	1.119	29.2	0.0036
Clay	XGBoost	Augmentation	0.50	0.042	1.308	7.8	1.0000
EC	RF	Augmentation	1.00	0.686	28.13	38.5	0.0010
N	TabPFN	Aggregation	0.25	-0.069	0.235	-0.4	0.8710
pH	CNN	Augmentation	1.00	0.609	0.115	91.3	0.00024

Table 6. Best-performing neighbourhood-based configuration for each soil property. RMSE improvement is calculated relative to the corresponding closest-pixel baseline using the same regression model.

Property	Model	Strategy	Radius (m)	${\mathit{R}}^2$	RMSE	RMSE improvement (%)	${\mathit{p}}_{\mathbf{Holm}}$
CaCO3	TabPFN	Augmentation	1.00	0.625	1.119	29.2	0.0036
Clay	XGBoost	Augmentation	0.50	0.042	1.308	7.8	1.0000
EC	RF	Augmentation	1.00	0.686	28.13	38.5	0.0010
N	TabPFN	Aggregation	0.25	-0.069	0.235	-0.4	0.8710
pH	CNN	Augmentation	1.00	0.609	0.115	91.3	0.00024

4. Discussion

The results show that the way UAV spectra are associated with laboratory measurements can substantially influence soil property prediction. Both the neighbourhood aggregation and the pixel-level augmentation improved the SOC prediction results in the favourable bare soil conditions in Flanders. Most neighbourhood-based configurations performed significantly better than their corresponding single-pixel baselines after Holm correction, supporting the overall benefit of incorporating spatial context. There was not however a consistently superior strategy across all models and different neighborhoud sizes. This indicates that the main finding is not optimizing a universal model but the value of incorporating spatial context surrounding each location.

This benefit can largely be explained by the difference in the area that the imagery and the laboratory measurements represent. Each laboratory value originates out of a composite sample collected from several subsamples within a local area, whereas a single UAV pixel represents only a small portion of that area. Neighbourhood-based extraction made the spectral representation of the imagery more consistent with the soil measurements. Aggregation may additionally reduce sensitivity to positional errors, micro scale variability and residual radiometric or mosaicking artifacts. The general improved performance at larger neighbourhood sizes supports this interpretation although the relationship was not strictly monotonic and remained dependent on the model and the extraction strategy.

Aggregation and augmentation incorporated spatial context in different ways. Aggregation reduced local spectral variability into a single representation for each sampling location, while the inclusion of per-band standard deviations retained some information about neighbourhood heterogeneity. Augmentation instead preserved the individual pixel spectra and allowed the models to learn from their local variation. The positive results obtained through aggregation also indicate that increasing the number of training observations is not always preferable to construct a more stable representation of each soil sample. Thus the aggregation offers practical benefit in terms of training time, memory requirements, and energy efficiency. Augmentation may nevertheless be particularly valuable for future large-scale deep-learning approaches capable of exploiting the increased number and diversity of pixel-level training instances.

The considerably weaker SOC results in Central Macedonia suggest that the extraction strategies could not overcome the limitations affecting SOC prediction at this site. This dataset was smaller and covered a narrower SOC range (Table 5) providing the models with less independent information and less target variability in order to learn a stable relationship. The hyperspectral acquisitions being in winter introduced further difficulties. Several fields were humid and contained growing vegetation, while the olive orchards had permanent canopies. Soil moisture can strongly modify soil reflectance and potentially mask SOC related spectral response while vegetation reduces the availability of representative exposed soil pixels. These factors probably contributed to the contrasting results with the optimal post harvest bare soil conditions in Flanders.

The relative NDVI-ranking procedure provided a consistent way of identifying the least vegetation-affected pixels under heterogeneous field conditions. Spectral-index thresholding can be highly sensitive to the selected threshold [28], and no single value could be applied across all samples without excluding an excessive number of pixels in some cases. Relative ranking was therefore less dependent on site-specific reflectance conditions. However, in neighbourhoods already dominated by bare soil, it may discard some valid bare-soil pixels. This represents a trade-off between retaining all potentially valid pixels and applying a consistent selection strategy less sensitive to local heterogeneous conditions.

These limitations did not affect all soil properties equally. Neighbourhood-based modelling improved the prediction of EC, CaCO₃, and pH, which indicates that local spatial context captured useful soil-related information even under challenging acquisition conditions. This could further indicate that these properties exhibited stronger spatial variation than SOC or stronger associations with the measured spectra than SOC. While properties such as pH and EC are commonly estimated through associations with mineral composition, moisture and other soil characteristics [29,30], several important mineral and carbonate absorption features occur at longer wavelengths than those used in this study, which may explain why clay content and total nitrogen remained difficult to predict. Therefore, the successful predictions likely reflect a combination of indirect spectral information and spatial covariation within the surveyed fields.

Radiometric calibration, cross-sensor harmonisation and pixel screening are preconditions for the reflectance values to carry physical meaning and their role is to establish the conditions under which the spectral support is reliable. The validation design and use of fixed model configurations enabled a consistent comparison of the extraction strategies which was the main objective of this study. All pixels associated with one location were kept within the same cross-validation fold, and their predictions were combined at each point, preventing information leakage from pixels sharing the same laboratory measurement. Since the strategies were evaluated at the same test locations, local factors such as soil conditions and illumination were also shared, focusing the comparison on the effect of spectral extraction. It must be emphasized that the performed evaluation does not provide any evidence for model generalization or transferability to other geographical areas. To achieve this, a different test setup is needed; this, however, is outside the scope of this study. Per-field sample counts (Table 1) were in any case too small to test the models reliably on unseen fields. Finally, part of the performance improvement observed for the neighbourhood-based strategies likely reflects spatial coherence within a field, such as local homogeneity in soil condition, alongside a physical relationship between reflectance and SOC.

The findings should be interpreted within the experimental scope of the study. The number of independent soil samples was limited, particularly in Central Macedonia, and augmentation increased the number of pixel-level training instances rather than the number of laboratory measurements. Although these six bands were deliberately selected to target primary SOC-related features, this restricted the spectral information available for modelling some soil properties. Because acquiring more bands increases data volume and acquisition time, this band configuration was a deliberate trade-off to maximize flight endurance (flight time per battery) for a UAV-borne FPI system operating within the VNIR range; extending the wavelength range into the SWIR, where several mineral and carbonate absorption features are more pronounced, would require a different sensor class and is outside the spectral range of the instrument evaluated here. Future work should extend the workflow through additional campaigns covering a wider range of soils, crops, seasons and surface conditions. Data-fusion strategies linking centimetre-resolution UAV observations with satellite imagery could support the upscaling of field-level predictions to regional applications. Federated learning could further support collaborative modelling across campaigns, institutions and geographical regions while reducing the need to centralise raw field data or precise sampling coordinates [31].

Taken together, these findings show that the spatial relationship between UAV imagery and field measurements is a central component of soil-property modelling. Neighbourhood-based representations improved upon single-pixel extraction under favourable bare-soil conditions and captured useful information for several properties under more challenging conditions. Spatial context should therefore be considered as an integral part of the modelling workflow along radiometric calibration, pixel screening and statistically sound validation.

5. Conclusions

This study developed and evaluated a quality-controlled workflow for converting UAV imagery (acquired using six selected bands of a tunable FPI hyperspectral sensor) and laboratory measurements into field-scale predictions of SOC and related soil properties. The workflow combined radiometric calibration and harmonisation, pixel screening, alternative spectral extraction strategies and statistically sound validation. In Flanders, both neighbourhood aggregation and pixel-level augmentation improved SOC prediction over single-pixel extraction across multiple regression models, achieving an R² of up to 0.784. In Central Macedonia, none of the evaluated configurations produced a reliable SOC model, although neighbourhood-based extraction improved the prediction of EC, CaCO₃ and pH, while clay and total nitrogen remained difficult to predict.

These results demonstrate that incorporating spatial context can improve the model performance when using of centimetre-resolution UAV imagery, but it cannot compensate if the imagery or reference dataset do not contain sufficient information about the target property. Reliable UAV-based soil mapping therefore depends on the complete workflow rather than on regression-model selection alone, including image quality, the correspondence between spectral extraction and field sampling, and validation that preserves the independence of sampling locations. The proposed workflow provides a reproducible basis for field-scale soil-property mapping with tunable UAV hyperspectral sensors and for further evaluation across additional campaigns and broader spectral ranges.