Introduction

Based on the widely accepted topoclimatological framework established by Oke (1987), local-scale climate analysis focuses on atmospheric variations occurring over distances ranging from several hundred meters to tens of kilometres. This spatial scale is particularly relevant for studying climate variability within a vineyard or small appellation area. Climate variability at the local level can lead to substantial differences in water availability for vines during grape ripening due to the highly localised nature of storm events (Bois et al., 2020). This variability is also reflected in considerable temperature differences (Bonnefoy et al., 2010), leading to consequences for plant species diversity and spatial distribution in natural environments (Oldfather et al., 2016; Opedal et al., 2015) and for crop phenology as observed for grapevines (de Rességuier et al., 2020). Temperature variations often correlate with terrain characteristics (elevation, slope, exposure, topographical depressions) but also with other topographical features (Bois, 2024; Oke, 1987), such as land use (forest or lake) and, more rarely, with surface soil characteristics (Jacquet & Morlat, 1997). As topography remains relatively stable over time, a few years of studies are often sufficient to understand the spatial variability of temperatures at the local level.

Thermal zoning of wine-growing regions at the local level can therefore be achieved through the integration of topographical information with temperature measurements taken by sensors deployed across the study area (de Rességuier et al., 2020; Le Roux et al., 2017; Webb et al., 2016) . However, no formal approach to the climate measurement system implementation (installation sites, network density, and sensor positioning) has been proposed to depict climate variability at the local scale level. The World Meteorological Organization rightly recommends deploying measurement systems in open areas, avoiding unusual topographical conditions such as steep slopes (World Meteorological Organization, 2024). At the local scale, characterising climate variability requires sampling a range of contrasted situations in order to represent the diversity induced by landscape elements, as for example the impact of hedges (Guyot & Seguin, 1978). Yet the concept of spatial resolution and sampling density in environmental sciences is a widely documented issue (Gruijter et al., 2010; Hengl, 2006). A few papers assess the change in weather station density to climate prediction at macro (Njoku et al., 2023) or at mesoscale levels (Berndt & Haberlandt, 2018; Rosillon et al., 2024). While significant improvement in temperature spatial interpolation accuracy have been observed with considerable variation of sensor density over Sweden (Njoku et al., 2023), little benefits have been reported in Belgium (Rosillon et al., 2024) or Northern Germany (Berndt & Haberlandt, 2018).

The deployment of climate measurement networks should be guided by prior knowledge of the sources driving spatial variability in the variable of interest.

Precipitation sampling strongly relies on the type of precipitation event, the catchment size and the time-resolution targeted, assessed through geostatistical approaches (Berne et al., 2004; Lebel & Laborde, 1988). For temperature, sensor network deployment can be determined using an a priori terrain analysis allowing accurate sampling of the different topographical configurations over the area under study (Webb et al., 2016). This said, a practical question remains: how many sensors should be deployed to relevantly depict the climate at the local scale level. Tomislav Hengl (2006), in a paper that addresses the issues of map scale and resolution in a very practical and relevant way, offers some answers. Amongst other suggestions, he points out that sample density and raster resolution can be related by:

$p=0.0791\sqrt{\frac{A}{N}}$     (Equation 1)

where p is the pixel resolution in meters, A is study area in m² and N is the number of samples (i.e., weather stations). Conversely, the number of station N can be determined by:

$N=\frac{0.006257A}{p²}$     (Equation 2)

For example, if the target resolution (pixel size) is 10 m over an area of 10 km² (i.e., 10⁷ m), a minimum of 25 stations should be deployed.

However, this approach is purely Euclidean-distance-based and therefore fails to account for environmental factors such as topographic features, whose spatial structure can now be depicted at very high spatial resolution (typically 1 m or 0.5 m at the time of writing) from LiDAR-derived terrain and surface models. Station network density may be substantially limited when the studied climate variable is highly predictable by topographic features, as for local solar radiation or air temperature (Bois, 2024; de Rességuier et al., 2020; Joly et al., 2012; Oke, 1987).

In this article, we explore the practical and methodological aspects of implementing high-resolution spatial temperature monitoring within a wine-growing context. Specifically, we address two critical research questions: (1) what is the quantifiable impact of sensor network density – comprising air temperature sensors housed within naturally ventilated radiation shields – on the spatial interpolation accuracy of daily minimum (T_min) and maximum (T_max) temperatures? and (2) to what extent does the choice of interpolation algorithm – ranging from multiple linear regressions to machine-learning ensembles and hybrid geostatistical models – depend on the density of the available data? These questions are addressed through a detailed case study involving the topoclimatic zoning of two neighbouring wine appellations in the Côte de Nuits region of Burgundy (France). By evaluating these factors at the grape-cluster level (60 cm) within a traditional high-density (1 m × 1 m) vineyard architecture, this work leverages a spatial intensive network of 112 sensors deployed across an 800-ha study area to provide preliminary insights into the operational thresholds required for reliable climate zoning in precision viticulture.

Materials and methods

1. Study area

The investigation was conducted within the viticultural landscapes of Chambolle-Musigny and Morey-Saint-Denis, located in the Côte de Nuits subregion of Bourgogne, France (Figure 1). While the core Protected Designation of Origin (Appellation d’Origine Controlée, AOC) areas of these villages encompass an area of 509 ha, the instrumented study domain extends over approximately 800 ha (8 km²). This experimental perimeter provides a comprehensive representation of the DPO “hierarchy” of Bourgogne terroirs, including 66 ha across six Grand Cru appellations (considered the top-level appellation), 105 ha of Premier Cru appellations, 163 ha of Village appellations, and 166 ha of Regional appellations. The study area is characterised by pronounced topographical heterogeneity, a primary driver of local climates. Elevation ranges from 233 to 408 m above sea level (mean: 270 m). Although the average slope is 6 %, the terrain exhibits significant steepness in the upper-slope sectors, where gradients reach up to 40 %. Notably, 11 % of the vineyard areas are located on slopes exceeding 15 %. They are particularly prone to variations in solar radiation interception and nocturnal cold air drainage.

2. Terrain analysis and descriptors

Terrain morphology was characterised using a High-Resolution Digital Terrain Model (DTM) provided by IGN (RGE Alti®) at a spatial resolution of 5 m. Preliminary analysis revealed that the discrete nature of the integer elevation values in the raw DTM induced step-like artifacts when calculating elevation-based topographical features (e.g., slope). To ensure surface continuity and mitigate these discontinuities, the elevation data were resampled using a Multilevel B-Spline interpolation (SAGA GIS software, v. 9.2.1; Conrad et al., 2015).

The refined DTM was processed within the SAGA GIS environment to derive a suite of 19 digital maps (rasters) based on 10 terrain descriptors: slope, slope orientation (transformed into sine and cosine), upstream watershed size, relative position index in relation to slope, topographic humidity index, landscape openness index, landscape closure index, potential solar irradiance from April to September, valley depth index, exposure index to the main winds in Burgundy (NNE, E, SSW, W, and NW). A detailed account of these descriptors and their calculation parameters is provided in Table S1. These topographical variables were used both to guide the sensor network strategic deployment and as predictive covariates for the subsequent spatial interpolation of the temperature data.

To account for macro-topographical features influence on temperature, a smoothed version of the DTM was generated by averaging elevation values within a 100 meter radius of all nearby pixels, following the windowing approach established by Joly et al. (2012). This smoothed elevation was used to derive the same 19 terrain descriptors, resulting in 38 rasters characterising the topographical environment of the study area.

3. Air temperature measurement network

Based on a strategic deployment design informed by preliminary topographical analysis and consultations with local winegrowers, a network of 112 sheltered air-temperature sensors was installed across the study area. This experimental array was progressively deployed between mid-January and late March 2023, achieving a high-density of approximately 112/8 = 14 sensors per km² (1 sensor per 7 ha).

The network’s operational capacity increased in phases: 64 sensors were active by March 1st, 96 by March 16th, and the full complement of 112 sensors was operational by August 2, 2023. Continuous data acquisition (every 10 minutes) was maintained through September 29, 2024.

The instrumentation consisted of thermistors probes (Model PB-5005-0M6, accuracy of 0.4 °C between 0 and 60 °C, 0.05 °C reading resolution or better, 10 sec. response time in water, calibrated following Gemini – manufacturer – protocol, detailed information at https://assets.geminidataloggers.com/pdfs/original/6404-tk-4023.pdf) coupled with Tinytalk data loggers (TK-4023, Gemini Data Loggers, UK, powered with ½AA 3.6V Lithium batteries which longevity is typically 2 years with this device), each housed within a naturally ventilated solar radiation shield (RS3-B -Hobo©/Onset©, USA). Tinytalk data loggers (including thermistor probes) and solar radiation shields have been purchased in 2021 at a price of 75 €/unit and 112 €/unit, respectively.

To ensure the microclimatic data representativeness and accuracy, a standardised installation protocol was strictly followed: (1) sensors were placed within active vineyard plots (2) and secured to wooden trellis stakes, 60 cm above the ground (grape cluster zone height) to target fruit height climate (3) sensors were intentionally installed adjacent to dead/missing vines, or young replant, in order to avoid microclimatic bias induce by dense foliage in close proximity to the solar radiation shield. Such conditions can generate localised cooling or increased humidity during the day and night, potentially altering summer temperature observations by 1 to 3 °C (Madelin et al., 2014). (4) Whenever possible, the solar radiation shields were oriented toward the west or northwest of the trellis stake to minimise potential radiative interference from the support structure.

Temperature data has been controlled by comparing graphically 10 minutes time-step data of groups of 20 sensors, day by day. Abnormal data (time shifts, inconsistent pattern of one sensor in comparison with the rest of the sensors, for several hours, pike of high/low temperature several degrees higher or lower than the group) was removed. Note also that data was lost from February 29th to March 6th, 2024 due to sensors memory saturation (wrong timing of data collection by operators). Finally, missing values or abnormal data detection and elimination resulted in gaps from 1.1 to 25.5 % of missing data by sensor. 17 of the 112 sensors exhibited more than 10 % of missing/abnormal data, while 93 sensors had 1.6 % of missing/abnormal data or less. Log of data is available on the data depository (see section 6 below).

Missing data, when limited to less than 3 consecutive values (3 × 10 minutes = 30 minutes) have been replaced by linear interpolation from previous and following values (na.approx function of the zoo package in R, see section 6 below).

Finally, minimum, maximum and average daily temperatures have been calculated. Minimum daily temperature of day d is lowest value recorded from 18:00 UTC of day d-1 to 18:00 UTC of day d. Maximum temperature of day d is the highest recorded from 06:00 UTC of day d to 06:00 UTC of day d+1. These time range followed Méteo-France (French national weather service) Standards.

Average temperature was calculated of the average values of temperature data collected during day d.

4. Spatial interpolation of air temperature

Daily average, minimum and maximum temperature from March 1st, 2023 to September 29th 2024 were interpolated. Seven spatial interpolation methods were compared using leave-one-out cross validation (LOOCV). LOOCV was selected arbitrarily. We acknowledge that k-fold cross validation would also have been an option because of the large number of sensors available.

The methods compared were: nearest neighbour (NN), multiple linear regression (LM), partial least squares regression (PLS), and Random Forest decision trees (RF). Three additional methods, LMK, PLSK, and RFK, consisted of adjusting the LM, PLS, and RF prediction by applying ordinary kriging of the residuals. LMK, also known as regression-kriging, is widely used in spatial interpolation of environmental data.

LM, RF and PLS methods used the 38 terrain descriptors as predicting variables, although only up to 3 terrain descriptors were selected for each interpolation, to avoid prediction model overfitting.

This preliminary selection was carried out by adjusting the multiple regression model as follows:

• 1 – All multiple regression models involving 1 to 3 of 38 terrain descriptors were evaluated, except for those including at least two autocorrelated descriptors. Two descriptors are considered autocorrelated if the Student's t-test applied to the Pearson correlation moment product is significant at the alpha = 0.01 % threshold.
• 2 – The multiple regression model that minimises the Akaike information criterion (AIC) was then selected.
• 3 – The terrain descriptors of the selected model were used as predictors for the LM, PLS, RF, LMK, PLSK, and RFK methods.

This selection procedure is repeated for each interpolation (i.e., for each day, each variable, and each data cross-validation step).

Initially, the interpolation methods performance was compared for each day of the study period and each variable (T_avg, T_min, T_max) with the data collected by all the sensors using leave-one-out cross-validation (LOOCV).

5. Bootstrap subsampling experiment

We tested the performance of seven spatial interpolation methods using a reduced number of temperature sensors. Through this experiment, we attempt to partially answer two questions:

• What is the optimal density of temperature sensors across the study area to produce a sufficiently robust climate spatialisation?
• Which spatial interpolation method is best suited to lower sensor densities in the study area?

Several subsamples were taken from the 112 available sensors:

• Sensor batch A: A validation set of 30 sensors sampled and removed from the data set
• From the remaining 92 sensors, a semi-random sensor batch B was drawn: samples of 20, 40, 60, or 80 sensors.

These four sets of sensors (20, 40, 60, or 80) make it possible to evaluate the interpolation performance using a smaller, and therefore more economical, set of equipment.

All sampling were achieved through a semi-stratified approach: the sampling was constrained by selecting a given number of sensors based on specific topographical characteristics. The terrain features of each sensor's environment were described by 38 terrain descriptors (see paragraph 2 of this section and Table S1).

Sensors with similar topographical characteristics were grouped using K-Means clustering method. Preliminary analysis of intragroup dispersion (within-group sum of squares) led to the selection of eight clusters (see Figure 1 and Figure S1), as recommended in Makles (2012). Within each group, a random sample was taken in proportion to the number of samples selected (20, 40, 60, or 80 sensors), as well as a sample of 30 sensors to be used for validation (validation set A). A constraining rule was applied so that at least one sensor from each topographic class was present in each of the two sensor sets, A and B.

Once these two sets of sensors (set A for validation and set B for interpolation) have been sampled, the daily temperature data measured by set B are predicted by interpolation at the coordinates of each of the 30 sensors in set A.

The experiment was carried out 30 times (30 random draws for each sample size) for average minimum and maximum temperature spatial interpolation of a set of 150 days randomly selected from the entire study period (March 2023 to September 2024). 30 repetitions were chosen as a balance between computational efficiency and statistical stability. Days for which the number of sensors deployed was less than the sample size required for the experiment were discarded.

6. Statistical analysis and software used

The performance of spatial interpolation was evaluated a suite of statistical indicators, including Root Mean Squared Error (RMSE), Bias, and the Coefficient of Determination (R²). Additionally, modelling efficiency (hereafter abbreviated as EFF) was calculated according to the criteria established by Nash and Sutcliffe (1970).

All data processing and geostatistical analyses were performed using the R software (R Core Team, 2024). Spatial data and raster management were handled via terra package (Hijmans, 2025), while terrain descriptors were generated using SAGA GIS modules interfaced through the RSAGA package (Brenning et al., 2025). General data manipulation and visualisation were performed within the tidyverse framework (Wickham et al., 2019), while specific modelling tasks were executed using specialised libraries as following: replacement of missing data used the function na.approx of the package zoo (Zeileis & Grothendieck, 2005); statistical efficiency was computed using hydroGOF (Zambrano-Bigiarini, 2024); geostatistical modelling, including automated variogram fitting and ordinary kriging, was implemented via gstat (Pebesma, 2004); and machine learning and regression-based interpolations were performed using randomForest (Liaw & Wiener, 2002) and pls (Liland et al., 2024) packages.

All source codes and source data are available at https://kdrive.infomaniak.com/app/share/909599/582126ad-94c5-43a2-a7d2-a60370ba3546

Apart from English translation or English correction in this manuscript by the DeepL translation engine (www.deepl.com/fr/translator, between September 2025 and April 2026), artificial intelligence was not used in this research.

Results

1. General climate features during the growing season

Table 1 compares the climatic conditions from April to September for the years 2023 and 2024 with the 1995–2024 normals for the Dijon station, located 12 km northeast of the study area.

The climate of the study area is classified as temperate (HI = 1,908 °C.d), sub-humid (DI = 110 mm), with very cool nights (CI = 11.2 °C) according to the geoviticultural climate classification by Tonietto and Carbonneau (2004).

The two years under study exhibit quite contrasting climatic conditions. The year 2023 recorded average temperatures from April to September that were 1.52 °C above normal. The Huglin index (2,225 HDD) classifies this vintage as “warm temperate” according to the classification mentioned above. Precipitation was 32 mm below normal for the April–September period.

The 2024 vintage is characterised by temperatures very slightly above normal. The distinctive feature of this vintage is high cumulative precipitation from April to September, classifying it as “humid” according to the geoviticultural climate classification (Tonietto & Carbonneau, 2004). It is the third year with the highest number of rainy days since the start of the station’s climate series (1922). More than one out of every two days (52 %) was rainy from April-September 2024 while the average (1995-2024) rate of rainy days is 41 %.

In the study area, the spatial distribution of temperatures (April–September average) reveals the following patterns (Figure 2).

Minimum temperatures are observed at low elevations (240 m a.s.l.) or at higher elevations (400 m a.s.l.). Low temperatures are also observed at the valley floor south of the study area in 2023 (Figure 2a). Minimum temperatures range over the study area 2.42 °C and 2.23 °C in 2023 and 2024, respectively.

Maximum temperatures (Figures 2b and 2e) range from 1.44 °C and 1.22 °C in 2023 and 2024, respectively, over the study area. The highest values are recorded in the middle and lower parts of the hillside, and the lowest values at higher elevations.

The lowest average temperatures (April–September) are recorded at higher elevations and in (Figure 2f) the valley floor south of the study area (Figure 2c). They range from 17.81 to 19.29 °C in 2023 (a difference of 1.49 °C) and from 16.87 to 18.09 °C in 2024 (a difference of 1.22 °C).

2. Performance of spatial interpolation methods according to the sample size

Cross-validation based on the 112 sensors reveals that the nearest neighbour (NN) method produces the largest errors (RMSE), while the regression-kriging (LMK) method minimises interpolation errors (Figure 3), both for minimum temperatures (RMSE = 0.40 °C on average) or maximum temperatures (RMSE = 0.55 °C on average). Random Forest with kriging (RFK) performance is similar to LMK for minimum temperature, and PLSK is slightly higher (0.42 °C). Residual kriging increases the accuracy of LM, PLS and RF methods for minimum temperature. In comparison, residual kriging has little impact on predicted RMSE of T_max. The RMSE of the interpolations never exceeds 1.60 °C, and very rarely exceeds 0.80 °C. Interpolation performs much better for minimum temperatures (median EFF = 0.660) than for maximum temperatures (median EFF = 0.264).

Spatial interpolation efficiency varied substantially from day to day, with a strong season pattern for maximum temperatures (Figure 4). From March 2023 to October 2024, efficiency of T_max interpolation using regression-kriging varied between 0 and 0.4, suggesting poor interpolation performances during the grapevine growing season. From November to February, T_max interpolation performance was considerably more variable with an average efficiency frequently over 0.4.

Interpolation errors are not spatially homogeneous. Figure 5 shows the bias (mean error, i.e., predicted values minus observed values) for spatial interpolation using regression kriging of 2023 and 2024 minimum and maximum temperatures. Relatively low biases are observed for both minimum (Figures 5a and 5b) and maximum temperatures (Figures 5c and 5d), with the exception of a few sensors often located at the edges of the interpolation area. In particular, the sensor located in the valley south of the study area shows a negative bias (underestimation) for minimum temperatures and a positive bias (overestimation) for maximum temperatures. Maximum temperatures are also significantly underestimated at the sensors located at the northernmost edge of the study area (Figures 5c and 5d). It is also noted that some sensors show marked positive biases in maximum temperatures in the southern part of the area.

The performance of the seven spatial interpolation methods was evaluated for samples of 20, 40, 60, and 80 sensors (Figure 6). Unsurprisingly, interpolation is of poor quality with small sample sizes, particularly for maximum temperatures. For minimum temperature, the median of all methods but RF and RFK remain below 0.4, the reference we used for acceptable performance from Perreault et al. (2022). The Random Forest (RF) method appears to be the most robust. It produces higher median efficiency levels than the other methods for both minimum and maximum temperatures, even for small sample sizes.

In general, there is a real increase in performance when the sample size increases from 20 to 40 sensors. Beyond 40 sensors, residual kriging improves the interpolation of minimum temperatures. However, residual kriging degrades the quality of the interpolation when small sample sizes are available. For maximum temperatures, increasing the number of data points gradually improves the accuracy of spatial interpolation.

Discussion

1. Interpolation methods performance

Among the methods compared, regression-kriging (LMK) produced the most accurate interpolation of daily minimum and maximum temperatures when conserving the entire sensor network set. This approach, equivalent to kriging with external drift (Hengl et al., 2003), offers explicit integration of relief parameters, enabling analysis of the contribution of the various terrain descriptors selected (Joly et al., 2012). The RF(K) methods have been observed to be amongst the most accurate for interpolation of coastal mud data (Li et al., 2011) or average monthly temperatures (Appelhans et al., 2015). Otgonbayar et al. (2019) observed considerable benefits of using the Random Forest algorithm with topographical descriptors and satellite-sensed surface temperature in comparison to PLS regression for interpolation of monthly temperature climatologies over Mongolia. However, depending on the situation, machine learning-based methods performed better (Le Roux et al., 2017; Sekulić et al., 2020) or worse (Cho et al., 2020; Webb et al., 2016) than multiple linear regression for the interpolation of daily temperatures or daily temperature-based agroclimatic indices at the local scale level.

The lower performance of spatial interpolation for maximum temperatures has already been observed in a study of a slightly larger area, in Bordeaux (Le Roux et al., 2017), with error magnitudes (RMSE) similar to those observed in the present study. However, a critical distinction lies in the measurement height: while Le Roux et al. (2017) and de Rességuier et al. (2020) used an intermediate height of 1.2 m, our sensors were positioned at 60 cm to directly monitor the grape-fruiting zone. This choice represents a significant departure from the World Meteorological Organization (WMO) guidelines, which mandate a standard height of 1.25 to 2.0 m to capture regional atmospheric signals. By descending to 60 cm, our measurements capture the full complexity of the canopy-ground boundary layer. At this height, the T_max signal is heavily "influenced" by localised energy fluxes, such as soil albedo and restricted air circulation within the vine rows. This explains why our T_max interpolation performance is lower than that of studies conducted at 2 m or 1.2 m: the "thermal noise" at cluster level is physically more intense. de Rességuier et al. (2020) observed a clear difference between the maximum temperatures measured inside the vineyard and those measured by a weather station installed in accordance with WMO standards. This reported bias was also seasonally dependent. Other studies show much better accuracy in spatial interpolation of daily maximum temperatures compared to minimum temperatures at the local scale. In the present study, the residuals from LM, PLS, or RF interpolation generally did not display any spatial structure that could have significantly improved the interpolation of maximum temperatures. In addition, substantial differences in maximum temperatures were observed over short distances. The quality of spatial interpolation of T_max shows marked seasonality: in winter, interpolation efficiency is significantly higher than in spring or summer (Figure 4b). It is likely that vegetation affects air temperature on a very local scale: despite the precaution of installing sensors in the vineyard where no vine grew (location of a dead or missing vine), it is possible that the vigour of the surrounding vines or the method of soil management affects maximum temperatures. In similar environments (sensors installed within the vineyard), Madelin et al. (2014) observed typical differences of 2 °C between two sensors installed under cover in summer, one outside the vegetation (missing vine location) and one located less than 50 cm away within the vine canopy. de Rességuier et al. (2023) found differences between both minimum and maximum daily temperatures measured over a cover-cropped zone compared to a tilled-soil zone within the same vineyard plot, with a marked seasonal effect. Finally, the seasonal differences could also be explained by a higher occurrence of synoptic-scale climatic situations leading to weaker relationships with topographic parameters.

2. Optimising sensor network for temperature zoning

The number of sensors to be deployed for climate mapping is an important consideration, as it is constrained by the costs and time available for installation and data management. At the local level, this study suggests using at least 40 temperature sensors to map temperatures locally (Figure 6). Increasing the number of sensors significantly improves the quality of maximum temperature interpolation and robustness of interpolation methods, resulting in lower dispersion in RMSE, as shown in Figure 7.

For minimum temperature, little improvement in spatial interpolation accuracy was reached with more than 60 stations. This corresponds to a density of 7.5 sensors/km² considering the surface of the study area (800 ha).

For maximum temperature, the recommended network deployment would be more than 80 sensors (10 sensors/km²). Yet, as mentioned here above, the poor accuracy of spatial interpolation of temperature is potentially due to the vegetation microclimate around the sensor. In such cases, no strong conclusion should be made for maximum temperature density issues.

If fewer than 80 sensors are installed, the Random Forest method(s) (+ kriging of residuals) RF(K) offer the best performance. The use of the SVM method, based on RF in machine learning, has shown good performance at the local level in the Bordeaux region (Le Roux et al., 2017). For large populations (80 sensors and more in this study), RFK and LMK offer similar performance. These conclusions, based on a single case, need to be tested in other wine-growing areas to draw more solid conclusions about the density of climate sensors to be deployed to capture climate variations.

The performance of spatial interpolation is affected by specific topographic features. For example, high errors (bias and RMSE) are observed at the sensor located at the bottom of a valley located south of the study area, for both minimum and maximum temperatures (Figure 5 and Figure S3). Thus, in areas with contrasting topography and land cover, it is likely more appropriate to increase the number of sensors. Conversely, in areas with relatively homogeneous topography, a lower sensor density will suffice to achieve similar spatial interpolation performance. The optimal pixel size of the maps produced, assessed using a simple method linking sensor density to map scale proposed by Hengl (2006) (see Eq. (1) in this paper), varies for the study area from 50 meters (20 sensors) to 25 meters (80 sensors) and up to 21 meters (for all 112 sensors deployed). In other words, the effective climate signal is about 25 m. Note that these dimensions are rather indications related to map scales rather than precise values that should be retained. They are indeed significantly higher than the resolution of the digital terrain model (5 meters) used to produce daily temperature rasters. This finer resolution makes it possible to consider microclimatic impacts caused by landscape features that would be poorly accounted for at the recommended resolution of about 25 m. For example, road and railway embankments can generate cold air lakes that are conducive to winter or spring frost (Bois et al., 2011; Quénol, 2002). However, the formula (9) in Hengl (2006) can be used to give an order of magnitude on the cartographic scale of the maps produced (here on the order of 1:26,000), providing users with caution for analysing spatial variability evidenced for shorter distances.

3. Applications of spatial interpolation for terroir analysis

Spatial interpolation of daily temperatures can be used to assess the potential impact of climate diversity on vine development and climate-related risks.

Figure 8 shows an example of its application in the study area. The number of frost days (T_min < 0 °C) from March to May in 2023 and 2024 is shown. Despite the relatively gentle topography, there are observed differences (measured by sensors) of 9 frost days in 2023: 3 days in the area at an elevation of 320 m a.s.l. north of the study area, compared to 12 days recorded by the sensor located in the plain, east of the study area, at an elevation of less than 230 m a.s.l. Spatial interpolation predicts an even wider variation of 3 to 18 days of frost in the same year. There is a certain recurrence in the spatial patterns of the number of frost days between the two vintages: frost occurs in enclosed areas (valley bottoms) in higher elevation areas, and especially at the bottom of hillsides, along transport infrastructure such as the railway line shown here. This information makes it possible to identify the areas most exposed to climate risks in order to apply appropriate management strategies (Snyder et al., 2005; Snyder & Melo-Abreu, 2005).

Finally, are two years sufficient to describe the spatial variability of temperatures at this scale? Previous studies have shown that, due to the strong dependence of the spatial structure of thermal fields on terrain, the average temperature distribution tends to repeat itself from one year to the next (Bois et al., 2018; de Rességuier et al., 2020). However, from one day to the next, the structure may differ, which justifies evaluating the links between terrain characteristics and measured temperature on a daily basis. In the present study, we observe a high degree of consistency in the spatial structure of minimum, mean, or maximum temperatures averaged over the month. This is particularly true for monthly minimum temperatures: those from 2023 are spatially correlated with those from 2024. This suggests that two years, or even a single year of study, is sufficient to characterise the average spatial structure of minimum temperatures for a given month (Figure S4). This is also the case for maximum temperatures, even though the values from certain sensors differ significantly from one year to the next (Figure S5). The correlation between the spatial structure of average temperatures is weaker, particularly in September (Figure S6).

Unfortunately, our dataset does not allow us to answer the following question: “How many years are needed to reliably describe the spatial distribution of temperatures and their potential impact on crops (in this case, grapevines)?” That is not the purpose of this study. However, the modest monthly comparison and the analysis of the scientific literature presented above seem to suggest that a few years (2 to 5, or even 1 year to study the impacts of minimum temperatures) should suffice.

4. Concluding remarks

The characterisation of wine-growing terroirs relies on rigorous description of the physical environment, particularly the fine-scale variability of the climatic conditions. This study focused on the question of the sensor density required to accurately map the thermal potential of vineyards through the interpolation of daily minimum (T_min) and maximum (T_max) at a local scale. While no universal rule dictates the relationship between station network design and the resulting accuracy of temperature spatial interpolation at a local scale, our findings suggest that the required sampling density might be determined by the intrinsic spatial scale variability of each thermal variable. The choice of the number of sensors, often constrained by economic considerations, must be evaluated on a case-by-case basis according to the targeted spatial accuracy, the topographical diversity of the study area, and the state of scientific knowledge in topoclimatology. However, a few useful points can be drawn from the results of this study and the literature cited in this article:

• The use of terrain descriptors calculated from high-resolution Digital Terrain Models (resolution < 5 m) makes it possible to consider phenomena such as the shading effect or the nocturnal cold air drainage that can lead to thermal differences with significant consequences for crops.
• The practical spatial resolution of a map can be assessed based on the size of the study area and the number of sensors deployed (Hengl, 2006).
• Across varying geographic scales, the predictive skill of daily temperature models remains relatively consistent, with RMSE values typically ranging from 0.1 to 1.5 °C. Our findings align with established literature at both the local (Cho et al., 2020; Le Roux et al., 2017) and mesoscale (Berndt & Haberlandt, 2018; Bois et al., 2018) levels, confirming the robustness of the geostatistical approach.
• In scenarios where few sensors are available, the use of Machine Learning frameworks, specifically Random Forest, could improve accuracy. By resolving non-linear interactions between geomorphometric variables and localised thermal signals, these algorithms maintain higher stability than traditional linear frameworks in data-sparse environments.

Ultimately, the challenge in implementing a high-resolution climate monitoring lies in a strategic trade-off: accurately sampling the topographical features that shape the landscape of the study area while accounting for the “thermal noise” of vineyard management and the influence of microclimatic factors. At a height of 60 cm, microclimatic factors such as localised soil management or canopy architecture introduce non-topographical signals that can significantly alter the representativeness of individual measurements. This study demonstrates that for terroir-level precision, the network must be sized not only to the area being studied but to the complexity of the thermal signal being tracked.

Acknowledgements

This study was funded by the Organisme de Défense et de Gestion (ODG) des Grands Crus de Chambolle-Musigny et Morey-Saint-Denis, the ODG of the AOC Morey-Saint-Denis, the ODG of the AOC Chambolle-Musigny et Vougeot, the Comité Bourgogne (BIVB) with the support of the French National Research Agency (ANR) as part of the France Relance program and the call for projects for the “Recovery Plan: Preservation of R&D Jobs” initiative funded by the European Union—NextGenerationEU. The authors would like to thank all the organisations that have financially supported the MusiMoClim project mentioned above, as well as all those who have contributed to the work carried out as part of this project. Thanks also to Myriam Kryshinets, Corinne Quintar, and Beatrice Casas for their administrative support, as well as to all the winegrowers who made their vineyards available for the installation of sensors and weather stations, for the observation of phenological stages, and for the collection of grape samples