Using Digital Photography to Track Understory Phenology in Mediterranean Cork Oak Woodlands

Abstract

Monitoring vegetation is extremely relevant in the context of climate change, and digital repeat photography is a method that has gained momentum due to a low cost–benefit ratio. This work aims to demonstrate the possibility of using digital cameras instead of field spectroradiometers (FS) to track understory vegetation phenology in Mediterranean cork oak woodlands. A commercial camera was used to take monthly photographs that were processed with the Phenopix package to extract green chromatic coordinates (GCC). GCC showed good agreement with the normalized difference vegetation index (NDVI) and normalized difference water index (NDWI) obtained with FS data. The herbaceous layer displayed a very good fit between GCC and NDVI (coefficient of determination, represented by r² = 0.89). On the contrary, the GCC of shrubs (Cistus salviifolius and Ulex airensis) showed a better fit with NDWI (r² = 0.78 and 0.55, respectively) than with NDVI (r² = 0.60 and 0.30). Models show that grouping shrub species together improves the predictive results obtained with ulex but not with cistus. Concerning the relationship with climatic factors, all vegetation types showed a response to rainfall and temperature. Grasses and cistus showed similar responses to meteorological drivers, particularly mean maximum temperature (r = −0.66 and −0.63, respectively). The use of digital repeat photography to track vegetation phenology was found to be very suitable for understory vegetation with the exception of one shrub species. Thus, this method proves to have the potential to monitor a wide spectrum of understory vegetation at a much lower cost than FS.

1. Introduction

Tracking seasonal patterns in vegetation dates back centuries [1] and was traditionally conducted through human observation of the most common metrics of plant structural change, for example, leaf unfolding—the time of year when leaves begin to sprout—and senescence—the time when they dry at the end of the growing season [2]. The timing of these events has tremendous implications on ecosystems, since plants are the base of the food chain. For example, birds may become desynchronized from arthropods, their food supply, if leafing dates change [3]. In addition, understanding phenology is also relevant in the context of process-based forest models, since it determines how ecosystems respond to different climate scenarios that are used to predict the impact of climate change [4]. While human observation is still the most accurate technique to monitor phenology, it is a time-consuming process, generally limited to a small area and number of sample individuals.

More recent remote sensing methods that rely on optical sensors aboard satellites orbiting the earth allow us to grasp changes in vegetation structure at a global and regional scale and have largely taken hold due to their convenience [5]. Despite their practical advantages, sensors are at a great distance from the target of interest, producing several potential impediments to viewing the object of study that should be considered, such as illumination and observation geometry and atmospheric or even ground conditions [6].

The spatial resolution of satellites is also much coarser than traditional methods, creating a gap between space-born remote sensing and actual field observations, which brings limitations to phenological research at the ecosystem level. To bridge the gap and achieve a broader and more detailed perspective of ecosystem phenology, it is possible to integrate automated up-close measurements of selected points of the study area [7,8]. Some study cases install spectral sensors on eddy covariance towers [5,9], but this solution can be costly and simply represents too much investment up-front. As an alternative to these automated methods, other studies propose the use of digital photographs and the green chromatic coordinates index (GCC) to monitor plant phenology [5,10]. This index quantifies the amount and weight of green pixels in a photograph.

The use of digital photography can offer a higher spatial resolution at relatively low cost and effort [2]. It is important to note, however, that temporal resolution will depend on the amount of automation achieved. Furthermore, it has been shown that phenological transitions can be well characterized by digital repeat photography [11] despite making use of simple commercial cameras, but also that GCC is able to monitor vegetation dynamics at lower temporal resolutions [12].There are even studies that evaluate camera choice and image file format to verify if there is significant impact on the results, finding that, in general, the method is “sufficient to characterize canopy development for phenological research” [13] and that GCC is the most effective index to minimize the effects of different scene illumination. Furthermore, some research focused on comparing time series of picture-based indices with spectral vegetation indices from near-surface data, finding good correlations [5,10,14].

The major advantage of using digital photography over more sophisticated spectral sensors is that it requires little setup and investment while still allowing us to collect direct observations of diverse vegetation that can be linked to ecosystem functioning [15]. In the context of a fast-paced changing climate, this feature is an essential one.

In recent years, there has been an increase in temperatures and in extreme weather and climate events [16] along with a decrease in precipitation, which affects the regular activity of many life forms and is leading to changes in phenological patterns [17]. Some regions are particularly vulnerable, such as the Mediterranean, where precipitation is expected to decrease during the summer [18,19,20,21]. An increased aridity could have long-lasting impacts on groundwater availability, since there is reduced recharge, and it may lead to inevitable consequences on tree mortality and understory structure [22]. For this reason, a practical approach to phenological research is of great importance, as it allows us to easily obtain an indication of changing climate.

Mediterranean cork oak woodlands are extremely biodiverse agrosilvopastoral ecosystems that provide valuable environmental services, such as carbon storage and fire prevention [23], which are rooted in the structure of vegetation [24]. This is a structure characterized by a heterogenous understory and cork oak trees (Quercus suber). The understory has a herbaceous layer comprised of C3 species with an annual life cycle where green-up occurs in the fall and senescence in the dry season [25]. Additionally, a layer of shrubs exhibit drought resistant traits that reflects their phenological patterns [26]. This diverse mixture of vegetation has starkly different responses to climate.

The growing season of grasses is governed by precipitation [27]; the onset and duration of the dry summer periods can cause shifts in life-cycle timing [28]. Shrubs provide the understory more resilience in warm temperatures and low rainfall [28]. However, the copying mechanisms of Mediterranean shrubs to altered patterns of precipitation can be inconsistent, either advancing or delaying flowering [27]. Furthermore, the timing of the understory phenological events affects the interactions between tree canopy and understory [22,29,30], which are particularly important in the ability of these ecosystems as a carbon sink [31].

Thus, the main goal of this work is to verify the applicability of digital photography and GCC as a tool to monitor understory vegetation in Mediterranean cork oak woodlands. The first specific objective is to understand how well the method correlates to other vegetation indices typically used to assess changes in green biomass. In order to do so, vegetation was organized into three plant functional types (PFTs) based on their response to climate: trees, shrubs, and grasses [7]. The second specific objective is to investigate if grasses, cistus (Cistus salviifolius), and ulex (Ulex airensis) exhibit different relationships between GCC and spectral vegetation indices. The third and final specific objective is to identify to what extent and time lag climatic variables, namely temperature and precipitation, influence GCC variability.

2. Materials and Methods

2.1. Site

The study site is located in Central Portugal (latitude 39°8′18.35″ N, longitude 8°19′57.56″ W) where the climate is typically Mediterranean, with a hot and dry summer followed by a rainy season that lasts until the early spring. The annual mean air temperature from 1981 to 2010 was 17.3 °C and the average precipitation was 652 mm [32].

The site is characterized by heterogeneous vegetation comprised of cork oak trees (Quercus suber L.) and an understory of shrub and grass species. The understory is mostly represented by the shrub species, cistus (Cistus salviifolius), and ulex (Ulex airensis). The herbaceous layer is dominated by C3 species, mainly grasses and legumes, and represents 39.2% of the soil cover, while the mixture of shrubs represents 27.6%. The remainder is accounted for by litter and bare soil [33].

2.2. Data

2.2.1. Climatic Variables

A tower of 22 meters was installed in 2009 as part of a wider research project and takes measurements of meteorological variables continuously. In this study, we considered temperature (CS215, Campbell Scientific, Inc., Logan, UT, USA) and precipitation (ARG100, Environmental Measurements Ltd., Gateshead, UK). The micrometeorological data were averaged at a time-step of 30 minutes using a datalogger (CR10X, Campbell Scientific) [34].

2.2.2. Digital Photographs and Reflectance Data

Data were collected with a monthly frequency from 2011 to 2018, except for shrubs, where measurements stop in March of 2013 due to an understory cut (see Figure S1 of the Supplementary Materials for a map of the locations). Measurements were taken when vegetation was fully exposed to sunlight, within two hours around solar noon and during days of clear sky, to minimize the influence of atmospheric conditions. Sample individuals of each vegetation type were selected randomly within the footprint of the tower: eight shrubs (four cistus and four ulex) and four plots of the herbaceous layer, defined by a 50 × 50 cm quadrat.

Photographs were taken at nadir view using a regular digital camera (Fujifilm FinePix S5600). Immediately after, reflected radiation was measured using a field spectroradiometer, Fieldspec 3 (ASD Inc., Boulder, CO, USA) [35], providing hyperspectral reflectance data in the range of 350–2500 nanometers. Measurements were taken at about 90 cm from the target with both the camera (no zoom) and the fieldspec. In regard to the fieldspec, a pistol grip was used, and calibration with a white reference (Spectralon panel, Labsphere, Inc., North Sutton, NH, USA) was utilized to allow conversion of the radiation signal into absolute reflectance (R). Each target was measured five times, and each reflectance value was the average of 25 continuous spectra measurements, covering a circular sampling area with a diameter of 44 cm (the field of view of the fiber optic cable is 25° degrees).

2.3. Vegetation Indices

2.3.1. Picture-Based Indices

All pictures were analyzed using the R-package Phenopix [36]. Phenopix allows the user to mark a selected area of the photograph as the region of interest (ROI). This area is computed to extract the value of the digital number (DN) contained in each pixel, as seen in Figure 1.

The DN values that were sampled were averaged over the entire ROI, resulting in a single digital number value for red (R), green (G), and blue (B) within any given digital number of a photograph. Richardson [2,37] describes with detail a very similar ROI method. The relative index of green DN is quantified by the green chromatic coordinates (GCC) equation [38]:

$$GCC=\frac{GreenDN}{RedDN+GreenDN+BlueDN}$$

(1)

where GCC is a positive value ≤ 1 that expresses the weight of green pixels in the total number of ROI pixels. The package offers options to filter out inconsistent data, but no filter was applied.

2.3.2. Reflectance-Based Indices

Reflectance measurements were used to calculate the Normalized Difference Vegetation Index (NDVI) [39] and the Normalized Difference Water Index (NDWI) [40]. NDVI is sensitive to the structure and amount of vegetation, while NDWI uses the short-wave infrared region (SWIR) of the spectrum, which has been shown to be sensitive to radiation absorption by water. Their equations are shown in

Table 1. Spectral vegetation indices and their corresponding equations.

Vegetation Index	Formula
Normalized Difference Vegetation Index	$NDVI=\frac{R_{800nm}-R_{670nm}}{R_{800nm}+R_{670nm}}$	(2)
Normalized Difference Water Index	$NDWI=\frac{R_{860nm}-R_{1240nm}}{R_{860nm}+R_{1240nm}}$	(3)

.

These spectral vegetation indices (spectral VIs) were selected due to their suitability as a proxy for green biomass and vegetation structure. In addition, they have been proven successful by other studies to track vegetation changes at the site [34].

2.4. Data Analysis

2.4.1. Vegetation Indices

Spectral VIs and GCC were aggregated by fieldwork day: NDVI and NDWI were averaged with an arithmetic mean and GCC with a weighted mean. The weight is the ratio of green pixels to total pixels of a sample individual.

Statistical correlations between GCC and spectral VIs required R package PerformanceAnalytics [41] and were calculated with Pearson coefficients (r). The results were several matrixes of bivariate scatterplots with a regression line, a histogram describing frequency, and the corresponding r (see Figures S2–S4 in the Supplementary Material).

The first modelling approach was used to fit six simple linear regressions using base R software functions, where GCC was the dependent variable. The regressions were defined for each study group—grasses, cistus, and ulex—and for each spectral VI (NDVI/NDWI) in order to determine the best predictor of GCC for each group.

In the second approach, we considered the PFT effect on the relationship between GCC and spectral VIs by building four multiple linear regression models with a dummy variable. In two of these models, the dummy indicates whether data belong to grasses or shrubs, while in the other two, only shrub data were used, and the dummy indicates which species of shrubs the data belong to—cistus or ulex. The independent variables in these four models are as follows: NDVI, NDWI, and PFT or PFT _shrubs, and the term NDVI × PFT or NDVI × PFT _shrubs refers to an interaction effect.

2.4.2. Climate Variables

Precipitation and temperature were averaged at different time-lags calculated backwards from measurement day: 0, 30, and 90 days, selected for their correlation with the spectral indices, as described by Cerasoli et al. [34].

Three multiple linear regressions—one for each vegetation group—were tested in order to understand the effect of meteorological variables on GCC. Each model includes one variable for temperature: either average of the minimum, mean, or maximum daily temperature, and one variable for total precipitation. Hence, 12 different predictor variables were obtained and selected based on the value of their correlation (carried out as described in Section 2.4.1.) with GCC.

3. Results

3.1. Climate

The monthly mean temperature follows a fairly consistent seasonal pattern (Figure 2). Additionally, every year, a tendency for increased temperatures in the warmest months was observed except for one year, 2014, in which the monthly mean and temperature amplitude decreased.

Monthly total precipitation is generally more concentrated in autumn and winter than spring—notice the vertical dotted lines in Figure 2. Summer records almost no precipitation. These observations are persistent, but there are differences in the amount and timing of precipitation. For example, in the autumn of 2011, the rainiest month was December, while in 2012, it was November.

The driest year was 2018 (291 mm) and the rainiest was 2014 (876 mm). It is important to note that the year of 2018 is incomplete, with no data after October. However, it is likely it would still remain a relatively dry year, because the total precipitation up the point of the measurements is only 291mm. There was also a severe drought in the growing season of 2012.

3.2. Time Series

There is a very good match between GCC, NDVI, and NDWI time series of the herbaceous layer, as shown on Figure 3a. Similarly, shrubs also show very encouraging results (Figure 3b,c).

Grasses displayed a higher seasonality than shrubs, as every year at the start of the rainy season, there was an increase in the expression of the VIs until May, and then a decrease, reaching the lowest values in summer. This is related to plant annual growth activity, which is lower in the driest months of the year and increases when rain occurs, particularly in the herbaceous layer. The growing season of 2012, in particular, shows a distinctively lower peak when compared with other years.

The temporal series of cistus showed a similar trend; however, in ulex, seasonal trends were not as evident. In the months after October, greater GCC and NDVI values were observed because the weather conditions are favorable for vegetative development. Correspondingly, values for December, January, and February are frequently the highest, because the plant has fully grown and the canopy has most likely hit its fullest vegetative growth. Despite this, the temporal series of ulex stays relatively constant in all VIs throughout the year.

In terms of inter-annual variability, there is no obvious differences, which suggests that the relationship between GCC and spectral VIs is very similar between each year (see Figure S5 in the Supplementary Materials).

3.3. Association among Indices

Table 2. Equations ($y=b_1{x}_1+b_0$) ¹ and r² of the linear regressions for grasses, cistus, and ulex.

	Simple Linear Regressions
	NDVI
	Equation	r	r2	p-value
grasses	y = 16x1 + 0.30	0.945	0.89	<0.05
cistus	y = 0.20x1 + 0.24	0.793	0.62	<0.05
ulex	y = 0.24x1 + 0.24	0.564	0.30	<0.05
	NDWI
	Equation	r	r2	p-value
grasses	y = 0.45x1 + 0.43	0.917	0.84	<0.05
cistus	y = 0.41x1 + 0.40	0.885	0.78	<0.05
ulex	y = 0.51x1 + 0.42	0.754	0.55	<0.05

¹ where $y=GCC;x_1=NDVIorNDWI$.

indicates that the shrubs’ time series of GCC, NDVI, and NDWI are not as highly correlated as grasses, especially ulex. However, it is possible to see that there are higher correlations within each group (grasses, cistus, and ulex) than between each group. This cohesiveness stems from the similarity of the structural changes that are being depicted by all indices. Furthermore, shrubs are more correlated with each other than with grasses.

Shrubs show better correlation between GCC and NDWI than NDVI. However, GCC and NDVI are still significantly related. Cistus displays some seasonality but ulex is very constant throughout the growing season (see Figure 3). There is one exception in 2012, when the trend decreases across all indices and shrub species.

Grasses show a better coefficient of determination for GCC vs. NDVI (r² = 0.89) than GCC vs. NDWI (r² = 0.84). However, shrubs show better results in the regression with NDWI than NDVI. We then tested if the relationship between GCC vs. NDVI and GCC vs. NDWI was similar across all groups, searching if one model could explain the observed variability. To define the groups, we created dummy variables (PFT and PFT _shrubs), where PFT separates grasses from shrubs and PFT _shrubs separates cistus from ulex.

Table 3. Model ($y=b_1{x}_1+b_2{x}_2+b_3{x}_3+b_0$) ¹ where plant functional type (PFT) considers two groups—grasses and shrubs (cistus and ulex are classed together as one group)—and assumes interaction between predictor variables NDVI/NDWI and PFT.

Models Separating Vegetation by PFT
NDVI
	Coefficients:
	Estimate	Std. Error	t value	Pr (>|t|)	Residual standard error: 0.01885 on 124 DF
(Intercept)	0.296693	0.005656	52.452	< 2 × 10−16 ***
NDVI	0.161242	0.013046	12.360	< 2 × 10−16 ***	Adjusted R-squared: 0.7274
PFT	−0.073573	0.017942	−4.101	7.4 × 10−5 ***	F-statistic: 113.9
NDVI×PFT	0.080054	0.027463	2.915	0.00422 **	p-value = 0.05
NDWI
	Coefficients:
	Estimate	Std. Error	t value	Pr (>|t|)	Residual standard error: 0.01562 on 124 DF
(Intercept)	0.434724	0.005529	78.628	<2 × 10−16 ***
NDWI	0.445858	0.030808	14.472	<2 × 10−16 ***	Adjusted R-squared: 0.8128
PFT	−0.018052	0.006100	−2.960	0.00369 **	F-statistic: 184.8
NDWI×PFT	0.060766	0.046307	1.312	0.19186	p-value: < 2.2 × 10−16

¹ where $y=GCC;x_1=NDVIorNDWI;x_2=PFT;x_3=NDVI\times PFTorNDWI\times PFT$. Note: three asterisks (*) indicates p < 0.001 and two indicates p < 0.01.

refers to the models where observations were divided by PFT—grasses and shrubs.

The variables NDVI, PFT, and NDVI × PFT are all significant predictors of GCC. The same can be said for NDWI, but the interaction effect is not significant (NDWI × PFT). Hence, we conclude that PFT separation is significant, and that one single model without the variable PFT for either NDVI or NDWI cannot describe the observed variability for grasses and shrubs.

We tested further to see if separating shrubs by species, as we initially did in the simple linear regressions, was appropriate.

Table 4. Model ($y=b_1{x}_1+b_2{x}_2+b_3{x}_3+b_0$) ¹ that separates shrubs by species and assumes interaction between predictor variables NDVI/NDWI and shrubs species (cistus or ulex).

Models with Shrub Distinction
NDVI
	Coefficients:
	Estimate	Std. Error	t value	Pr (>|t|)	Residual standard error: 0.02009 on 62 DF
(Intercept)	0.2395414	0.0205342	11.665	<2 × 10−16 ***
NDVI	0.1972340	0.0299884	6.577	1.15 × 10−8 ***	Adjusted R-squared: 0.6406
PFT shrubs	0.0007734	0.0461891	0.017	0.987	F-statistic: 39.62
NDVI×PFT shrubs	0.0386738	0.0645447	0.599	0.551	p-value: 2.01× 10−14
NDWI
	Coefficients:
	Estimate	Std. Error	t value	Pr (>|t|)	Residual standard error: 0.01577 on 62 DF
(Intercept)	0.402185	0.004186	96.084	< 2 × 10−16 ***
NDWI	0.414499	0.044312	9.354	1.84 × 10−13 ***	Adjusted R-squared: 0.7787
PFT shrubs	0.022588	0.005405	4.179	9.34 × 10−5 ***	F-statistic: 77.25
NDWI×PFT shrubs	0.098773	0.085586	1.154	0.253	p-value: < 2.2 × 10−16

¹ where $y=GCC;x_1=NDVIorNDWI;x_2=PF{T}_{shrubs};x_3=NDVI\times PF{T}_{shrubs}orNDWI\times PF{T}_{shrubs}$. Note: significance (*) as labelled in Table 3.

shows the model that analyses if each shrub species is not significantly different to justify the use of a single model with the PFT variable.

The model (Table 4) GCC vs. NDVI has an overall p-value = 2.01 × 10⁻¹⁴ and an adjusted R² = 0.64. PFT _shrubs is not significant (p-value < 0.05); hence, grouping shrub species together appears to not add relevant information for the model. Regardless, fitting an equation to all three groups shows higher R-squared for shrubs than the individual regressions (r² = 0.62 and 0.30 for cistus and ulex, respectively). When using NDWI as a proxy for GCC, results are more ambiguous. The separation of shrubs shows a slightly lower agreement (R² = 0.78 in Table 4) compared to the model where shrubs are together (R² = 0.81 in Table 3), but the p-value of PFT _shrubs is more significant than that of PFT. Figure 4 shows the graphical representation of these four models.

There is a clear difference in the slope of GCC vs. NDVI when grouping by PFT and quite some overlap when separating by shrubs (Figure 4a). For GCC vs. NDWI, the separation between cistus and ulex data points is greater than in GCC vs. NDVI (Figure 4b). Furthermore, the data points cluster closer to the regression lines in GCC vs. NDWI, regardless of the organization of samples. This can be traced back to the residual standard errors in Table 3 and Table 4 (0.01562 and 0.01577, respectively). It indicates a lower difference between the observed and the predicted values for GCC vs. NDWI when compared to GCC vs. NDVI.

3.4. Association between GCC and Climate Variables

When it comes to the climate signal in GCC, the accumulated precipitation of 90 days (Prec90) is always the highest correlated precipitation variable across all groups (Figure S6 of the Supplementary Materials). The time lag factor largely affects grasses more than shrubs (Figures S7–S9).

Mean temperature is more variable. In grasses and cistus, the highest correlation is with mean maximum temperature at different time lags, indicating that grasses react quicker than cistus to temperature changes. Shrub ulex did not show a clearly stronger correlation with any of the measures of mean temperature tested.

The highest correlations of mean temperature and total precipitation grouped at different time windows (see Section 2.4) were selected for the multiple linear regressions that are shown in

Table 5. Multiple linear regression between GCC, mean temperature, and total precipitation.

Grasses	Selected Variable	t Value	Pr (>|t|)
	Tmax30 (r = −0.66)	−3.793	0.000383 ***	Adjusted R-squared: 0.5419
	Prec90 (r = 0.66)	3.813	0.000360 ***	p-value: 3.882 × 10−10
cistus	Selected Variable
	Tmax90 (r = −0.63)	−3.393	0.00196 ***	Adjusted R-squared: 0.3786
	Prec90 (r = 0.44)	1.063	0.29611	p-value: 0.0003021
ulex	Selected Variable
	Tmin90 (r = −0.37)	−1.739	0.0924	Adjusted R-squared: 0.1215
	Prec90 (r = 0.31)	1.177	0.2483	p-value: 0.05441

Note: Significance of asterisks (*) as labelled in Table 3 and the dot indicates p < 0.1.

.

The results for grasses indicated that Tmax30 (p-value = 0.000383) and Prec90 (p-value = 0.000360) were both significant in explaining GCC variability. The results for cistus showed that only Tmax90 (p-value = 0.00196) was significant, while precipitation was not. Additionally, for ulex, neither Tmin90 nor Prec90 was significant at α = 0.05.

4. Discussion

4.1. The Applicability of GCC

4.1.1. Time Series Analysis

Despite relatively high variability, the picture-based index GCC was able to represent seasonal trends in understory vegetation, as indicated by how well it correlates with spectral VIs. Furthermore, these results are in agreement with the recent work on grasses of other authors [5,42]. Spectral VIs have the ability to capture the variability of several characteristics. NDVI has higher sensitivity to changes in density of green leaves [40,43] and, thus, an increased ability to follow variations at different stages of plant development. Photographs are able to capture a large range of visual changes in vegetation, but GCC focuses only on the green channel. Thus, the lower amount of information associated with the GCC index could be related to the lower amplitude of values in the time series (Figure 3).

Changes in photosynthetic structures are connected with seasonality and should account for GCC variability, since it is chlorophyll pigments that mostly confer vegetation color. However, the relationship between foliar area and GCC is not straightforward [11,44], and further work to understand if GCC is more sensitive to LAI or chlorophyll content is needed.

The shrub ulex showed decreased GCC in March/May likely related to flowering despite the consistently high leaf cover observed in photographs. Shrub cistus has higher GCC attributed to pictures with a lesser amount of small sized leaves with more chlorophyll per unit of dry mass [26]. In cistus, GCC appears more sensitive to foliar area, but in ulex, it is more sensitive to chlorophyll content.

Peak chlorophyll content occurs more towards the end of the growing season, and if GCC captures chlorophyll content well, this could be associated with a reduction of signal in the green channel as leaves begin to dry and brown [44]. Contrarily, we observed higher values in the beginning of the growing season, suggesting that GCC is sensitive to leaf area. Despite this ambiguity in what GCC is truly measuring at a certain time in the phenology of vegetation, the use of GCC to study the heterogenous vegetation of Mediterranean woodlands provided a rather similar overall trend to the more conventional NDVI and NDWI.

The close correlations between GCC and NDVI/NDWI in grasses and cistus (Table 2) suggest that the observed seasonal variability in values is due to a similarity in inherent biophysical properties and response to climate factors [42,45]. Grasses and cistus are more similar with each other than ulex, sharing the strategy to renew photosynthetic structures frequently and adapting to rainfall.

Water availability is an important factor in timing the start of the growing season. It is essential to enable germination and biological processes resulting in vegetation growth [46]. At the end of the growing season, as senescence progresses, tissues loose water progressively. Thus, in annual species such as grasses, a high correlation between GCC vs. NDWI is not surprising. In semi-deciduous evergreen shrubs, the good agreement observed between GCC and NDWI is congruent with the premise that density of leafy structures is also very related to water availability [47,48], despite a slow reaction to rainfall due to a more effective water use by the deeper root system.

4.1.2. Method Usefulness

The ease with which good quality photographs can be taken nowadays is indisputable, and the use of cameras in phenological research poses a very quick and relatively cheap solution to recording phenological changes [2].

The Mediterranean woodlands are a fairly static ecosystem in terms of overall carbon stock [24]. Shrub species on site have almost no leaf renewal regardless of whether there is more water availability (see Figure 2 and Figure 3) [26,49]. We captured measurements during a drought year (2012) that did not affect shrubs greatly, as they stayed relatively intact throughout the growing season. Thus, the authors believe that a low amount of sample individuals is not a major concern when the signal is quite constant. However, for grasses, which are definitely more impacted by water, applying the method as described in the paper implies a priori knowledge of the study site [50].

The method has the potential to grow with little effort overtime into a comprehensive database that provides a visual archive and does not require a trained individual to be maintained. It can be capable of offering a temporal resolution tailored to necessity, but most importantly, allows for an efficient way to access highly heterogeneous understory phenology. The combination of these characteristics is a key feature of the method, since often times, the most dispendious and time-consuming part of small -scale phenological research is data collection in the field [51,52].

On-the-ground monitoring offers the best spatial resolution for heterogenous areas such as the understory of the Mediterranean woodlands, but it also has limitations. To the knowledge of the authors, issues with lighting conditions, low temporal resolution, and area coverage have not yet been fully addressed, and few studies express concerns in this area [2,53]. For this reason, such a study should perhaps be the focus of future research since it could either affirm the use of GCC and digital photographs by helping to explain variability in the data, or fully define the extent to which it could be used.

In this work, the time series used to derive GCC was very small compared to other studies [5,10] while still yielding useful results that are in agreement with other phenological observations in similar environments [26,50]. Hence, it goes in favor of the concept of using digital photography as a means to track phenology, enlarging the range of ecosystems where the method has shown success. However, more importantly, it provides information on the ability of the method to deliver usable results with limited data.

4.2. Differences among Grasses, Cistus, and Ulex

The multiple linear regression models are clearly different (Table 3 and Table 4). GCC of grasses has a better correlation with NDVI, and GCC of shrubs correlates better with NDWI. This difference extends further into how data should be organized since our results show that separating by PFT is essential when using GCC as a proxy for NDVI. However, for NDWI, the most effective way to organize vegetation would be to consider the two shrub species separately, because shrubs are sufficiently different in terms of water usage (Figure 4).

In practice, to estimate NDVI from GCC it is advisable not to use photographs or a single value that accounts for all three groups, since one single simple model does not explain the full variability of data. In the case of NDWI, photographs should not have both cistus and ulex, as the values should be treated separately, i.e., using individual photographs. Despite that, the models GCC vs. NDWI have a higher predictive ability than GCC vs. NDVI, so it is preferable to use NDWI when attempting to fit a model with more than one vegetation type.

4.3. Responsiveness of GCC to Climate

Large differences were found in the relationship between GCC and climate variables and their time lag according to the vegetation group, which emphasizes the need to carefully monitor each component of the understory. Furthermore, the magnitude of response to climatic variables is different according to species.

Grasses display the greatest changes in GCC during their annual life cycle, and both Tmax30 and Prec90 influence GCC values (p-value = 3.882 × 10⁻¹⁰). As previously demonstrated in many other studies [22,30,49], grasses are largely influenced by rainfall. This dynamic is more evident during the extreme drought of 2012, when precipitation was as low as 132 mm in the January–August period and all indices (GCC and VIs) had consistently low values (see Figure 2 and Figure 3).

The higher importance of Prec for grasses than shrubs is in agreement with the biophysical properties of each group. Grasses have an annual life cycle that adjusts every year according to the timing of environmental conditions [28]. Shrubs are more drought-resistant species that invest in the root system [49], making them able to withstand longer periods of low water availability; thus, their GCC did not appear to be as influenced by precipitation.

Tmax30 is equally important but showed an opposite trend (r = −0.66). The results suggest that GCC can be reasonably predicted by Tmax30 and Prec90 (R² = 0.54); however, further analysis on the relative importance of each variable to GCC is necessary.

The semi-deciduous cistus was significantly influence only by Tmax90 (r = −0.63). On the other hand, the evergreen ulex showed no significant impact of the two most highly correlated climate variables (Tmin90 and Prec90).

Ulex always keeps its foliage and did not respond to the extreme drought in 2012, likely because it has deep roots and leaves modified into spines that prevent high transpiration. Cistus has a shallow root system that translates into some intra-annual seasonality that is similar to grasses, and accordingly, both have higher correlation with Tmax and Prec90. These differences could explain the observed lower amplitude of VIs of ulex when compared with grasses and cistus and also the progressively lower p-values of each multiple regression between GCC and climate variables (Table 5).

The higher correlation with higher time lags in shrubs agrees with a much less marked seasonal behavior than grasses and shrubs’ ability to endure unfavorable conditions since both species seek to maintain productive structures, despite different strategies.

Such differences are a reflection of distinct canopy dynamics. The foliar structure of cistus is not completely lost every season, even under prolonged summer drought conditions [54]. It is also visible in lower GCC values at the start of the rainy season that increase moderately about two months later. This slow response to autumn rains is supported by other authors with studies in the same location [26]. On the contrary, ulex always maintains a similar canopy density [34].

5. Conclusions

The complex dynamic between vegetation and climate in Mediterranean cork oak woodlands, namely, the understory, needs to be monitored at a fine resolution due to its sensitivity to climate change, and the GCC index allows for a cost-effective method to do so. However, more work is required to pinpoint what phenological changes can be represented by GCC.

The different seasonality patterns observed in grasses and shrubs can be accounted for by intrinsic biophysical properties of the vegetation and the various survival strategies of each vegetation group and do not represent an inability of the method to properly depict variations. The herbaceous layer was highly related to rainfall regime, while cistus and ulex, especially the latter, demonstrated less sensitivity to climatic conditions. Grasses were also impacted by temperature with a 90 days lag time. The shrub cistus displayed an intermediate response between grasses and ulex, which exhibited no impact of either temperature or precipitation. These responses were well represented by the GCC trend, despite the uncertainty of whether it is more sensitive to leaf area or chlorophyll content.

Thus, the findings of this study suggest that digital cameras might not substitute the need for more sophisticated devices such as field spectroradiometers and more comprehensive data. However, in many cases, they can be an alternative that offers a good compromise of cost and benefit that allows vegetation to be adequately characterized.