phenopix: a R package for image-based vegetation phenology Gianluca Filippaa,∗, Edoardo Cremonesea, Mirco Migliavaccab, Marta Galvagnoa, Matthias Forkelb, Lisa Wingatec, Enrico Tomellerid, Umberto Morra di Cellaa, Andrew D. Richardsone aEnvironmental Protection Agency of Aosta Valley, ARPA Valle dAosta, Climate Change Unit (Italy) bMax Planck Institute for Biogeochemistry, Department Biogeochemical Integration (Jena, Germany) cINRA, UMR ISPA (Bordeaux, France) dInstitute for Applied Remote Sensing, EURAC (Bolzano, Italy) eHarvard University, Department of Organismic and Evolutionary Biology (Cambridge, MA, USA) Abstract In this paper we extensively describe new software available as a R package that allows for the extraction of phenological information from time-lapse digital pho- tography of vegetation cover. The phenopix R package includes all steps in data processing. It enables the user to: draw a region of interest (ROI) on an image; extract red green and blue digital numbers (DN) from a seasonal series of im- ages; depict greenness index trajectories; fit a curve to the seasonal trajectories; extract relevant phenological thresholds (phenophases); extract phenophase un- certainties. The software capabilities are illustrated by analyzing one year of data from a se- lection of seven sites belonging to the PhenoCam network (phenocam.sr.unh.edu/), including an unmanaged subalpine grassland, a tropical grassland, a deciduous needle-leaf forest, three deciduous broad-leaf temperate forests and an evergreen needle-leaf forest. One of the novelties introduced by the package is the spa- tially explicit, pixel-based analysis, which potentially allows to extract within- ecosystem or within-individual variability of phenology. We examine the rela- ∗Corresponding author Email address: g.filippa@arpa.vda.it (Gianluca Filippa) Preprint submitted to AFM December 29, 2015 © 2016. This manuscript version is made available under the Elsevier user license http://www.elsevier.com/open-access/userlicense/1.0/ http://phenocam.sr.unh.edu/ http://ees.elsevier.com/agrformet/viewRCResults.aspx?pdf=1&docID=7591&rev=1&fileID=226875&msid={DD76000D-D1E4-4811-8A6D-2033ADD1F633} tionship between phenophases extracted by the traditional ROI-averaged and the novel pixel-based approaches, and further illustrate potential applications of pixel-based image analysis available in the phenopix R package. Keywords: image analysis, community ecology, pixel-based analysis, phenology 1. Introduction Traditional monitoring of plant phenology relies on direct human observa- tions of discrete phenological events, or phenophases, such as bud-burst, flow- ering, autumn decoloring, and leaf-fall (e.g. Lechowicz, 1984; Richardson et al., 2006; Galvagno et al., 2013; Migliavacca et al., 2008; Filippa et al., 2015). Such5 observations are typically made on a limited number of individual organisms, across a limited geographic area (i.e., often at a specific research site). On the other hand, satellite remote sensing allows observing land surface phenology on regional to global scales but has a limited representativeness for phenolog- ical changes at ecosystem or species-level (White and Nemani, 2006; Delbart10 et al., 2005; Busetto et al., 2010; Hufkens et al., 2012; Forkel et al., 2015). At an intermediate scale, near-surface remote sensing of phenology makes use of radio- metric instruments or imaging sensors. Near-surface remote sensing quantifies, at high temporal resolution, and with a flexible degree of spatial integration (i.e., the potential to look across the canopy as a whole, but at the same time15 focus on individual organisms), seasonal changes in the optical properties of veg- etated surfaces (e.g Jenkins et al., 2007; Richardson et al., 2007; Soudani et al., 2012). Recent studies have demonstrated the soundness of digital cameras as multi-channel imaging sensors (Richardson et al., 2009; Klosterman et al., 2014; Wingate et al., 2015; Migliavacca et al., 2011).20 In the recent literature, different approaches have been used to process dig- ital images of the vegetation canopy. Several researchers have exhaustively addressed the issues of data quality and data filtering in order to reduce noise in the seasonal trajectories of greenness (e.g. Sonnentag et al., 2012; Julitta 2 et al., 2014; Migliavacca et al., 2011). Other authors focused on curve fit-25 ting/smoothing methods for extracting dates from phenological time-series (i.e. phenophases, see e.g. Zhang et al. 2003; Beck et al. 2006; Gu et al. 2009; Elmore et al. 2012; Klosterman et al. 2014 and the greenbrown R package). The full exploitation of the existing plethora of methods and approaches would benefit from a comprehensive framework in which to compare them across ecosystem30 types and climate conditions. Phenological networks based on the analysis of digital images are growing world- wide. In the US, PhenoCam network includes over 200 sites that use a standard camera and follow a standard protocol (phenocam.sr.unh.edu/webcam/). In Europe, the EUROPHEN network consists of about 60 flux sites equipped with35 digital cameras (Wingate et al., 2015). Similar networks are also present in Aus- tralia (www.aceas.org.au, Brown et al. in preparation) and Asia (Nasahara and Nagai, 2015). The increasing use of repeated digital photography for phenologi- cal research highlights the need of easy-to-use, open source and flexible software tools for the processing of the images. Additionally, the existing networks are40 poorly coordinated in terms of camera types, settings and sampling protocols, thus presenting a big challenge to build a flexible tool capable of facing the large diversity of ecosystems, image quality and setups. In this paper we present a collection of functions packed in a software available as a R package (R Core Team, 2015), called phenopix (r-forge.r-project.org/projects/phenopix/).45 We will first show the main features of the package and then illustrate its ap- plication on a selection of sites belonging to the PhenoCam dataset. Lastly, a section focused on pixel-based analysis will: 1) examine the relationship between phenological thresholds (phenophases) extracted from the average seasonal tra- jectory of greenness over a region of interest (ROI-averaged approach) and from50 each pixel of such region (pixel-based approach); 2) illustrate potential appli- cations of pixel-based image analysis to discriminate between subtly different phenological seasonal trajectories within the same image scene. 3 http://r-forge.r-project.org/projects/phenopix/ 2. PhenoCam sites The study sites used to illustrate the functionality of the phenopix pack-55 age belong to the PhenoCam network (phenocam.sr.unh.edu/webcam/) and are illustrated in table 1 and fig.1. Table 1: Main characteristics of the selected PhenoCam sites. PFT: plant functional type, GRA: grassland, DBF: deciduous broadleaf forest, DNF: deciduous needle-leaf forest, ENF: evergreen needle-leaf forest Site ID Coords (deg) Elev (m) PFT Camera type Dominant species torgnon-ND 45.8N 7.6E 2160 GRA Nikon D5000 Nardus stricta kamuela 20.0N -155.7E 850 GRA Stardot C3/C4 grasses bartlett 44.1N -71.3E 268 DBF Axis 211 Acer rubrum; Fagus grandifolia harvard 42.5N -72.1E 340 DBF Stardot Quercus rubra; Acer rubrum mammothcave 37.2N -86.1E 226 DBF Stardot Quercus sp.; Carya sp. torgnon-LD 45.8N 7.6E 2091 DNF Nikon D5000 Larix decidua harvardhemlock 42.5N -72.2E 345 ENF Stardot Tsuga canadensis 3. Main functions The typical work-flow of the phenopix package is summarized in the flowchart shown in fig. 2. First, one (or more) region(s) of interest (ROI) is (are) cho-60 sen, then digital color numbers are extracted from the ROI of each image, and processed to obtain a seasonal time series. After filtering the time series, data are fitted with either a double logistic equation or a smoothing curve, on which phenological thresholds (phenophases) are extracted. Finally, uncertainty of the fit and of phenophases can be computed.65 3.1. Regions of interest (ROIs) The scene of the picture rarely includes only the targeted vegetation canopy, thus the user will want to choose a particular region within the scene for analysis. Even more often, more than one region may be of interest, for example in a 4 mixed forest one might independently analyze different deciduous species and70 evergreen trees (e.g. Ahrends et al., 2008). The function DrawROI() allows the user to draw one or more regions of interest on-screen, using the mouse cursor on a chosen reference picture. 3.2. Extract vegetation indices From the digital color values of each image the green chromatic coordinate75 (GCC ) is computed. GCC is a vegetation index derived from photographic image and quantifies the greenness relative to the total brightness. GCC is computed as follows: GCC = GDN RDN + GDN + BDN (1) where GDN , RDN , BDN are the green, red and blue digital numbers, re- spectively (Gillespie et al., 1987). Similarly, chromatic coordinates of red and80 blue (RCC and BCC ) are also computed. Several indices based on RGB colors have been developed in the last years, including for example the green excess index (GEI) (Woebbecke et al., 1995; Mizunuma et al., 2013). Some authors also used a combination of GCC and RCC to extract autumn phenophases (e.g. Klosterman et al., 2014). For simplicity, all subsequent analyses will be focused85 on GCC but phenopix allows the analysis of the whole variety of color indices, including the computation of new ones. The function extractVIs() extracts raw red, green and blue digital numbers from each pixel in the ROI, and computes color chromatic coordinates as in eq.1. Vegetation indices can be computed on ROIs with 2 approaches (fig. 2):90 1) the ROI-averaged approach: color chromatic coordinates are computed for each pixel and then averaged over the whole ROI, and (2) the pixel-based ap- proach, where each pixel belonging to the ROI is analyzed separately (section 5). The procedure is repeated for each image in the archive. A time-stamp is retrieved from the file name of the image and a time series of the computed95 indices is returned. Specific rules must be followed in naming the image files, 5 and the reader is referred to the package help pages for more details. 3.3. Data filtering Data retrieved from images often need robust methods for filtering the time100 series. Bad weather conditions, low illumination and dirty lenses are among the most common issues that determine noise in the time series of vegetation indices. Previous studies (e.g. Sonnentag et al., 2012; Julitta et al., 2014; Migli- avacca et al., 2011; Papale et al., 2006) provide the background upon which the filtering techniques implemented in phenopix were chosen. We designed a105 function autoFilter() based on 5 different approaches. (1) A filter called night which removes GCC records lower than a certain thresh- old, by default 0.2. This filters removes night data or images with very scarce illumination (e.g. early morning or sunset images, cloudy days, etc.). (2) A filter called blue, based on blue chromatic coordinate (BCC ) (Julitta et al.,110 2014). First, hourly BCC data is aggregated at a daily resolution and daily means and standard deviations are computed. A quantile (by default the 5th percentile) is then calculated on standard deviations. This threshold is then added and subtracted from the daily mean BCC to generate a seasonal enve- lope. Raw data falling outside this envelope are discarded. The blue filter was115 designed to remove images with dominant clouds and/or snow on the canopy, because the blue chromatic coordinate was found to be very sensitive to such conditions. (3) A filter called mad, following the method of Papale et al. (2006), an outlier detection based on the double-differenced time series, using the median of ab-120 solute deviation about the median (MAD) that is a robust outlier estimator. (4) A filter called spline, following the method of Migliavacca et al. (2011), based on recursive spline smoothing and residual computation followed by removal of outliers falling outside a given residual envelope. (5) A filter called max, following the method of Sonnentag et al. (2012), based125 on the identification of the 90th percentile values in three-days moving windows. 6 Filters can be applied alone or in sequence so that the user can choose the best combination of filters that suits the image archive to be processed. Addition- ally, each filter has its own discarding criteria that may be tuned by the user according to the quality of the input GCC time series. The default behavior of130 the filtering function is to use a sequence of night, spline and max filters. In general, this sequence will be effective enough to properly filter the GCC time series. The user is advised to apply the max filter to effectively minimize the impact of changes in scene illumination (Sonnentag et al., 2012). 3.4. Fit a curve to the GCC seasonal course and extract phenophases135 The extraction of phenophases is done in two steps (following the approach of Forkel et al. 2015; Klosterman et al. 2014). First, a curve is fitted to the GCC seasonal curve to reduce the influence of single observation and to better capture the seasonal behaviour. In a second step different extraction methods can be used to retrieve phenophase metrics from the fitted seasonal curve. We selected140 four different double logistic equations to be included in the package (table 2). The equations differ in the number of parameters to be optimized and hence in the flexibility of the fitting curves. For a thorough examination of the fittings and explanation of curve parameters, see the correspondent publications. In addition to the double logistic equations, an approach based on a smoothed145 cubic spline is also available. There are cases, e.g. with very noisy time series or weak signal (low seasonal amplitude in greenness) where the double logistic fits fail. In such cases, the spline method is the only possibility to extract a seasonal trajectory. Table 2: Double logistic model options available in phenopix, along with their option names as implemented in the package. The first two methods were first implemented in the greenbrown R package and described in Forkel et al. (2015) Equation Reference Option name f(t) = mn + (mx - mn)·( 1 1+e(−rsp∗(t−sos)) + 1 1+e(rau∗(t−eos)) ) Beck et al. (2006) ’beck’ f(t) = m1 + (m2 − m7t) · ( 1 1+e (m′ 3 −t)/m′ 4 − 1 1+e (m′ 5 −t)/m′ 6 ) Elmore et al. (2012) ’elmore’ f(t) = (a1t + b1) + (a2t 2 + b2t + c) · ( 1 [1+q1e −h1(t−n1)]v1 − 1 [1+q2e −h2(t−n2)]v2 ) Klosterman et al. (2014) ’klosterman’ f(t) = y0 + a1 [1+e−(t−t01)/b1 ]c1 − a2 [1+e−(t−t02)/b2 ]c2 Gu et al. (2009) ’gu’ 7 After curve fitting, a suite of different approaches is available to extract150 phenophases, as summarized in table 3, and illustrated in figure 3. The phenophase extraction methods lead to the definition of dates that may have markedly dif- ferent ecological meanings and be useful for several applications in the field of environmental sciences. However the discussion of such meanings is beyond the objectives of this presentation paper, and the reader is referred to the publica-155 tions listed in table 2 for further information. Here it is important to point out that the availability of very simple methods (e.g. trs ) and more sophisticated ones (e.g. Klosterman method) allow the user to choose the extraction method better suited to the phenological trajectory of the investigated ecosystem and/or based on the quality of the input data.160 Table 3: Methods for phenophase extraction available in phenopix. sos: start of season, eos: end of season, los: length of season, pop: peak of season position, peak: maximum seasonal GCC , mgs: mean growing season GCC , msp: mean spring GCC , mau: mean autumn GCC , rsp: rate of spring greenup, rau: rate of autumn senescence, UD: upturn date, SD: stabilization date, DD: downturn date, RD: recession date, prr: peak recovery rate, psr: peak senescence rate Name Phenophases Description trs sos, eos, los, pop, mgs, peak, msp, mau Based on a user-defined thresh- old of seasonal development of GCC derivatives as TRS plus rsp and rau Based on local extremes in the first derivative klosterman Greenup, Matu- rity, Senescence, Dormancy Based on local extremes in the rate of change of curvature k (Kline, 1998) gu UD, SD, DD, RD, prr, psr Based on a combination of lo- cal maxima in the first derivative (Gu et al., 2009) The combination of four fitting methods plus the cubic spline smoothing 8 with the four phenophase extraction methods produces 20 different available combinations in output from a single yearly GCC time series. A specific function (greenExplore()) is designed to give an overview of all fittings and phases for a given dataset (fig. 4). The RMSE for each fit is also shown, so that the user can165 use it as a criterion for fit selection. Once the fit is chosen, by examining the plot along a row the user can evaluate which phenophase method is more suitable to the GCC time series. The possibility to combine different fitting and phenophase methods to a GCC time series provides a framework in which to compare in detail the different methods currently in use. This need was highlighted by170 Keenan et al. (2014) in a recent attempt to link forest phenology as described by time-lapse photography and physiology as described by traditional plant trait measurements. In addition to the above described methods, phenopix also implements a phenophase extraction method (function PhenoBP()) based on linear piecewise175 regression and correspondent break points in the time series (Wingate et al., 2015). This method was designed to accommodate multiple greening peaks during the same season, which is typical, for example, of water limited ecosys- tems such as Kamuela (fig 5) or managed grasslands and croplands. 4. Estimation of the uncertainty180 Traditionally, the analysis of digital images for phenology has rarely included the estimation of uncertainty on phenophase extraction, despite its paramount importance (e.g. for the evaluation of method robustness and as input data for the optimization of process-based or land surface models, where phenology is currently poorly represented (Richardson et al., 2013)).185 The phenopix package provides two approaches for the estimation of uncer- tainty, one for the fitted double logistic equations and one for the smoothing spline method. For the fitted equations the residuals between the model fit and the observed data are used to generate random noise and fitting is ap- plied recursively to randomly-noised original data. This procedure results in190 9 an ensemble of curves and curve parameters on which phenophase extraction is performed, thus providing an uncertainty estimate of phenophases. For the spline smoothing method, the uncertainty is estimated by combining a random noise as for the fitted equations with multiple changes in the spline’s degrees of freedom, to account for the arbitrariness in their choice. In particular, spline195 degrees of freedom are set at 5% of the time series length and then let vary re- cursively between 1% and 5% of the time series length to generate the ensamble of smoothing curves on which the uncertainty is estimated. In addition to the residuals method, for the fitted equations only, we are also developing a method based on the hessian matrix of the curve parameters.200 Figure 6 shows the uncertainty analysis for Harvard deciduous 2013 GCC sea- sonal trajectory, with 500 replicates. The fitting method applied in this case was Klosterman. By default, the confidence interval is represented by the 10th and 90th percentile of the distribution of all extracted phenophases, but other op- tions are available. In any case, the user can access the uncertainty data tables205 of (1) curve equation parameters and (2) extracted phenophases to customize the computation of the uncertainty envelope. 5. Phenopix application: spatial analysis Based on an analysis run on all sites included in this paper except for Ka- muela (i.e. six year-sites), we present in this paragraph the application of the210 pixel based analysis to investigate the spatial variability of phenology within a ROI (pixel-based approach). To date, few studies have explored the possibility of analyzing a set of images pixel by pixel rather than averaging color values on the overall region of interest (Julitta et al., 2014; Ide and Oguma, 2013). Therefore, an emergent question215 is how the ROI-averaged analysis is representative of all ROI pixels in terms of extracted phenophases? And, is this relationship consistent across fitting meth- ods and ecosystem types? The phenopix package allows the user to perform pixel based image analysis and provides a number of functions to display and 10 analyze the results. It must be noted that pixel-by-pixel analysis is computa-220 tionally intense. We report in table 4 the computation times required to analyze a ROI containing 10000 pixels on a seasonal time series of 5000 images with all fits available in phenopix for 3 different computers, using one single processor. Computation time may be as high as 70 hours, with very low values for spline smoothing and higher for data filtering and fitted equations. However, parallel225 computation available in the spatial functions of phenopix enable considerably reduced computation times compared to those reported in table 4. Table 4: Estimated computation time (hours) required to complete different steps of the pixel- based analysis processing 10000 pixels from a seasonal time series of 5000 images using one single processor on three different computers step 3.07 GHz 2.60 GHz 2.13 GHz filtering 28 49 52 spline 0.02 0.02 0.04 beck 36 63 69 elmore 0.7 0.8 1.5 klosterman 10 16 19 gu 22 37 42 We evaluate the relationship between approaches (i.e. ROI-average vs pixel- based) and fitting/phenophase methods across sites in fig. 7. Two different fitting methods (namely Klosterman and spline) applied to each pixel of the230 ROIs and then averaged are in good agreement with each other (fig. 7, first above diagonal panel), with a slightly higher correlation for spring than for au- tumn phases. The error bars for this relationship suggest a higher variability for autumn than for spring phenophases (first below diagonal panel). When comparing the ROI-averaged and the pixel-based approach, the best relation-235 ship is found with the Klosterman fitting method applied to both approaches, with correlation coefficients of 0.97 and 0.75 for spring and autumn phenology, respectively. However, the relationship between spline fit applied pixel-by-pixel 11 and Klosterman fit applied to the average ROI also shows very good agreement. In contrast, the spline method applied to the average ROI (last column in fig.240 7) has a consistently weaker relationship with all other methods, resulting in earlier autumn phases and a much higher variability. Across phenophase methods (different symbols in fig. 7), there is no evidence of a systematic bias between the two approaches for any particular phase extraction method. Across ecosystem types, needle-leaf forest (ENF, i.e. harvardhemlock245 site) shows a consistent departure from the 1:1 line even in the best relation- ships (third panel from left in the first row of fig. 7), probably due to the lower signal to noise ratio in the seasonal trajectory of evergreen trees. However, the most striking difference is the consistently lower correlation for autumn than for spring phases, suggesting that in the same ecosystem and possibly also within250 a single species autumn phases are more variable than spring phases. In summary, the faster spline smoothing applied pixel-by-pixel leads to the iden- tification of phenophases in substantial agreement with the more computationally- intense Klosterman fit. Additionally, being more flexible, the spline fitting leads to a much lower number of unfitted pixels compared to the Klosterman fit.255 Hence, from this analysis, spline smoothing is preferred over curve fitting for pixel-by-pixel analysis. For the ROI-averaged approach, it is strongly prefer- able to perform curve fitting over spline smoothing, provided that the GCC seasonal trajectory does not show multiple peaks, related to e.g. water stress or management practices.260 One interesting application of the pixel-based analysis stems from the ex- amination of the frequency distribution of a given phenophase across all pixels. When a phenophase shows a bi- (or multi-) modality and a consistent spatial distribution, pixels may be grouped according to the values of one (or more) phenophase(s) and separate seasonal trajectories can be obtained (fig. 8). By265 applying this procedure to the torgnon-ND grassland site we were able to iden- tify a bimodal distribution in maturity phase (i.e. the end of spring growth), as extracted after Klosterman curve fitting. We then selected pixels with maturity onset dates falling around the two modes (i.e. DOY 180±3 and DOY 205±3) 12 and computed an average trajectory for these two groups. This was possible270 because in the pixel-based approach each pixel has associated curve parameters to reconstruct the seasonal trajectory. The resulting average trajectories are markedly distinct around the seasonal peak, and reflect the different ecology of forbs and grasses occurring over very short distances (<20 cm) at this grassland site (Julitta et al., 2014). Further interpretation of this behavior is beyond the275 objective of this paper, but this preliminary analysis highlights the potential use of spatially explicit image processing to identify different phenology related to different plant functional types. More sophisticated clustering methods can be applied to pixel values in order to distinguish different phenological patterns. For example, at the Bartlett site (a280 mixed deciduous forest) we might expect a different phenology in early or late flushing species (fig. 9a). With a cluster analysis (Hartigan and Wong, 1979) using phenophases extracted with all available methods as an input matrix, we were able to clearly distinguish two portions of the canopy (fig. 9b). In this example the cluster analysis was used to define two sub-ROIs, that reentered the285 processing chain (ROI-averaged approach, fig. 2) and led to the identification of two markedly different seasonal trajectories (fig. 9c,d). Beech-dominated portions of the canopy (cluster 2) show a later yellowing/browning compared to maple, prevailing in cluster 1 (Richardson et al., 2009). 290 6. Summary and outlook We presented a new R package called phenopix that collects the most up- to-date processing techniques for repeated digital photography of vegetation canopy in the context of phenological and ecological research. Together with the basic and well established processing steps, the package implements several new295 features, including the possibility to combine different fitting and phenophase methods, to compute uncertainty on extracted phenophases, and to conduct pixel-by-pixel analysis. Phenophase maps obtained by the pixel-based approach 13 may represent a promising tool for a variety of ecological analyses, including for example: 1) the study of species-specific phenology; 2) the examination of300 spatial variability of apparently homogeneous ecosystems such as grasslands; 3) the potential to discriminate between overstory and understory or age-related differences in phenological patterns of open canopy forests. One major future step in the evolution of phenopix package includes the com- putation of the so called camera-NDVI, in the formulation proposed by Petach305 et al. (2014). The phenopix R package offers a suite of standardized and reproducible process- ing code that may be adopted, deployed and further developed by the existing phenocameras networks worldwide. Acknowledgments310 ADR acknowledges support from the National Science Foundation, through the Macrosystems Biology program (grant EF-1065029). The package was developed in the framework of the ALCOTRA Interreg Project n 227 e-PHENO Reti fenologiche nelle Alpi. We thank Morgan Furze, Harvard University, for assistance with editing the315 manuscript. References Ahrends, H.E., Brügger, R., Stöckli, R., Schenk, J., Michna, P., Jeanneret, F., Wanner, H., Eugster, W., 2008. Quantitative phenological observations of a mixed beech forest in northern Switzerland with digital photography.320 Journal of Geophysical Research 113, G04004. URL: http://doi.wiley. com/10.1029/2007JG000650, doi:10.1029/2007JG000650. Beck, P.S.A., Atzberger, C., Hø gda, K.A., Johansen, B., Skidmore, A.K., 2006. Improved monitoring of vegetation dynamics at very high latitudes: A new method using {MODIS} {NDVI}. Remote Sens-325 ing of Environment 100, 321–334. URL: http://www.sciencedirect. 14 http://doi.wiley.com/10.1029/2007JG000650 http://doi.wiley.com/10.1029/2007JG000650 http://doi.wiley.com/10.1029/2007JG000650 http://dx.doi.org/10.1029/2007JG000650 http://www.sciencedirect.com/science/article/pii/S0034425705003640 http://www.sciencedirect.com/science/article/pii/S0034425705003640 http://www.sciencedirect.com/science/article/pii/S0034425705003640 com/science/article/pii/S0034425705003640, doi:http://dx.doi.org/ 10.1016/j.rse.2005.10.021. Busetto, L., Colombo, R., Migliavacca, M., Cremonese, E., Meroni, M., Galvagno, M., Rossini, M., Siniscalco, C., Morra Di Cella, U., Pari,330 E., 2010. Remote sensing of larch phenological cycle and analysis of relationships with climate in the Alpine region. Global Change Bi- ology URL: http://doi.wiley.com/10.1111/j.1365-2486.2010.02189.x, doi:10.1111/j.1365-2486.2010.02189.x. Delbart, N., Kergoat, L., Toan, T.L., Lhermitte, J., Picard, G.,335 2005. Determination of phenological dates in boreal regions using normalized difference water index. Remote Sensing of Environment 97, 26–38. URL: http://www.sciencedirect.com/science/article/ pii/S0034425705001288, doi:http://dx.doi.org/10.1016/j.rse.2005. 03.011.340 Elmore, A.J., Guinn, S.M., Minsley, B.J., Richardson, A.D., 2012. Landscape controls on the timing of spring, autumn, and growing season length in mid-Atlantic forests. Global Change Biology 18, 656–674. doi:10.1111/j. 1365-2486.2011.02521.x. Filippa, G., Cremonese, E., Galvagno, M., Migliavacca, M., Morra di345 Cella, U., Petey, M., Siniscalco, C., 2015. Five years of phenological monitoring in a mountain grassland: inter-annual patterns and evalua- tion of the sampling protocol. International Journal of Biometeorology URL: http://link.springer.com/10.1007/s00484-015-0999-5, doi:10. 1007/s00484-015-0999-5.350 Forkel, M., Migliavacca, M., Thonicke, K., Reichstein, M., Schaphoff, S., Weber, U., Carvalhais, N., 2015. Codominant water control on global interannual variability and trends in land surface phenology and greenness. Global Change Biology 21, 3414–3435. doi:10.1111/gcb.12950. 15 http://www.sciencedirect.com/science/article/pii/S0034425705003640 http://www.sciencedirect.com/science/article/pii/S0034425705003640 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2005.10.021 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2005.10.021 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2005.10.021 http://doi.wiley.com/10.1111/j.1365-2486.2010.02189.x http://dx.doi.org/10.1111/j.1365-2486.2010.02189.x http://www.sciencedirect.com/science/article/pii/S0034425705001288 http://www.sciencedirect.com/science/article/pii/S0034425705001288 http://www.sciencedirect.com/science/article/pii/S0034425705001288 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2005.03.011 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2005.03.011 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2005.03.011 http://dx.doi.org/10.1111/j.1365-2486.2011.02521.x http://dx.doi.org/10.1111/j.1365-2486.2011.02521.x http://dx.doi.org/10.1111/j.1365-2486.2011.02521.x http://link.springer.com/10.1007/s00484-015-0999-5 http://dx.doi.org/10.1007/s00484-015-0999-5 http://dx.doi.org/10.1007/s00484-015-0999-5 http://dx.doi.org/10.1007/s00484-015-0999-5 http://dx.doi.org/10.1111/gcb.12950 Galvagno, M., Wohlfahrt, G., Cremonese, E., Rossini, M., Colombo, R., Filippa,355 G., Julitta, T., Manca, G., Siniscalco, C., Morra di Cella, U., Migliavacca, M., 2013. Phenology and carbon dioxide source/sink strength of a subalpine grassland in response to an exceptionally short snow season. Environmental Research Letters 8, 025008. Gillespie, A.R., Kahle, A.B., Walker, R.E., 1987. Color enhancement of highly360 correlated images. II. Channel ratio and chromaticity transformation tech- niques. Remote Sensing of Environment 22, 343–365. Gu, L., Post, W., Baldocchi, D., Black, T., Suyker, A., Verma, S., Vesala, T., Wofsy, S., 2009. Characterizing the seasonal dynamics of plant com- munity photosynthesis across a range of vegetation types, in: Noormets, A.365 (Ed.), Phenology of Ecosystem Processes. Springer New York, pp. 35–58. URL: http://dx.doi.org/10.1007/978-1-4419-0026-5_2, doi:10.1007/ 978-1-4419-0026-5_2. Hartigan, J.A., Wong, M.A., 1979. Algorithm as 136: A k-means clustering al- gorithm. Journal of the Royal Statistical Society. Series C (Applied Statistics)370 28, pp. 100–108. URL: http://www.jstor.org/stable/2346830. Hufkens, K., Friedl, M., Sonnentag, O., Braswell, B.H., Milliman, T., Richard- son, A.D., 2012. Linking near-surface and satellite remote sensing measure- ments of deciduous broadleaf forest phenology. Remote Sensing of Environ- ment 117, 307–321. URL: http://linkinghub.elsevier.com/retrieve/375 pii/S0034425711003543, doi:10.1016/j.rse.2011.10.006. Ide, R., Oguma, H., 2013. A cost-effective monitoring method using dig- ital time-lapse cameras for detecting temporal and spatial variations of snowmelt and vegetation phenology in alpine ecosystems. Ecological Informat- ics 16, 25–34. URL: http://www.sciencedirect.com/science/article/380 pii/S1574954113000290, doi:10.1016/j.ecoinf.2013.04.003. Jenkins, J.P., Richardson, A.D., Braswell, B.H., Ollinger, S.V., Hollinger, D.Y., Smith, M.L., 2007. Refining light-use efficiency calculations for 16 http://dx.doi.org/10.1007/978-1-4419-0026-5_2 http://dx.doi.org/10.1007/978-1-4419-0026-5_2 http://dx.doi.org/10.1007/978-1-4419-0026-5_2 http://dx.doi.org/10.1007/978-1-4419-0026-5_2 http://www.jstor.org/stable/2346830 http://linkinghub.elsevier.com/retrieve/pii/S0034425711003543 http://linkinghub.elsevier.com/retrieve/pii/S0034425711003543 http://linkinghub.elsevier.com/retrieve/pii/S0034425711003543 http://dx.doi.org/10.1016/j.rse.2011.10.006 http://www.sciencedirect.com/science/article/pii/S1574954113000290 http://www.sciencedirect.com/science/article/pii/S1574954113000290 http://www.sciencedirect.com/science/article/pii/S1574954113000290 http://dx.doi.org/10.1016/j.ecoinf.2013.04.003 a deciduous forest canopy using simultaneous tower-based carbon flux and radiometric measurements. Agricultural and Forest Meteorology385 143, 64–79. URL: http://www.sciencedirect.com/science/article/ pii/S0168192306003613, doi:http://dx.doi.org/10.1016/j.agrformet. 2006.11.008. Julitta, T., Cremonese, E., Migliavacca, M., Colombo, R., Galvagno, M., Siniscalco, C., Rossini, M., Fava, F., Cogliati, S., di Cella, U.M., Menzel,390 A., 2014. Using digital camera images to analyse snowmelt and phenol- ogy of a subalpine grassland. Agricultural and Forest Meteorology 198199, 116 – 125. URL: http://www.sciencedirect.com/science/article/ pii/S0168192314001981, doi:http://dx.doi.org/10.1016/j.agrformet. 2014.08.007.395 Keenan, T.F., Darby, B., Felts, E., Sonnentag, O., Friedl, M., Hufkens, K., O’ Keefe, J.F., Klosterman, S., Munger, J.W., Toomey, M., Richardson, A.D., 2014. Tracking forest phenology and seasonal physiology using dig- ital repeat photography: a critical assessment. Ecological Applications , 140206175103002URL: http://www.esajournals.org/doi/abs/10.1890/400 13-0652.1, doi:10.1890/13-0652.1. Kline, M., 1998. Calculus: An Intuitive and Physical Approach (Second Edi- tion)(Google eBook). Dover Publications. URL: http://books.google.com/ books?id=tZy8AQAAQBAJ&pgis=1. Klosterman, S.T., Hufkens, K., Gray, J.M., Melaas, E., Sonnentag, O., Lavine,405 I., Mitchell, L., Norman, R., Friedl, M.A., Richardson, A.D., 2014. Evaluating remote sensing of deciduous forest phenology at multiple spatial scales using PhenoCam imagery. Biogeosciences 11, 4305–4320. Lechowicz, M.J., 1984. Why do temperate deciduous trees leaf out at different times? adaptation and ecology of forest communities. The American Natu-410 ralist 124, pp. 821–842. URL: http://www.jstor.org/stable/2461303. 17 http://www.sciencedirect.com/science/article/pii/S0168192306003613 http://www.sciencedirect.com/science/article/pii/S0168192306003613 http://www.sciencedirect.com/science/article/pii/S0168192306003613 http://dx.doi.org/http://dx.doi.org/10.1016/j.agrformet.2006.11.008 http://dx.doi.org/http://dx.doi.org/10.1016/j.agrformet.2006.11.008 http://dx.doi.org/http://dx.doi.org/10.1016/j.agrformet.2006.11.008 http://www.sciencedirect.com/science/article/pii/S0168192314001981 http://www.sciencedirect.com/science/article/pii/S0168192314001981 http://www.sciencedirect.com/science/article/pii/S0168192314001981 http://dx.doi.org/http://dx.doi.org/10.1016/j.agrformet.2014.08.007 http://dx.doi.org/http://dx.doi.org/10.1016/j.agrformet.2014.08.007 http://dx.doi.org/http://dx.doi.org/10.1016/j.agrformet.2014.08.007 http://www.esajournals.org/doi/abs/10.1890/13-0652.1 http://www.esajournals.org/doi/abs/10.1890/13-0652.1 http://www.esajournals.org/doi/abs/10.1890/13-0652.1 http://dx.doi.org/10.1890/13-0652.1 http://books.google.com/books?id=tZy8AQAAQBAJ&pgis=1 http://books.google.com/books?id=tZy8AQAAQBAJ&pgis=1 http://books.google.com/books?id=tZy8AQAAQBAJ&pgis=1 http://www.jstor.org/stable/2461303 Migliavacca, M., Cremonese, E., Colombo, R., Busetto, L., Galvagno, M., Ganis, L., Meroni, M., Pari, E., Rossini, M., Siniscalco, C., Morra di Cella, U., 2008. European larch phenology in the Alps: can we grasp the role of ecological factors by combining field observations and inverse modelling? International415 journal of biometeorology 52, 587–605. URL: http://www.ncbi.nlm.nih. gov/pubmed/18437430, doi:10.1007/s00484-008-0152-9. Migliavacca, M., Galvagno, M., Cremonese, E., Rossini, M., Meroni, M., Son- nentag, O., Cogliati, S., Manca, G., Diotri, F., Busetto, L., Cescatti, A., Colombo, R., Fava, F., Morra di Cella, U., Pari, E., Siniscalco, C., Richard-420 son, A.D., 2011. Using digital repeat photography and eddy covariance data to model grassland phenology and photosynthetic CO2 uptake. Agricultural and Forest Meteorology 151, 1325–1337. Mizunuma, T., Wilkinson, M., L. Eaton, E., Mencuccini, M., I. L. Morison, J., Grace, J., 2013. The relationship between carbon dioxide uptake and canopy425 colour from two camera systems in a deciduous forest in southern England. Functional Ecology 27, 196–207. URL: http://doi.wiley.com/10.1111/ 1365-2435.12026, doi:10.1111/1365-2435.12026. Nasahara, K., Nagai, S., 2015. Review: Development of an in-situ observation network for terrestrial ecological remote sensing - the phenological eyes net-430 work (pen). Ecological Research 30, 211–223. URL: http://dx.doi.org/ 10.1007/s11284-014-1239-x, doi:10.1007/s11284-014-1239-x. Papale, D., Reichstein, M., Aubinet, M., Canfora, E., Bernhofer, C., Kutsch, W., Longdoz, B., Rambal, S., Valentini, R., Vesala, T., Yakir, D., 2006. To- wards a standardized processing of Net Ecosystem Exchange measured with435 eddy covariance technique: algorithms and uncertainty estimation. Biogeo- sciences 3, 571–583. Petach, A.R., Toomey, M., Aubrecht, D.M., Richardson, A.D., 2014. Mon- itoring vegetation phenology using an infrared-enabled security camera. Agricultural and Forest Meteorology 195-196, 143–151. URL: http:440 18 http://www.ncbi.nlm.nih.gov/pubmed/18437430 http://www.ncbi.nlm.nih.gov/pubmed/18437430 http://www.ncbi.nlm.nih.gov/pubmed/18437430 http://dx.doi.org/10.1007/s00484-008-0152-9 http://doi.wiley.com/10.1111/1365-2435.12026 http://doi.wiley.com/10.1111/1365-2435.12026 http://doi.wiley.com/10.1111/1365-2435.12026 http://dx.doi.org/10.1111/1365-2435.12026 http://dx.doi.org/10.1007/s11284-014-1239-x http://dx.doi.org/10.1007/s11284-014-1239-x http://dx.doi.org/10.1007/s11284-014-1239-x http://dx.doi.org/10.1007/s11284-014-1239-x http://linkinghub.elsevier.com/retrieve/pii/S0168192314001257 http://linkinghub.elsevier.com/retrieve/pii/S0168192314001257 http://linkinghub.elsevier.com/retrieve/pii/S0168192314001257 //linkinghub.elsevier.com/retrieve/pii/S0168192314001257, doi:10. 1016/j.agrformet.2014.05.008. R Core Team, 2015. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. URL: http://www. R-project.org/.445 Richardson, A., Jenkins, J., Braswell, B., Hollinger, D., Ollinger, S., Smith, M.L., 2007. Use of digital webcam images to track spring green-up in a deciduous broadleaf forest. Oecologia 152, 323–334. URL: http://dx.doi. org/10.1007/s00442-006-0657-z, doi:10.1007/s00442-006-0657-z. Richardson, A.D., Bailey, A.S., Denny, E.G., Martin, C.W., O’Keefe, J., 2006.450 Phenology of a northern hardwood forest canopy. Global Change Biology 12, 1174–1188. Richardson, A.D., Braswell, B.H., Hollinger, D.Y., Jenkins, J.P., Ollinger, S.V., 2009. Near-surface remote sensing of spatial and temporal variation in canopy phenology. Ecological Applications 19, 1417–28.455 Richardson, A.D., Keenan, T.F., Migliavacca, M., Ryu, Y., Sonnentag, O., Toomey, M., 2013. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agricultural and Forest Meteo- rology 169, 156–173. Sonnentag, O., Hufkens, K., Teshera-Sterne, C., Young, A.M., Friedl, M.,460 Braswell, B.H., Milliman, T., OKeefe, J., Richardson, A.D., 2012. Digital repeat photography for phenological research in forest ecosystems. Agricul- tural and Forest Meteorology 152, 159–177. Soudani, K., Hmimina, G., Delpierre, N., Pontailler, J.Y., Aubinet, M., Bonal, D., Caquet, B., de Grandcourt, a., Burban, B., Flechard, C., Guyon, D.,465 Granier, a., Gross, P., Heinesh, B., Longdoz, B., Loustau, D., Moureaux, C., Ourcival, J.M., Rambal, S., Saint André, L., Dufrêne, E., 2012. Ground- based Network of NDVI measurements for tracking temporal dynamics of 19 http://linkinghub.elsevier.com/retrieve/pii/S0168192314001257 http://linkinghub.elsevier.com/retrieve/pii/S0168192314001257 http://dx.doi.org/10.1016/j.agrformet.2014.05.008 http://dx.doi.org/10.1016/j.agrformet.2014.05.008 http://dx.doi.org/10.1016/j.agrformet.2014.05.008 http://www.R-project.org/ http://www.R-project.org/ http://www.R-project.org/ http://dx.doi.org/10.1007/s00442-006-0657-z http://dx.doi.org/10.1007/s00442-006-0657-z http://dx.doi.org/10.1007/s00442-006-0657-z http://dx.doi.org/10.1007/s00442-006-0657-z canopy structure and vegetation phenology in different biomes. Remote Sens- ing of Environment 123, 234–245. URL: http://dx.doi.org/10.1016/j.470 rse.2012.03.012, doi:10.1016/j.rse.2012.03.012. White, M.A., Nemani, R.R., 2006. Real-time monitoring and short-term forecasting of land surface phenology. Remote Sensing of Environment 104, 43–49. URL: http://www.sciencedirect.com/science/article/ pii/S0034425706001660, doi:http://dx.doi.org/10.1016/j.rse.2006.475 04.014. Wingate, L., Ogée, J., Cremonese, E., Filippa, G., Mizunuma, T., Migliavacca, M., Moisy, C., Wilkinson, M., Moureaux, C., Wohlfahrt, G., Hammerle, A., Hörtnagl, L., Gimeno, C., Porcar-Castell, A., Galvagno, M., Nakaji, T., Mori- son, J., Kolle, O., Knohl, A., Kutsch, W., Kolari, P., Nikinmaa, E., Ibrom,480 A., Gielen, B., Eugster, W., Balzarolo, M., Papale, D., Klumpp, K., Köstner, B., Grünwald, T., Joffre, R., Ourcival, J.M., Hellstrom, M., Lindroth, A., Charles, G., Longdoz, B., Genty, B., Levula, J., Heinesch, B., Sprintsin, M., Yakir, D., Manise, T., Guyon, D., Ahrends, H., Plaza-Aguilar, A., Guan, J.H., Grace, J., 2015. Interpreting canopy development and physiology us-485 ing the EUROPhen camera network at flux sites. Biogeosciences Discussions 12, 7979–8034. URL: http://www.biogeosciences-discuss.net/12/7979/ 2015/, doi:10.5194/bgd-12-7979-2015. Woebbecke, D., Meyer, G., Von Bargen, K., Mortensen, D., 1995. Color in- dices for weed identification under various soil, residue, and lighting condi-490 tions. Transactions of the ASAE 38, 259–269. URL: http://cat.inist.fr/ ?aModele=afficheN&cpsidt=3503524. Zhang, X., Friedl, M.a., Schaaf, C.B., Strahler, A.H., Hodges, J.C.F., Gao, F., Reed, B.C., Huete, A., 2003. Monitoring vegetation phenology using MODIS. Remote Sensing of Environment 84, 471–475. doi:10.1016/S0034-4257(02)495 00135-9. 20 http://dx.doi.org/10.1016/j.rse.2012.03.012 http://dx.doi.org/10.1016/j.rse.2012.03.012 http://dx.doi.org/10.1016/j.rse.2012.03.012 http://dx.doi.org/10.1016/j.rse.2012.03.012 http://www.sciencedirect.com/science/article/pii/S0034425706001660 http://www.sciencedirect.com/science/article/pii/S0034425706001660 http://www.sciencedirect.com/science/article/pii/S0034425706001660 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2006.04.014 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2006.04.014 http://dx.doi.org/http://dx.doi.org/10.1016/j.rse.2006.04.014 http://www.biogeosciences-discuss.net/12/7979/2015/ http://www.biogeosciences-discuss.net/12/7979/2015/ http://www.biogeosciences-discuss.net/12/7979/2015/ http://dx.doi.org/10.5194/bgd-12-7979-2015 http://cat.inist.fr/?aModele=afficheN&cpsidt=3503524 http://cat.inist.fr/?aModele=afficheN&cpsidt=3503524 http://cat.inist.fr/?aModele=afficheN&cpsidt=3503524 http://dx.doi.org/10.1016/S0034-4257(02)00135-9 http://dx.doi.org/10.1016/S0034-4257(02)00135-9 http://dx.doi.org/10.1016/S0034-4257(02)00135-9 Figure 1: A sample image from all PhenoCam sites included in this study masked on the chosen region of interest (ROI), along with the seasonal trajectories of filtered GCC in year 2013. Ecosystem type abbreviations are as in table 1. 21 Figure 2: The work-flow of the processing chain in the phenopix R package. Red and blue objects denote specific features of pixel-based and ROI-averaged approaches, respectively. Pixel-based approach produces a phenophase map (left), whereas ROI-averaged approach results in a GCC seasonal trajectory and a certain number of phenophases, depending on the method chosen (right). The example is from torgnon-ND grassland site in year 2013. 22 Figure 3: Illustration of the methods used to extract phenological thresholds (phenophases) in a seasonal GCC trajectory (f(t)). a) trs method: Phenophases are defined as the day of year (DOY) when the 50% GCC is reached either during greenup (sos) and autumn (eos), the DOY of maximum GCC is called pop; b) derivatives method: Phenophases are defined as the DOY when f’(t) shows the absolute maximum and minimum; c) Klosterman method: Phenophases are defined as the two local maxima (greenup) and two local minima (autumn) in the rate of change in curvature k’ (Kline, 1998; Klosterman et al., 2014) and are named after Zhang et al. (2003); d) Gu method: Maximum and minimum of f’(t) are used to define slopes of recovery and senescence lines tangent to the curves (green and brown lines). The intersection between these lines and baseline and maxline define the four phenophases in the original formulation (Gu et al., 2009). To account for the mid season decrease in GCC we have further defined a plateau line as a linear fit to GCC values between SD and DD, in order to adjust the definition of phenophase DD. For phenophase abbreviations and details see table 3. 23 Figure 4: All fitting and phenophase methods applied to Harvard deciduous forest GCC seasonal trajectory in year 2013. This plot is the output from the functions greenExplore() and plotExplore(). Plots in the same row share the same fitting method, whereas those in the same column share the same phenophase extraction method. The RMSE for each fitting is also annotated in the first column. 24 Figure 5: Phenophases as determined by the break point analysis (PhenoBP() function) applied to Kamuela grassland GCC seasonal trajectory in year 2013. The thickness of the vertical dashed lines is proportional to the uncertainty of the extracted break points (same scale as the x-axis). Uncertainty is estimated using the function confint() as the 95% confidence interval. 25 Figure 6: The uncertainty analysis applied to Harvard deciduous forest GCC seasonal trajec- tory in year 2013. The grey lines represent the 500 simulated curves on which phenophases are extracted. Shaded areas on phenophases represent 10th and 90th percentiles of the 500 replications. Data are fitted using the Klosterman method. Phenophases are extracted with a) the trs method; b) the derivatives method, c) klosterman method, d) gu method. 26 Figure 7: Scatterplot matrix showing: a) above the diagonal, the relationship between phenophases extracted with two different approaches (pixel-based and ROI-averaged) and two different fitting methods (spline smoothing and Klosterman fit), and b) below the diag- onal, errors associated. Error bars for pixel approach represent the mean absolute deviation of all pixels, whereas error bars for ROI-averaged approach represent the uncertainty as de- scribed in section 4. Note that the scale for error bars is provided in the first below-diagonal panel. Different symbols denote different phenophase extraction methods. Close symbols are for spring phenology and open symbols for autumn phenology. Different colors represent dif- ferent PFTs (abbreviations for PFTs are as in table 1). Data from all sites except for Kamuela were used in this analysis. 27 Figure 8: Results of the spatial analysis conducted on torgnon-ND grassland site. a) Map of the spatial distribution of maturity (Klosterman fitting method, Klosterman phenophases method), with darker colors indicating earlier occurrence of the phase (in DOYs). b) Aver- age seasonal trajectory of pixels clustered according to the bimodal distribution of maturity. Patches with forbs dominating show an earlier maturity compared to grasses (dashed lines). In the inset, density distribution of Maturity phase across all pixels. Shaded areas denote the maturity intervals used to subset ROI pixels. The resulting subsets were in turn used to extract the seasonal trajectories of forbs- and grass-dominated portions of the ROI. 28 Figure 9: Results of the cluster analysis run on Bartlett pixel-based phenophases. a) Sample image from DOY 275, b) Distribution of clusters in the ROI, c) Density distribution of Greenup and Dormancy phases across all pixels, and d) Averaged seasonal trajectories for the two clusters. All extracted phases (from the four methods presented in section 3) were used in input to clustering. 29