Introduction

During the last decade, the cultivation of almond trees (Prunus dulcis) has experienced a major boom worldwide, particularly in Europe where it has been pushed by both higher demand and price (Lorite et al., 2020; Torres-Sánchez et al., 2018). Spain is the main producer of almonds in Europe and the third largest in the world after the United States and Australia, with a 31.1% increase in cultivated area between 2010 and 2020 (from 547 822 to 718 540 ha). Much of this increase was due to new intensive or super-intensive orchards under irrigation (in Spain from 40 855 ha in 2010 to 118 202 ha in 2020), which represented a 290% increase in such types of orchard (source: Spanish Ministry of Agriculture, Fisheries and Food). Such numbers show how almond has evolved from being an almost marginal crop to become a focal point of attention amongst woody crops. In this regard, an important aspect to highlight is the advances made in training and pruning systems that have allowed a process of continuous intensification and greater mechanization of tasks such as pruning and harvesting (Iglesias, 2020). In addition, the almond has been recognized as a fruit species with the application of increasingly specialized production technologies. Nevertheless, although mechanization for intensive crop management has increased (Bechar & Vigneault, 2016) most modern orchards are still uniformly managed and generally no consideration is given to possible intra-field spatial variability. Understanding this variability is key to achieving greater efficiency and sustainability in food production and forms part of the roadmap of the Common Agricultural Policy of the European Union within the 'Green Deal' and 'From Farm to Fork' strategies (European Commission, 2019).

The measurement of geometrical and structural tree parameters is an essential but labour-intensive activity in fruticulture. Canopy size, tree height, porosity and crop vigour are key factors to take into account for orchard management tasks such as training or pruning, which are major contributors to crop production costs (Díaz-Varela et al., 2015). To achieve a more sustainable and precise management of orchards, different precision agriculture technologies have started to be developed at research level in the last decade (Arnó et al., 2013; Colaço et al., 2018; Escolà et al., 2017). These include LiDAR (Light Detection And Ranging), a technology that is acquiring growing importance because it allows very precise geometric and structural parameters to be obtained, including height, width, volume, leaf density and porosity (Escolà et al., 2017). One example of its application can be found in the works of Gené-Mola et al., (2019, 2020) who used a mobile terrestrial laser scanner (MTLS) to detect fruits in an apple orchard. They achieved success rates higher than 80% due to the higher infrared light reflectance of apples in relation to leaves and trunks. Other works have addressed the characterization of tree volume and canopy density for precision spraying (Gu et al., 2021; Mahmud et al., 2021). LiDAR point clouds have also been applied to the precise mapping of tree branch structures and their automatic classification to support precision pruning (Zhang et al., 2020). In this line, a procedure for automated tree pruning using LiDAR scans was also proposed by Westling et al. (2021). However, LiDAR point clouds still constitute voluminous data sets, which are difficult to convert into useful information to support management decisions (Gregorio & Llorens, 2021). This constitutes a particular bottle neck for the application of this technology in orchard management, where monitorization of canopy geometry and structure at different times along the annual crop cycle is necessary for precise management.

Another alternative to map orchard geometric parameters is the use of structure-from-motion photogrammetry from unmanned aerial vehicles (UAVs). This provides a large capacity to monitor the development and dynamics of tree growth and structure over time (Johansen et al., 2018). Zhang et al. (2019) reviewed the use of UAVs in orchard management, describing sensing devices, platforms and data processing approaches for the measurement of biophysical and geometric parameters, the detection of diseases and health status, or yield estimation, among others. Their main conclusion was that, although UAVs have numerous potential uses in orchard management, they are still largely underutilized for several reasons. For example, the vast amounts of data require significant processing and specialized operators are needed to derive, clean and obtain useful information from the point clouds. It should also be noted that UAV operating restrictions may be applicable and licensed pilots required in some countries and/or production regions. In addition, there are commonly differences in the quality of the UAV images associated with the speed, height or direction of the flight, and there may be additional problems related to the presence of shadows depending on the azimuth and inclination of the sun at the time of acquisition (Tu et al., 2020). This has implications in terms of the alteration of the radiometric values of the sensed covers. Therefore, UAV imagery still presents low field efficiency for large farm operations (Sozzi et al., 2021). However, although the use of UAVs in the monitoring of commercial orchards has not been extensive, some important research studies in this regard have been conducted. These include the quantification of pruning impacts on olive tree architecture and annual canopy growth using UAV-based 3D modelling (Jiménez-Brenes et al., 2017), the mapping of the 3D structure of olive trees using point clouds and object-based image analysis (Torres-Sánchez et al., 2018), and a study of the phenotypic variability of different almond varieties by mapping tree height, volume, flowering dynamics and flower density from photogrammetric 3D point clouds and object-oriented based images analysis (OBIA) (López-Granados et al., 2019).

A third alternative for data acquisition and orchard monitoring is satellite imagery, which offers an increasingly high temporal and spatial resolution. More specifically, the Sentinel-2 mission (Copernicus Program, European Union) provides a spatial resolution of 10 m and a time revisit of 5 days considering the two satellites of the mission (Sentinel-2A and Sentinel-2B) (European Space Agency (ESA), 2015). Other missions or satellite constellations include PlanetScope (Planet Labs Inc.), which offers daily images with very high resolutions of up to 3 m (Planet Labs Inc, 2021). However, when considering crops with discontinuous patterns of vegetation, such as hedgerow orchards, the inter-row presence of weed vegetation and/or bare soil can affect the pixel values when computing spectral indices. To solve this drawback, some works have focused on extracting pure canopy pixels from UAV images using crop-specific masks (Khaliq et al., 2019; Pastonchi et al., 2020). However, although the best correlations with satellite vegetation indices (VIs) were obtained for NDVI (normalized difference vegetation index) maps computed with mixed UAV pixels, as reported by Pastonchi et al. (2020), neither satellite nor UAV images were able to correctly detect the spatial variability of vigour within the vineyard.

Despite the above-mentioned drawbacks, the use of satellites to monitor crops and support management actions is increasingly attracting the attention of scientists, commercial service companies and farmers. From the scientific point of view, one of the latest developments can be found in the work of Abdelmoula et al. (2021) who established a method to assess the spatio-temporal dynamics of biophysical parameters in olive orchards using Sentinel-2 images fused with PlanetScope and radiative transfer models. This allowed parameters such as leaf area index (LAI), chlorophyll and water content, and mesophyll structure to be retrieved. Barajas et al. (2020) evaluated the organoleptic characteristics and the agronomic, phenological and nutritional quality of pistachios (Pistacia vera L.) based on the NDVI, calculated in the nut-filling phenological stage from Sentinel-2 images. No significant differences were observed between vigour levels, although the most vigorous trees had nuts with a higher percentage of fibre and protein. In a study on almonds, Chen et al. (2019) developed an enhanced bloom index to quantify flowering status in orchards of the Central Valley of California from series of PlanetScope and Sentinel-2 imagery. There has also been growing interest in recent years in mobile and web applications that offer vegetation vigour monitoring based on Sentinel-2 and/or PlanetScope images. Such applications, examples of which include Auravant (Buenos Aires, Argentina), E-STRATOS (Barcelona, Spain), and SatAgro (SatAgro Sp., Warsaw, Poland), usually integrate different data sources to provide a detailed description of crop status which allows the user to make appropriate management decisions. In some cases, the application is customized according to the needs of the individual client, making interpretation of the information easier as a support for the decision-making process.

Nonetheless, despite the advantages offered by satellite imagery and the increasing existence of services based on continuous monitoring capabilities, there remains a gap that needs to be bridged in terms of the estimation of geometrical and structural vegetation parameters of orchards on the basis of the applications and services available to farmers. More specifically, in-depth knowledge is required about the relationships between geometric and structural parameters of the fruit tree architecture and canopy density, which are used in the orchard management decision-making process, and the vegetation indices that can be derived from satellite images.

In this line, the present work aimed to study the relationships between ground-derived geometric and structural parameters of super-intensive almond orchards and vegetation indices derived from satellite images. The goal was to use the information obtained to support management actions, such as differential pruning and the application of plant protection products, and to facilitate the delimitation of within-field management zones on the basis of VIs, potentially reducing the dependence on more time-consuming processing techniques such as 3D LiDAR scans or UAV photogrammetry. Moreover, the results could serve as a basis for the use of this type of relationship in commercial applications based on the monitoring of fruit orchards from detailed satellite images.

Accordingly, the objectives of the present work were: (a) to explore the spatial relationships between LiDAR-derived canopy parameters and vegetation indices derived from multispectral satellite images in a super-intensive almond orchard; (b) to define potential management zones at different times of the crop cycle along with a campaign based on satellite images that reliably reproduce the within-field variability in canopy characteristics; and (c) to compare the cluster maps based on vegetation indices with the cluster maps derived from the interpolation of geometric and structural canopy parameters in order to evaluate their consistency.

Materials and methods

Study area

The study area comprised a 0.75 ha block of a super-intensive almond orchard (Prunus dulcis cv Lauranne Avijor) located in Raimat (Catalonia, NE Spain, X 288 260, Y 4 616 100, ETRS89 UTM 31 Zone T, Fig. 1), with a plantation pattern of 3.2 × 1.5 m (2800 trees ha−1). The tree rows form a continuous and vertical wall (Fig. 1). The orchard was planted during the 2016/17 winter. The orchard is pruned each year in winter to create a hedge of low vigour branches (fine, short with a certain horizontal tendency). Then, at the end of the spring, the orchard is mechanically pruned to maintain an efficient and active exposed leaf area to develop the maximum productive potential, and to facilitate the work of harvesting machines. The width of the vegetation hedge is usually adjusted to a distance of about 0.35 m on each side of the central axis of the three rows. With that, lighting and aeration conditions are improved. In the present case study, the spring mechanical pruning was done on June 19th, 2019.

Fig. 1
figure 1

Location of the case study area. Right: almond field and block where LiDAR data was acquired. Bottom left: view of the super-intensive almond orchard

Full size image

Methodological scheme

LiDAR data, PlanetScope and Sentinel-2 satellite images from two dates (June and September 2019) were acquired and analysed. The data analysis procedure is shown in Fig. 2 and the acquisition and processing methods are explained in detail in the following sections.

Fig. 2
figure 2

Overview of the method workflow

Full size image

Canopy orchard parameters were obtained from the LiDAR point clouds. The LiDAR datasets were georeferenced and summarized every 0.5 m along the rows. After cleaning the data for possible outliers, a spatial analysis of LiDAR-derived parameters was carried out to know their spatial autocorrelation. These parameters were later mapped using ordinary kriging to be compared with VIs from satellite images.

A canonical correlation analysis (CCA) was then carried out between two sets of variables (a LiDAR-derived geometric and structural parameter set and a satellite multispectral VI set). The objective was to identify the main inter-relationships between geometric and/or structural LiDAR-derived canopy parameters and satellite VIs.

Once the most inter-related parameters and VIs had been identified, a cluster analysis was carried out to delineate potential management zones. For that, a series of analysis of variance (ANOVA) tests was carried out to discriminate significant differences in the canopy parameters from the vegetation indices cluster zones (k = 2). Finally, to complement this cluster analysis, a map class comparison of the reclassified cluster VIs and LiDAR maps was performed to evaluate their consistency by means of the fraction correct (FC) and Kappa index (KI) of agreement.

LiDAR data acquisition and processing

LiDAR data was acquired on June 22nd and September 20th, 2019, with a self-developed mobile terrestrial laser scanner using a VLP-16 LiDAR sensor (Velodyne Lidar Inc., Silicon Valley, USA). This sensor has a range of 100 m and dual return capability. It uses 16 simultaneous laser beams with a 30° field of view to acquire ~ 300 000 points s−1 with a 360° scanning window. These data were georeferenced with a Leica GPS 1200 GNSS-RTK system (Leica, Wetzlar, Germany). The LiDAR sensor was mounted in a self-propelled mobile platform moving at a constant speed of 2 km h−1.

Data from 24 tree rows (84 m long each) were acquired during the survey, with a total number of 3D points of 469 million. This point cloud was processed with RStudio software (Version 1.4.1717 using R version 4.0.3 as calculation engine), with an own-developed R code described in Llorens et al. (2019). By applying this code, different geometric and structural parameters of the canopy were summarized every 0.5 m along the tree rows (Table 1). In the processing algorithm, two regions of interest (ROIs) are considered (Fig. 3): ROI-A are rectangular parallelepipeds along the almond rows with centroids each 0.5 along the rows and ROI-B are parallelepipeds of 0.10 m in height along the vertical axes passing by the ROI-A centroids. ROIs A and B were used consecutively to extract the LiDAR crop parameters from all the rows of the field. To reduce the point cloud to be analysed, only the central LiDAR beam was used to extract information from the point cloud.

Table 1 Description of the geometric and structural vegetation parameters extracted from the LiDAR point cloud along the tree rows
Full size table
Fig. 3
figure 3

Graphical description of the different regions of interest A, B used in the process of crop parameters extraction (adapted from Llorens et al., 2019). ROI-A is shown as a zenithal view and ROI-B as a cross-sectional view

Full size image

Satellite data acquisition and processing

To compute the different VIs used in this study, six satellite images were acquired (Table 2).

Table 2 Satellite images dataset used in the present study
Full size table

All images were of the bottom-of-atmosphere (BoA) surface reflectance type and projected onto the WGS84 UTM 31 N co-ordinate system, as provided by the Copernicus and Planet Labs services. The PlanetScope satellite constellation consists of multiple individual nano-satellites called Doves (Planet Labs Inc, 2021). To ensure homogeneity with respect to the calibration of the different sensors of the PlanetScope satellites, the VIs calculated from the June 25th and 27th images were averaged. The same procedure was applied to calculate the VIs from the September 19th and 20th images.

The first step for image analysis was the calculation of different VIs. Four of the most quoted VIs in the scientific literature were calculated (Xue & Su, 2017): the enhanced vegetation index (EVI) (Huete et al., 1997), the NDVI (Rouse et al., 1973), the green normalized difference vegetation index (GNDVI) (Gitelson et al., 1996), and the green–red vegetation index (GRVI), which is a good index to detect changes in crop phenology (leaf-colour change) (Motohka et al., 2010).

Mapping of LiDAR-derived canopy parameters on satellite pixels

As mentioned above, LiDAR data were summarized in points along the tree rows every 0.5 m (ROI-A) (Fig. 3). Then, for the purpose of comparison with the VIs from the Sentinel-2 and PlanetScope images, geostatistical interpolations of each LiDAR-derived parameter were carried out. These interpolations were performed using QGIS and its Precision Agriculture Tools (PAT) plugin (Ratcliff et al., 2019), which uses VESPER software (Variogram Estimation and Spatial Prediction Plus Error) (Minasny et al., 2005). The central co-ordinates of the image pixels covering the study area were used as the grid points for variable interpolation. Only pixels that fully included ROI-A centre points were considered (Fig. 4).

Fig. 4
figure 4

LiDAR data along the tree rows with an indication of the centroids of the interpolation pixels (PlanetScope and Sentinel-2 grids) where the LiDAR-derived geometrical and structural orchard parameters were interpolated. The background image is UAV-sourced and was taken on June 26th, 2019

Full size image

For geostatistical interpolation, Webster and Oliver (2007) suggest that, when more than 100 sample points are considered, the best method to apply is ordinary local block kriging. This type of kriging estimates the values of the variable at the grid points taking into account the value of the original sampled points included in a mobile window defined by the user. Ratcliff et al. (2019) suggested that it should be approximately 5 times the pixel size which, in this case, is 15 m for PlanetScope images and 50 m for Sentinel-2 images. Before kriging, and following the criteria of Taylor et al. (2007), values above or below ± 2.5 standard deviations (SD) were considered outliers and removed from the LiDAR ROI-A data set. Therefore, in this study, a total of 1595 ROI-A points were considered as input for local block kriging using the PAT QGIS plugin (high-density kriging).

Canopy orchard parameters can be distributed unevenly at reduced distances. In such cases, exponential or spherical models are the most suitable to fit semivariograms for kriging interpolation (Millán et al., 2019). Therefore, in this study exponential models were used. As mentioned above, two interpolation grids (block grids in the PAT QGIS plugin) were defined: one with 697 central points for 3 m resolution of PlanetScope images and one with 54 central points for 10 m resolution of Sentinel-2 images.

Analysis of the spatial structure of LiDAR-derived canopy parameters

Changes in local variability are reflected in the semivariogram parameters (Oliver, 2010). The spatial dependency of each LiDAR-derived parameter was tested using the Cambardella Index (CI) (Cambardella et al., 1994) (1) and range (A1) which is the distance at which some variograms reach the sill. Point samples that are at a higher distance than the range are supposed to be spatially independent (Oliver, 2010).

$$CI = { }\frac{C0}{{C0 + C1}}{*}100$$
(1)

where C0 is the nugget variance, C1 is the estimate of the spatial structural variance, allowing C0 + C1 (sill) to be obtained as a measure of the total variance contained in the data.

According to Cambardella et al. (1994), if the CI is less than or equal to 25%, the variable has a strong spatial structure (or spatial autocorrelation), if the ratio is between 25 and 75%, it has a moderate spatial dependency, and if it is greater than 75% the variable has a weak spatial dependency. In other words, the higher the CI, the worse the spatial structure of distribution of a parameter. When C0 is high, parameters change a lot in short distances and, therefore, the semivariogram cannot explain correctly the total variance along the orchard. Therefore, a higher probability exists of obtaining a randomized interpolated map.

Another important parameter to determine how those properties vary spatially is the range (A1), which evaluates the adequacy of the interpolated maps in terms of reflecting the spatial distribution of variability. The higher the range, the more homogeneous distribution a property can have.

Since ordinary kriging was made based on local semivariograms, the Cambardella Indices (1) and range (A1) were calculated from the average value of the semivariogram parameters, referred to as the ROI-A centroids (Table 4).

Correlation between geometrical and structural orchard parameters and vegetation indices

To study the relationship between the two considered groups of variables (LiDAR-derived parameters and VIs), a canonical correlation analysis was carried out. The aim was to select the VIs that best explain the variability of the geometry and structure of the vegetation for the definition of potential management zones.

A CCA is a generalization of multiple regression used to highlight associations between two datasets U and V containing a defined number of variables Y1, Y2, …, Yn and X1, X2, …, Xn, respectively. In this method, several Y variables are simultaneously related to several X variables. In the present case, set U was composed of the four multispectral indices (Y1-EVI, Y2-NDVI, Y3-GNDVI, Y4-GRVI) and set V included the four LIDAR-derived parameters (X1-MaxHeight, X2-MaxWidth, X3-SMaxWidth, X4-Porosity-Avg). After the estimation of linear combinations among the parameters corresponding to the same set, the method searches first for a pair of canonical variables (U1 and V1) that, each being a linear combination of the corresponding set variables, show the highest possible correlation between them. The procedure continues looking for new canonical variables (U2 and V2), and so on, until all the variability is explained (Fig. 5). This procedure allows hidden relationships between variables to be discovered (Manly & Navarro, 2017). This analysis was implemented using RStudio (Version 1.2.5001 using R version 3.6.1) and the CCA R package, freely available from the Comprehensive R Archive Network (C RAN, https://cran.r-project.org/).

Fig. 5
figure 5

Canonical statistical procedure. Set U includes the four multispectral indices (Y1-EVI, Y2-NDVI, Y3-GNDVI, Y4-GRVI) and set V the four LIDAR-derived parameters (X1-MaxHeight, X2-MaxWidth, X3-SMaxWidth, X4-Porosity-Avg)

Full size image

Definition and analysis of management zones

The variables selected from the CCA were clustered using a k-means algorithm (Hamerly & Elkan, 2003) to delineate potential management zones. The k-means algorithm implemented in Google Earth Engine was used. Two clusters (classes) were created. To determine whether the clusters were able to effectively differentiate almond tree canopies with different characteristics, a series of one-factor ANOVA tests was carried out, with the LiDAR-derived parameters as the dependent variable and the clusters the analysis factor of the model. Three t-student tests were carried out to assess statistical differences among the groups. To minimize errors of Type I the Bonferroni adjustment was applied (Manly & Navarro, 2017).

Cluster maps comparison

To assess whether LiDAR-derived parameter clusters could be estimated from the cluster maps defined from the multispectral indices, a comparison of cluster maps was carried out.

For this purpose, the Map Comparison Kit (MCK version 3.2.3 http://mck.riks.nl/) software was used. This program provides a good understanding of the differences between pairs of maps (Visser & De Nijs, 2006). In this case, the map pairs were compared using two statistics: the FC and the KI. The FC (fraction correct), or percentage of agreement, is calculated as the number of equal cells in both maps divided by the total number of cells in the map. This statistic is considered flawed as an overall measure of similarity because it tends to consider maps with few or unevenly distributed categories to be more similar than those with many, more equally distributed categories (Visser & De Nijs, 2006). For a better-balanced measure of similarity, the Kappa index of agreement can be used. It is based on a contingency table that summarizes the cross-distribution of categories over the two maps. This index is the result of multiplying two statistics (2): Kappa Histo (KHisto), a measure of the agreement of the total number of pixels considered for each class category among the two maps, and Kappa Location (KLoc), which is sensitive to differences in the spatial distribution of classified pixels in the two maps.

$$Kappa = K_{Histo} \times K_{Loc}$$
(2)

Results

Spatial structure of the orchard canopy geometric and structural parameters

Table 3 shows the main descriptive statistics for the original LiDAR data aggregated at 0.5 m along the orchard rows. According to these results, the geometric parameters (MaxHeight, MaxWidth and SMaxWidth) presented considerably lower coefficients of variation (CV) than the structural one (Porosity-Avg) in both June and September. Table 4 shows the Cambardella parameters (C0 and C1), and the Cambardella indices (CI) and ranges (A1) calculated from the parameter average of the local semivariogram after the local kriging interpolations to PlanetScope and Sentinel-2 grids. Figure 6 offers an example of the LiDAR-derived parameters maps corresponding to June interpolated to the PlanetScope and Sentinel-2 grids.

Table 3 Descriptive statistics of the original normalized LiDAR-derived geometrical and structural orchard parameters
Full size table
Table 4 Averaged local semivariogram parameters from interpolated ROI-A data to PlanetScope and Sentinel-2 grids
Full size table
Fig. 6
figure 6

Interpolated LiDAR-derived parameters to the satellite grids: PlanetScope (top), Sentinel-2 (bottom) in June 2019

Full size image

Regarding the CI values (Table 4), the geometrical parameters reflected a better spatial structure than the porosity parameter. In addition, spatial structure for all parameters was higher after spring pruning (June) than before harvesting (September). Porosity presented the poorest spatial autocorrelation (CI values around 70%) due to high values of C0 and C1 (low explained variance). This property presented higher range values, meaning that this poor spatial autocorrelation is continuous throughout the orchard. Maximum height (MaxHeight) was the geometrical parameter with the best spatial autocorrelation after mechanical pruning (CI values ≤ 25%). This parameter also presented the most important changes from June to September (from 20–24% to 50–56%). The width-related parameters (maximum width and cross-section) showed a moderate-to-high spatial autocorrelation (CI between 33 and 40%).

CI values were similar for both PlanetScope and Sentinel-2 semivariograms. Therefore, as reflected in Fig. 6, the spatial structure of distribution of the different parameters was very similar in both sets of interpolated maps. There existed a high correspondence between the width-related variables and the inverse correspondence of them with porosity. In addition, maximum height presented a different gradient of distribution along the rows.

Relationships between orchard canopy geometrical and structural parameters and vegetation indices

To study the associations between the groups of variables considered, a CCA was performed considering Set U (MaxHeight-MH; MaxWidth-MW; SMaxWidth-SMW; Porosity-Avg-P) and Set V (EVI-(E); NDVI-(N); GNDVI-(GN); GRVI-(GR) (Fig. 5).

The first pair of canonical variables CanR1 (Xcan1 and Ycan1) showed correlation coefficients between 0.71 (PlanetScope image) and 0.81 (Sentinel-2 image) (Table 5). The correlation coefficients of the second canonical variable (CanR2) were lower but still significant, except for the Sep S-2 dataset. Figure 7 shows the two-dimensional graphs that represent the weight of the original variables in the highest significant canonical variables CanR1 (Xcan1, Ycan1) and CanR2 (Xcan2, Ycan2). The first dimension reflected direct correlations of NDVI and GNDVI with the geometrical parameters (highlighting MaxWidth and SMaxWidth), as well as inverse relationships of those indices with Porosity-Avg. This trend was identified in both dates and both satellite datasets. The second dimension was more significant and stronger in June for both the PlanetScope and Sentinel-2 datasets. MaxHeight was the geometrical parameter with the worst relationship in general with all the VIs.

Table 5 Canonical correlation analysis significance
Full size table
Fig. 7
figure 7

Canonical correlation analysis graphs. CanR1 means the correlation (Xcan1, Ycan1). CanR2 means the correlation (Xcan2, Ycan2). Xcan are the canonical variables of the geometrical and structural orchard parameters and Ycan the canonical variables of the vegetation variables. The triangles (MH (MaxHeight), MW (MaxWidth), SMW (SMaxWidth) and P (Porosity-Avg)) represent the (Xcan1, Xcan2) and the points (EVI (E), NDVI (N), GNDVI (GN), and SAVI (S) the (Ycan1, Ycan2)

Full size image

Analysis of potential management zones

Regarding the CCA analysis, the first canonical variable reflected direct and strong correlations of geometrical parameters with both NDVI and GNDVI. These VIs also reflected an inverse, but also significant, relationship with Porosity-Avg. Therefore, the NDVI and GNDVI from PlanetScope and Sentinel-2 images were clustered using the k-means algorithm and a cluster analysis performed. In addition, only the canopy orchard parameters that presented the best performance in the CCA (MaxWidth, SMaxWidth and Porosity-Avg) were considered for the ANOVA tests.

Figure 8 shows a basic statistics summary corresponding to the NDVI and GNDVI vigour classes (clusters) from the PlanetScope and Sentinel-2 images of June and September 2019.

Fig. 8
figure 8

NDVI and GNDVI from PlanetScope and Sentinel images in June (upper part of the figure) and September (bottom part of the figure) in the experimental almond orchard. Cluster a = High values; Cluster b = Low values. N is the number of pixels in the class

Full size image

From Fig. 8, it can be observed that the NDVI and GNDVI cluster maps showed similar spatial patterns and, in both cases, small differences were found in the mean and standard deviation values. However, LiDAR-derived parameters were significantly differentiated (p-value < 0.001), as reflected in Table 6. This table presents the ANOVA t-student tests (three for NDVI and three for GNDVI) carried out to assess the variability of the LiDAR-derived parameters in the delimited potential management zones According to the Bonferroni correction, given the total number of ANOVA tests (three) with an overall significance level of 0.05, the probability of declaring each univariate ANOVA test as significant was 0.016.

Table 6 Analysis of variance (ANOVA) of LiDAR-derived parameters according to NDVI and GNDVI (2 classes)
Full size table

As Table 6 suggests, both NDVI and GNDVI were useful to separate different classes of the LiDAR-derived parameters. However, NDVI seems to have a better performance than GNDVI in general. For that reason, for the rest of the cluster analysis, only the NDVI was utilised. In order to check the correspondence between the NDVI and the LiDAR-derived cluster classification (k = 2), a map comparison was carried out. The results of this analysis are presented in Table 7 and Fig. 9.

Table 7 Statistics calculated for the 2 cluster maps from NDVI and LiDAR selected parameters in June and September 2019: FC (Fraction correct), KI (Kappa index), KLoc, KHisto
Full size table
Fig. 9
figure 9

Within-field maps of the NDVI and LiDAR-derived parameters

Full size image

As can be seen in Table 7, the FC was higher than the KI. The reason was explained in the methodology section. FC is a method which only considers the percentage of change of pixels, with the KI being more restrictive. Despite this, both indices showed the same trends. In accordance with Landis and Koch (1977) criteria, width-related parameters (MaxWidth and SMaxWidth) showed moderate KIs (0.41–0.60) and, slight to fair agreements were found between the VI cluster maps and Porosity-Avg (around 0.30). These trends can be appreciated in Fig. 9, where the correspondence of the spatial patterns of the NDVI and the LiDAR-derived parameters interpolated to the satellite image grids is shown. As for the interpolated maps, despite their differences in spatial resolution, the cluster distributions of PlanetScope and Sentinel-2 were very similar. However, higher correspondences were found between the LiDAR clustered maps and the NDVI from Sentinel-2.

Regarding the KHisto and KLoc indices, high agreements were found. KHisto was high in general, which means that the total number of pixels of each class did not experience significant changes between pairs of maps: MaxWidth (0.71–0.78), SMaxWidth (0.78–0.97), porosity-Avg (0.85–0.97). With respect to the KLoc values, which refers to the spatial location agreement, moderate correspondences were obtained when width-related parameters (MaxWidth and SMaxWidht) and NDVI within-field zones were compared (around 0.65). Slight to fair agreements were reflected between both maps in relation to porosity (around − 0.37).

Discussion

The aim of the present work has been to determine the suitability and accuracy of PlanetScope and/or Sentinel-2 satellite images to estimate canopy-related geometric and structural parameters derived from terrestrial LiDAR surveys in hedgerow almond orchards and to use this information to define potential within-field management zones. Such research could represent an advance in the continuous monitoring of super-intensive orchards, facilitating the interpretation of VIs in terms of different key geometric and structural parameters and hence the work of farmers and application services in fruticulture.

One of the challenges in this research was related to the differences in data resolution between the LiDAR point clouds (0.5 m) and the VIs computed from satellite images (3 or 10 m pixel resolution). However, results are encouraging, with significant differences found in the relationships between some of the variables. From the CCA analysis, higher NDVI and GNDVI values were strongly related to higher values of canopy width and cross-sectional area (MaxWidth, SMaxWidth, respectively). In contrast, higher values of the same VIs corresponded with lower values of canopy porosity (Porosity-Avg). This trend was identified in both satellite datasets. These findings are relevant because the geometric parameters (particularly the width-related ones) are important to determine canopy volume and because porosity is related to leaf density and, hence, to the optimization of dose rates of plant protection products (Arnó et al., 2013; Escolà et al., 2013). In this respect, and in relation to canopy volume and in particular in almond orchards, Underwood et al. (2016) presented a mobile terrestrial scanning system to predict yield for individual trees. Their approach obtained an R2 = 0.77 between canopy foliage volume and almond yield. In addition, Lampinen et al. (2011) and Maldera et al. (2021) established that harvesting efficiencies depended on tree canopy size and that yield production was strongly related to the previous year’s spur leaf area. However, canopy width is not the only parameter with a crucial impact on orchard management. Among other factors, porosity is considered a very important element to decide how to handle certain orchard management tasks (Castillo-Ruiz et al., 2016; Duga et al., 2015). Porosity constitutes a good indicator for total light interception and its distribution within the canopy. This is very important for spur health and is essential to estimate the application of plant protection products. Additionally, the balance between yield production and tree vegetation development is a key factor when planning an orchard management campaign. With the above in mind, being able to achieve a regular control (inter-season and intra-season) from VI-based estimation of width parameters, such as maximum width and cross-sectional area, and porosity distribution in the orchard can constitute an advantage for the adoption of appropriate management practices that can deliver higher spur productivity and yield.

However, not all the parameters have the same opportunity to be homogeneously spatially distributed in management zones. A key factor presented in this work was the study of the spatial distribution of the different LiDAR parameters within the orchard at two important times in the campaign. The spatial autocorrelation analysis of the variables through the Cambardella index showed the existence of different spatial structures between the geometric parameters and porosity. The local semivariograms related to canopy geometry were suitable to detect how these parameters varied with distance (low C0), showing that the spatial continuity of these variables within the orchard had a moderate-to-high spatial autocorrelation (CI between 33 and 40%). In contrast, porosity showed high variability at very short distances (high C0). This property had considerably higher coefficients of variation than the geometric parameters, and was the variable with the lowest spatial dependency (highest CI). As suggested by Manly and Navarro (2017), the high spatial variability of porosity could potentially degenerate into very fragmented and patchy spatially distributed management zones if based on this property. However, the porosity range showed that the spatial continuity of this high variability remained at a higher distance than in the case of the geometrical parameters and, as can be appreciated in both Figs. 6 and 9, the value of pixels was not overly haphazardly distributed. This could be due to use of the block kriging method in the interpolation, since it manages to smooth the pixel values.

From the ANOVA, it can be observed that both the NDVI and the GNDVI were able to detect small but significant differences in the geometrical and structural orchard parameters. However, in relation to both the ANOVA test and map comparison analysis, the NDVI can be considered more adequate for creation of potential management zones. Other authors have highlighted the suitability of using the NDVI for canopy characterisation in fruit orchards because of its relationship with canopy-related geometric and structural characteristics (Caruso et al., 2017; Johnson et al., 2003). In other horticultural crops, such as grapevine, Martínez-Casasnovas et al. (2012), in their analysis of differential management zones based on NDVI classes concerning vine development, grape maturity and quality, found that differentiation of vine and grape juice characteristics in two zones was optimal for the study target.

Evidence of favourable interpolation results was found when assessing the consistency between the potential management zones delimited from the NDVI and LiDAR parameters maps by means of analysis of similarity (Table 7, Fig. 9). The percentage of agreement (Fraction Correct) indicated a moderate-to-high global concordance among the comparisons with LiDAR interpolated maps, highlighting the parameters related to width (MaxWidth and SMaxWidth). Regarding the Kappa index and KLoc and KHisto statistics, good global concordances were found among the NDVI cluster maps from PlanetScope and Sentinel-2 datasets and the cluster maps obtained from the LiDAR interpolated maps. This is due to the effect of smoothing related to the interpolation and is very important for simplifying orchard management. Thus, Sentinel is also appropriate for the application of different orchard management actions.

In the present research, a point to highlight is that, in spite of the differences in spatial resolution, a good correspondence between spatial distribution patterns from the PlanetScope and Sentinel-2 datasets was found. Sentinel-2 proved to be highly suitable to reflect the spatial variability of orchard geometrical parameters, as well as canopy porosity. This concurs with the findings of Pastonchi et al. (2020) who developed a multi-temporal comparison between satellite and ground data with UAV information in an experimental vineyard. They concluded that Sentinel-2 images allowed within-vineyard vigour spatial variability to be detected correctly. Results of the present study indicate that management zones delineated from the VIs could be used to define crop management practices associated with canopy width-related parameters and porosity. In consequence, the results of the present work can add value to current knowledge.

As described in the introduction section, other authors have demonstrated the utility of remote sensing images, in particular from UAVs, to derive accurate geometric canopy parameters in orchards (Torres-Sánchez et al., 2018), or with the combination of geometric parameters and VIs (Jurado et al., 2020). These are mainly based on measurements using photogrammetric techniques and, in some cases, with LiDAR data as reference (Hobart et al., 2020). However, monitoring based on UAVs or LiDAR 3D point clouds instead of satellites has the limitation of the frequency of image acquisition and time-consuming data processing, respectively.

Conclusions

The methods evaluated in this work can complement the possibilities of canopy characterization in super-intensive orchards based on high temporal and spatial resolution satellite imagery, as they are able to detect geometrical and structural changes in almond orchard canopies. In this respect, vegetation parameters related to width, cross-sectional area and porosity of the canopy along the rows offered a high correlation, especially with NDVI but also with GNDVI from Sentinel-2 and PlanetScope platforms. Few differences were found when clustering parameters, with VIs seeming to work for the monitoring of crop status and evaluation of the results of certain agronomic actions. This methodology could be interesting as an input to building a model approach, based on the age and physiological characteristics of trees and also considering the water, energy and nutrient balance, in order to simulate crop growth and better estimate yield production. Additionally, it can be a useful tool to improve field management software applications based on satellite imagery by adding geometrical and structural parameter estimations and hence contribute to the management decision-making process. Although this study was based on data from a super-intensive almond orchard, the methodology employed can be applied to other crops with hedgerow cropping patterns and the mapping of canopy parameters can be extended to larger orchards.