sensors
Article
Observing Spring and Fall Phenology in a Deciduous Forest with Aerial Drone Imagery
Stephen Klosterman 1,* ID and Andrew D. Richardson 1,2,3 1 Department of Organismic and Evolutionary Biology, Harvard University, Cambridge, MA 02138, USA; Andrew.Richardson@nau.edu 2 School of Informatics, Computing and Cyber Systems, Northern Arizona University, Flagstaff, AZ 86011, USA 3 Center for Ecosystem Science and Society, Northern Arizona University, Flagstaff, AZ 86011, USA * Correspondence: steve.klosterman@gmail.com
Received: 29 September 2017; Accepted: 5 December 2017; Published: 8 December 2017
Abstract: Plant phenology is a sensitive indicator of the effects of global change on terrestrial ecosystems and controls the timing of key ecosystem functions including photosynthesis and transpiration. Aerial drone imagery and photogrammetric techniques promise to advance the study of phenology by enabling the creation of distortion-free orthomosaics of plant canopies at the landscape scale, but with branch-level image resolution. The main goal of this study is to determine the leaf life cycle events corresponding to phenological metrics derived from automated analyses based on color indices calculated from drone imagery. For an oak-dominated, temperate deciduous forest in the northeastern USA, we find that plant area index (PAI) correlates with a canopy greenness index during spring green-up, and a canopy redness index during autumn senescence. Additionally, greenness and redness metrics are significantly correlated with the timing of budburst and leaf expansion on individual trees in spring. However, we note that the specific color index for individual trees must be carefully chosen if new foliage in spring appears red, rather than green—which we observed for some oak trees. In autumn, both decreasing greenness and increasing redness correlate with leaf senescence. Maximum redness indicates the beginning of leaf fall, and the progression of leaf fall correlates with decreasing redness. We also find that cooler air temperature microclimates near a forest edge bordering a wetland advance the onset of senescence. These results demonstrate the use of drones for characterizing the organismic-level variability of phenology in a forested landscape and advance our understanding of which phenophase transitions correspond to color-based metrics derived from digital image analysis.
Keywords: phenology; Harvard Forest; leaf color; plant area index; drone; UAV

1. Introduction
Phenology, the study of recurrent biological events, has been a focus of plant science for centuries [1]. Phenology responds to interannual and spatial variability in environmental conditions, particularly temperature [2,3], and also mediates key ecosystem functions, including carbon assimilation and evapotranspiration [4–8]. As global environmental change becomes an increasingly public issue, the value of plant phenology as a sensitive indicator of the effects of change has kindled interest in creating and interpreting phenology records, such as those from digital repeat photography. Tower- or building-mounted phenocams, typically located at positions just above the canopy [9–11], preserve a valuable visual record of vegetation phenology in forests and other ecosystems. “Near-surface” methods such as phenocams complement the phenology records of satellite remote sensing, which extensively observe entire landscapes, but at a spatial resolution that typically makes it impossible to discern individual plants [12]. Unlike other near-surface methods

Sensors 2017, 17, 2852; doi:10.3390/s17122852

www.mdpi.com/journal/sensors

Sensors 2017, 17, 2852

2 of 17

based on radiometric sensors, e.g., [13,14], images from phenocams have the advantage that they can be visually interpreted by direct examination.
With phenocam imagery, it is also possible to objectively quantify phenology from time series of vegetation “greenness” or “redness”. Discrete phenophase transition dates can be derived from those time series, using techniques developed for analyzing satellite remote sensing data [15,16]. These methods have been successfully applied to the evaluation of satellite remote sensing phenology products [17,18] and exploration of the connections between ecosystem function, e.g., carbon assimilation, and the greenness metrics of phenology [5,6]. However the oblique angle of tower-mounted cameras has been suggested to result in biased estimations of the timing of canopy maturity, compared to the vertically integrated plant area index (PAI), as well as in situ observations of leaf expansion [19]. Phenocam imagery does make it possible to easily integrate across multiple organisms in the camera field of view and thus characterize landscape-level phenology. But, in mixed-species stands, these estimates may be biased because the organisms located closest to the camera dominate the field of view and hence are over-represented, while more distant organisms are under-represented.
Drones, also called UAVs, open up an exciting new possibility for the near-surface remote sensing of plant phenology. With drones, researchers can obtain aerial photography with a similar spatial resolution to tower-mounted phenocams, but at the landscape scale, similar to satellite remote sensing. Compared to traditional aircraft, drones can be operated at a fraction of the cost, making more frequent observations feasible. Additionally, using photogrammetry techniques with drone images facilitates significant advances over tower-mounted cameras. Orthomosaics simulate an undistorted nadir perspective of the canopy, with a consistent spatial resolution over landscapes [20]. Because of this feature, orthomosaics enable the identification and analysis of larger numbers of individual organisms than is typically possible using tower-mounted camera imagery.
Previous research has begun to explore the potential of drones for the study of forest phenology. The first study to use aerial drone imagery to explore forest structural and color properties demonstrated that the changing color of leaves at different times of the growing season can be quantified by applying digital image processing techniques to georeferenced orthomosaics [20]. Subsequently, researchers leveraged the high spatial resolution of near-surface aerial drone photography to delineate crowns and study the phenology of individual trees. It was found that the color of individual tree crowns could be used to identify their species, based on knowledge of their expected phenological status during different times of the season [21]. In another study, high temporal frequency aerial images taken during springtime were used to show that the timing of budburst of individual trees, observed in situ, appears to coincide with the beginning of an increase in a canopy greenness metric calculated from the images of these trees [22]. In our previous work, which used an earlier year of imagery from the same study site reported here, we used drone imagery to show how landscape scale variance in phenology could be attributed to plant species, and explored the nature of the spatial variability of phenology at different spatial resolutions [23].
Our goal in this study is to extend previous efforts to determine the leaf life cycle events of trees that correspond to digital image analysis metrics of those same trees in aerial drone imagery. We use imagery of the complete growing season to go beyond budburst and consider leaf expansion in spring, as well as color change and abscission in autumn. We also compare metrics derived from aerial images with the progression of PAI, to interpret color indices with reference to canopy structural development. Finally, we examine how contrasting air temperature regimes at microsites within the study area correlate with spatial variation in landscape phenology, after accounting for spatial variance induced by differences in tree species. These results demonstrate the use of drones for observing the complete seasonal cycle of deciduous canopy development at the landscape scale, with a high enough spatial resolution to discern organism-level phenomena.

Sensors 2017, 17, 2852

3 of 17

Sensors 2017, 17, 2852
2. Materials and Methods

3 of 16

2. Materials and Methods 2.1. Study Site

W2.1e. ScotunddyuScitteed our study at Harvard Forest in Petersham, MA, USA (42.5377◦ N, 72.1715◦ W).
The climaWteeiscotenmdupcetreadteo,uwr istthudaymateaHnaravnanrudaFloprreestciipnitPaettieornshoafm1,1M0 cAm, UaSnAd (m42e.a5n37a7n°nNu,a7l2t.e1m71p5°erWat)u. re of 7.1 ◦CT.hTehcelimstautde yisateremapiesraatem, iwxietdh ademceidanuoanuns-ueavleprrgerceipeintaftoiornesot,f w11i0thcmsoamnde mweoaondaynwnueatllatenmdsp.eDraetucirdeuous trees oinf 7th.1e °sCtu. dTyhearsetuadiyncaluredaeispraedmoimxeidnadnetcliyduroeduso-eavker(gQreueenrcufosrersutb, rwa)itahndsormeed wmoaopdlye (wAectelarnrdusb.rum), as weDrulelbcraiudsmuy)o,eualsslotwrweeelslbiainrsctyhheel(lsoBtuwedtubyliaracrheaal(lBeingecthulaulnadieaelnplesrgeihsda)o,nmiAenimnsiasen)r,tilAcyamrneedrbioceaaeknc(hbQeue(Fcerhacgu(uFsasrgugubrsraag)nradannifddoilrfioeald)ia,m),aaanpndldebb(Allaacccekkr oak (Querocauks(vQeuluertciunsa)v.eluTthinea)s. tTuhdeystaurdeyaairsea2.i8s 2h.8ahianisnizsiezeaannddccoonnttaaiinnssaapppprorxoixmimatealtyel1y90109t0re0estrweeitshwa ith a diamdetiaemr aettebrraetabsrtehasetighheitg≥ht1≥010cmcm(F(Figiguurree11))..

FigurFeig1u. reSt1u.dSytuadryeaaraetaHatarHvaarrvdarFdoFreosretsot non5/52/211//1177 ((DDOOYY11414)1.)L. oLcoatciaotnioonf moficmroiscirteostietme pteemratpuerreature loggelrosgginedrsicinadteicdataesdbalsubelu“eד”×”sysmymbboolsls..
2.2. Drone Image Acquisition and Processing 2.2. Drone Image Acquisition and Processing
We used methods described in an earlier study [23] to obtain and process aerial photography.
WBreiefulyse, dwme eutshedodas ddreosncerib(3eDdRinAarndueCaorpliteerr sQtuudady-C[23F]ratomoe,bt3aDinRaonbdotipcrso, cBeesrskealeeryi,alCpAh, oUtoSgAr)aphy. Briefleyq, uwipepuesdewdiathdarConaneo(n3DPoRwAerrsdhuoCt Aop33te0r0 Qcaumaedr-aC(3F5rmammef,il3mDeqRuoibvoalteicnst,fBocearlkleelnegyt,hC2A8 m, UmS,Aap) perqouxi.pped with a16CmanilloionnPpoiwxeelsr)s.hWoteAus3e3d00thceasmamereac(a3m5emramanfidlmimeaqgueisveattlienngtsffoocraallllefnligghthts2, a8nmdmto,oakpppicrtouxr.e1s 6ofmaillion pixelsg)r.aWy reeufesreednctehseqsuaamree(CcaomloerCrahaenckderimclaasgseics,eXtt-irnitge,sGforarnadllRflaipgihdtss,,ManI,dUtSoAok) bpeifcotruereeascohffaligghrta.yFlriegfhetrence squarfceroel(qoCur eocnhlocayrnCgweha,esdcerkopeuerngcdhlliaynsgesvoicen,rytXh-feirviwteeed,aatGhyersardncuodrniRdngiatipsopindrsisn(,dgMalteeIa,sfUsohSuoAtwa)nnbdineefFoviergreuyreewa2ce)he.kIflmdigaughrietns.gwFaeluirgteuhtmatknferneleqaauft ency was raoumgihnilmy uevmersyhufitvteer dspaeyesddoufri1n/1g0s0p0 rsi,nwg iltehafcoonusttaanntdexepvoesruyrewdeuekrindgureiancgh afluigtuhtm. Tnhleeasfamcoelocrolcohrange, depenbadlianngceownatshueswedefaotrhaelrl caocqnudisititiioonnsd(adtaest,esbeschaouwsencoinnsFisitgeuntreco2lo).r Ibmalaagneces iws enreecetsaskareynfaotr aremliaibnliemum shuttdeirgsitpael ecdamoefra1/o1b0se0r0vast,iownisthofcpohnesntaonlotgeyx[p1o1]s.uWreedcuonrdinugcteedacfhligflhitgshatt. mTihde-dsaaym(ebectwoleoernb1a0laan.mce. was used afnodr a3lpl .amcq.)uoinsietiitohnerdcaletaers,obr eecvaeunlsyeocvoenrscaissttednatycso, alonrdbnaelvaenrcdeuirsinngecpeesrsioadrys offovrarreilaibalbelceloduigdictoavl ecra,mera
obseravsaetxiopnossuoref pwhaes ncoonlostgaynt[1d1u]r.inWg eflicgohntsd. uFoctreedacflhigimhatsgeartymaciqdu-disaityion(bdeatwte,ewene c1r0eaat.emd.oartnhdop3hopt.oms .) on eitherofcltehaersoturdeyvaernelayuosvinegrcaabsotudta1y0s0, aJPnEdGnpehvoetrosdutarkienng wpietrhioadnsinotfevrvaarlioambleetecrloscurdipct o(Cvaenr,oansHexacpkosure was cDoenvsetlaonptmdeunrtinKgitfl, ighhtttps.:/F/cohrdeka.wchikiiam.caogme/rwyiakci/qCuHisDitKio),n wdaitthe, wthee crPehaotteodScoarnthopphohtootgorsamofmtehterystudy area usosfitnwgaraeb(oAugtis1o0f0t, SJPt.EPGeteprhsboutrogs, tRauksesnia)w, ainthd apnerfinortmerevdaflionmalegteeorrsecferripentc(iCnganinoEnRHDaAcSkIDMeAvGeIlNopEment Kit, hAHttuaprt:ov/Sa/yrdcnhcFdo(Irkne.tswetriDgkraiatapa.hcAo, rHmchu/invwtesiv[k2ili4l/e]C., AHLD, UKS)A, w). iTthhetohrethPohpohtootSoscaunsepdhinottohgisrastmudmyeatrreyasvoafiltawblaerein(tAhegisoft, St. Petersburg, Russia), and performed final georeferencing in ERDAS IMAGINE AutoSync (Intergraph,
Huntsville, AL, USA). The orthophotos used in this study are available in the Harvard Forest Data
Archive [24].

Sensors 2017, 17, 2852
Sensors 2017, 17, 2852

4 of 17
4 of 16

Figure 2F.ig(uAre–L2.):(AC–lLo)s: eC-luopse-(u3p0 (m30 bmyb3y03m0 m) )ofofaaeerriiaall iimmaaggees sfofroar saelseeclteiocntioofndoatfesdfartoems f2r0o1m5. E2a0c1h5. Each close-upclsohseo-wupssthhoewRsOthIesRuOsIesdusteodatnoaalnyazlyeztewtwoorreeddooaakkttrreeees.sS. oSloidlimdamgeangtaenRtOaIsRwOeIrsewuseerdetuo scaeldcutloatcealculate
solid magenta GCC and RCC values (circles symbols) and curves; similar with the dashed cyan ROIs
solid ma(tgrieanngtaleGsyCmCbaonlsd.).RImCCagvesalhuaevse (bceiercnleeqsusaylimzebdotlos)haanvedtchue rsvaemse; sbirmighiltanresws i(tshumthoefdRa, sGh,eadndcyBan ROIs (triangledisgyitmal nbuomlsb.)e.rsI)mfoargveissuahliazvaetiobnepeunrpeoqsuesailniztheids ftioguhrea;v(Me ,tNh)e: tshaemreesubltrinigghGtnCCeasnsd(RsuCCmtimoef sRer,ieGs, and B

digital n(suymmbboelrs)s;)cuforvrevfiitssuoarliinztaertpioonlatpiounrsp; aonsdesesitnimtahteisd pfihgeunroep;h(aMse,tNra)n:sitthioenrdeastueslt(ivnegrtiGcaCl lCinaesn).dThReCC time

series (svyemrtbicoalsl)in; ecsusrhvoewfistsixoarninnutaelrdpaotleastfioornbso; tahnGdCeCsttimimeasteerdiesp, hwehniloepshparisnegtdratnesictiaolcnuldaatetdesfr(ovmerRtiCcCal lines).

The vertsiecraielslainreeosnslhyoshwowsinxfaorntnhue atrlede awtieths froedr sbportinhgGleCaCvetsi.mRCeC-sEeOriFeso,ccwurhrielde osnptrhine gsadmaetdesatceafolcrubloatthed from RCC seritreeseas.re only shown for the tree with red spring leaves. RCC-EOF occurred on the same date for

both treWese. calculated green chromatic coordinate (GCC) [25] and red chromatic coordinate (RCC) [11]

time series from orthophotos using Matlab (R2017a), to account for changes in scene illumination due

Wetocadlifcfuerleantceeds ignraetemnocsphhreormic acotincdcitoioonrsdoirntaimtees(GofCdCa)y[b2e5t]waenendflriegdhtsc:hromatic coordinate (RCC) [11]

time series from orthophotos using Matlab (GRC2C0=1G7a/()R, t+oGac+cBo)u, nt for changes in scene illumin(1a)tion due

to differences in atmospheric conditions or times of day between flights:
RCC = R/(R + G + B),

(2)

where R, G, and B are the mean redG, gCrCee=n,Gan/d(Rbl+ueGd+igiBta)l, numbers, respectively, in regions of
interest from image files. Example GCC and RCC time series and additional processing details are

(1)

available in Figure S1. We note that thReCpChe=noRlo/g(yRs+ignGal+obBs)e,rved in the GCC and RCC time series of

(2)

vegetation was approximately 10 times greater in amplitude than the noise observed in the analysis

where Ro,f tGhe, garnady rBefearreenctehsequmareea, nwhriecdh,agrorseeednu,eatnodfacbtloursesducihgiatsavl anryuimngbiellrusm, irneastpioencctoivnedliyti,onins arnedgions of

interestthfreodmiffeimrenatgteimfielseosf. dEayxaomf flpiglehtGs. CWCe alnsdo eRxpClCorteidmoethseerriniedsicaens,dinacldudiintigoGnCaCl+pRrCoCc, eGsRsVinI,gEdxGet,ails are available in Figure S1. We note that the phenology signal observed in the GCC and RCC time series of vegetation was approximately 10 times greater in amplitude than the noise observed in the analysis of

the gray reference square, which arose due to factors such as varying illumination conditions and the

different times of day of flights. We also explored other indices, including GCC + RCC, GRVI, ExG, and Hue, as described in the caption of Figure S2, but found the most reliable results from GCC and RCC, and used them in our analysis.

Sensors 2017, 17, 2852

5 of 17

Color indices were calculated using different regions of interest (ROIs) depending on the goal of the analysis. For comparison between image-based metrics and in situ observations of individual tree phenology, we drew ROIs around the crowns of the trees that we observed, as shown in Figure 2. For comparison to upward photo-based estimates of PAI, we created ROIs representing the canopy projection of the image field of view from the ground perspective. To examine phenology across the landscape, we used a square grid of 10 m ROIs (grid cells).
2.3. Estimating Phenology Dates from Time Series Data
We used curve fitting methods detailed in Klosterman et al. [18] to estimate phenology dates from time series data, including GCC and additional data described below. Sigmoid models are commonly used to approximate the seasonal trajectory of vegetation indices and to facilitate the determination of phenological transition dates. The greendown sigmoid is a variant of the standard dual sigmoid model [16]; a key difference is that the greendown sigmoid allows for a linear decrease in greenness over the course of the growing season. We used the greendown sigmoid model to approximate the seasonal trajectory of GCC for each ROI. Phenology dates were calculated from curve fit parameters by finding the dates of extrema in the curvature change rate (CCR) [15]. We used these methods to calculate dates for the start, middle, and end of canopy development in spring (GCC-SOS, GCC-MOS, GCC-EOS), and similar dates for canopy senescence in fall (GCC-SOF, GCC-MOF, GCC-EOF). We estimated uncertainties using the Jacobian matrix of curve fit parameters to generate Monte Carlo ensembles of phenology dates, and calculated the inner 95% confidence interval of these ensembles.
To estimate phenology dates from RCC time series, we calculated the dates of crossing the 10th, 50th, and 90th percentiles of a linear interpolation of RCC values, while RCC increased to its spring maximum, to determine RCC-SOS, RCC-MOS, and RCC-EOS. We used the date of the maximum redness value in fall to represent RCC-EOF. We determined uncertainties as the interval between the nearest two observation dates, or for transitions that occurred on an observation date, as the interval between the midpoints of that date and the nearest two dates.
In addition to these automated curve fitting and interpolation approaches, we also visually identified the day of the first observable leaves in aerial images of the individual trees discussed below in Section 2.4.1.
2.4. In Situ Measurements
We identified eight microsites within the study area based on earlier results, using imagery from a previous year at this site. We had previously found that species composition explained most of the spatial variance in phenology across the Harvard Forest landscape, using a multiple linear regression [23]. However, we identified locations within the study area where regression residuals were relatively large; these were the same areas in regressions using both the 2013 data from the previous study and the 2015 data from the present study. We located microsites in these places, as well as locations where residuals were relatively small, to better understand the drivers of phenology throughout the study area. To provide context for the image-derived phenological transition dates, we made direct observations of phenology on individual trees, and measured canopy-level PAI using digital cover photography. To characterize the microsite environment, we measured air temperature.
2.4.1. Direct Observation of Trees
We observed the phenological status of 30 trees spread across the eight microsites within the study area during each drone photography acquisition date (Figure 1, three to five trees per microsite), and one additional date when the drone was not flown due to weather constraints (5/17/15, DOY 137). Similar to a protocol of established observations of tree phenology at Harvard Forest [26], we estimated the percentage of buds burst on trees, and average leaf size (length) as a percentage of mature leaf size. In autumn, we estimated the percentage of leaves that had changed color due to senescence, and that had fallen. From these observations, we created time series using linear interpolation, and found the

Sensors 2017, 17, 2852

6 of 17

point in time corresponding to the deciles of phenological status for each phenology event: day of year for 10%, 20%, . . . 90% of budburst, leaf size, leaf color, and leaf fall. We calculated uncertainties in the same way used for RCC time series.
We determined the leaf life cycle events corresponding to image-derived metrics of individual tree ROIs by finding the decile of progress in a specific event that had the lowest RMSD with respect to a given metric across the 30 trees under observation. In other words, we compared all deciles (10%, 20%, . . . ) of all observed events (e.g., budburst, leaf size) to each of the image-derived transitions (SOS, MOS, . . . ), and identified the decile that had the minimum RMSD across all trees (e.g., 10% budburst was closest to SOS). We also examined correlations of GCC and RCC values on each date with percentages of progress in each life cycle event on the same date, within and across trees. Since color indices and leaf transitions typically trend in the same direction, for example, increasing GCC and leaf size in spring, an ordinary least squares regression may yield spuriously low coefficient standard errors. Therefore, we used an econometric approach to time series regression: heteroscedasticity and autocorrelation consistent (HAC) regression. HAC regression calculates the variances of regression coefficients (i.e., square of standard error, used to calculate regression p-value) based on the inferred autocorrelation structure of regression residuals; the coefficient estimates themselves are the same as those of ordinary least squares. We used the ‘hac’ function in Matlab, with pre-whitening at a lag of one time step and the quadratic spectral weighting scheme, as we found that these options led to the most conservative estimates of coefficient variances in a sensitivity analysis [27,28]. We used t-tests to see if there were significant differences in regression parameters between trees. We also calculated ordinary least squares regressions of pooled values across all trees.
2.4.2. Digital Cover Photography (DCP)
Upward photos were taken with the same model of camera used for drone photography (non-fisheye lens). We took images manually for upward photos, on the same dates we acquired drone photography, using the same image settings (i.e., automatic exposure), except for white balance, which was set to “auto” for DCP [29]. We used a level to aim the camera directly upward and positioned the camera on the posts holding the temperature data loggers (1 m height, below most mid-canopy trees), in the same orientation on each date. The projected camera field of view had an area of 533 m2 in the canopy (at the average Harvard Forest canopy height of 22 m). We estimated PAI from upward photos at each microsite using a freely available DCP software tool obtained from Craig MacFarlane (contact information in references [29,30]). We calculated sigmoid curve fits of PAI values to estimate transition dates as with GCC time series. However, we observed that PAI values were generally stable from June to September, so we used the traditional sigmoid model as opposed to the greendown sigmoid, effectively setting the slope of the summertime greendown to zero. We also directly compared GCC and RCC values with PAI values, by performing HAC regressions at the microsite level, as both PAI and color indices typically trend in the same direction at each microsite. We also performed ordinary least squares regressions of data pooled across microsites for spring and fall, similar to the analysis described in Section 2.4.1.
2.4.3. Air Temperature Measurements and Effects of Microclimate on Phenology
Because the microsite locations included places with the largest residuals in a regression of species composition and phenology dates, we suspected that additional factors may contribute to the timing of microsite phenology. Temperature effects on phenology have been widely studied [18,31–33]. Therefore, we recorded air temperatures half hourly from 11 April 2015, approximately one month before the beginning of leaf out, through December 2015, after completion of leaf fall. To do this, we installed a HOBO U23-004 temperature data logger at a 1 m height with an RS-1 radiation shield (Onset Computer Corp, Bourne, MA, USA) at each microsite.
To determine the effect of microclimate on phenology, we correlated spring and fall temperatures with residuals from a statistical model that accounted for the effects of species variation and land cover

Sensors 2017, 17, 2852

7 of 17

type, but not microclimate, on the aggregated canopy phenology of 10 m grid cells. A 10 m grid was chosen based on previous analysis, taking into account spatial inaccuracies in the orthophotos and the nature of the species data [23]. More specifically, the model used tree species from a map of all woody stems ≥1 cm diameter at breast height [34]. This map covered 89% of the land area monitored in drone imagery. Woody species composition was determined using all species that appeared in at least 10 of 245 grid cells, and had at least a 20% basal area fraction in any one grid cell. The remaining species were lumped into one predictor to ensure that the fractional species composition of each grid cell summed to 1. Any grid cell that had <50% coverage in the aerial images was eliminated from this analysis, removing partially covered grid cells on the edge of the imaged area and ensuring all data points had adequate spatial representation. We used a no-intercept regression since fractional species compositions are mixture components that sum to 1 [35]. We then calculated the average residual from these spring and fall regressions for the three by three windows of grid cells centered on the microsite locations (30 m by 30 m areas), and regressed these microsite-average residuals on the monthly means of daily minimum, mean, and maximum temperatures for April (MOS) and September (MOF), as temperatures preceding phenology events are commonly used to predict the timing of those events [18,31–33]. We note that we used RCC-MOS for grid cells with red spring leaves, following a criterion discussed in the Results Section 3.1, as long as the amplitude in RCC was greater than the range of noise observed in the gray reference square (Figure S1). For all other grid cells, we used GCC-MOS. We used Bonferroni correction to account for making three temperature comparisons with both phenology transition dates.
3. Results
3.1. Choice of Color Index in Spring Time
While changes in canopy greenness have commonly been used to track spring budburst and leaf expansion in deciduous forests, we found that in some instances, these processes appeared to be more associated with changes in redness (Figures S2 and S3). For two of the 30 trees we observed in situ, both red oaks, we saw that leaves were various shades of red (pink, orange) in color during leaf expansion (Figure 1, close up in Figure 2D,E, and Figure S4). Additional oak trees, which were not under in situ observation at the microsites, can also be seen to display red spring leaves in Figure 1. We note that leaves higher in the canopy, i.e., those on the top of the crown and most visible from an aerial view, often appeared to be redder than leaves closer to the ground (Figure S4).
For trees with red spring leaves, the springtime GCC profiles showed a delayed increase compared to trees with green spring leaves, including nearby conspecifics that we observed to have nearly the same leaf expansion phenology (Figure 2M,N). However, the springtime RCC time series exhibited a marked peak for trees with red spring leaves. Then, as RCC decreased, GCC increased as leaves became greener in late spring (Figure 2F,G). We found that for these trees, the springtime amplitude in RCC (the increase from the dormant value to the spring maximum) was more than 45% of the spring amplitude in GCC (46% and 64% for the two trees), while for all other trees, it was less than 35%. The springtime increase in RCC was closer in time to observed leaf expansion than GCC, as well as a range of other color indices we considered (GCC + RCC, GRVI, ExG, Hue, Figure S2). We also note that the leaves of oak trees with red spring leaves had no apparent color difference in autumn from those of trees with green spring leaves (Figure 2K).
3.2. Leaf Life Cycle Events of Trees: Correspondence to Image Metrics
We determined the leaf life cycle events corresponding to image metrics by finding the closest decile of life cycle event (i.e., 10% budburst) to each metric (SOS, MOS, etc.) for the 30 trees under observation. Following our observation that RCC was a better indicator of spring leaf phenology than GCC for trees with red spring leaves, we used RCC metrics for the two trees where the spring amplitude in RCC was greater than 40% of spring amplitude in GCC, and GCC metrics for all other trees.

Sensors 2017, 17, 2852

8 of 17

We found that SOS was closest to 10% budburst (RMSD 4.7 days), although the date of 10% budburst was not significantly correlated to SOS across trees (r = −0.04, p > 0.05). This low (and not statistically
significant) correlation is likely related to the fact that there was relatively little variability in the timing
of budburst across trees: all trees under observation had 10% budburst within a period of about a week
(DOY 121-129), while SOS date uncertainties were 11 days on average (inner 95% confidence interval, i.e., ±5.5 days). However, we found that our determinations of the first observable leaves in aerial
image orthomosaics were correlated with SOS (r = 0.46, p < 0.05), indicating SOS dates derived from
color-based vegetation indices were associated with biologically meaningful phenological transitions
at the branch-to-canopy level.
We found that MOS was closest to 40% leaf size (RMSD 3.6 days), and EOS was closest to 70% leaf
size (RMSD 6.9 days; spring time GCC, budburst, and leaf size observations with associated transition datesSsenhsoorws 2n017i,n17F, 2i8g5u2 re S5 for an example tree). Later spring transitions were more variab9loef 1a6cross
trees, with a range of 14 days for 40% leaf size and 18 days for 70% leaf size (Figure 3A,B). Across specielmesa,afxthcimoeluotwrmpoaenrmcdeonlsattategcroevsma(amluveoesrnawgseeprerec=cie0os.r9r,4er,leaadtlel odtraewkeisathpn<lde0ar.fe0fd1a,lmTl apabeplrelceeS,n3ht)a,agwdehssiil(metrfeiaellalarRvlCeeCraavfgaeelxurpe=sain−n0sc.il9ou4nd, i9pnhgofethn16eology. Howetrveeers, wAimthesruicffaicniebneteocbhseesrvwateiornesapm<o0n.0g1t)h. SeimfirilsatrttroeethsetospmriankgteimtheeaMnaOlySsisa,nwdeEfoOuSndtrsaingsniitfiiocannst, and attaincoarr4e0la%tioanndwh7e0n%poleoalifnsgizaell. tTrehees: mleaefacnolionrnweras9c5o%rrecloantefiddteonGceCCi,nwtehricvhaldseacrceraossesdtirneeasutwumerne(sri=x and 12 day−0s.8fo9,rpM<O0.S00a1n),dasEwOeSllmasetRrCicCsa,srietsipnecrcetaivseedlyt.oMthOeSfailsl mtyapxiicmaullmy t(hr e= m0.8o4s,tpce<r0t.a0i0n1)m. Leetarifcfainll wsiagsmoid curvemfiottrienhgigahnlaylycoserrse[la1t8e]ddtuoeRtCoC tahseitcdleeacretaimsedinagftoefr the fsatlelempeaxstimpuarmt o(rf =th−e0.i8n9c,rpea<s0e.0i0n1s)pthrianngGtCimC (er GCC (Figur=e−20.M56),.pC<o0n.0s0e1q)u, seinntcley,GwCCedfeocurenadsetdhadtutehteo McolOorScfhraonmgei,mbeafgoerealneavlyesibsecgoaunlfdalrlienpgr.eAsesnint tsrperein-tgo,-tree variatwioenniontetdhesigtinmifiincagnot fdilfefaefreenxcpeasninsiroengroebsssieornvseldopienssaintud iwntiethrcegpretsabteertwsteaetnisttrieceaslosfigthneifisacmanecseptehcaiens EOS (MOSinr t=he0s.e52co, mp =pa0r.i0so0n3s;.EOS r = 0.38, p = 0.04, Figure 3A,B).

GCC- or RCC-MOS, DOY

GCC-SOF, DOY

GCC- or RCC-EOS, DOY

40% Leaf size , RMSD = 3.6 days 145
A 140
135

70% Leaf size , RMSD = 6.9 days B
160
150

20% Leaf color , RMSD = 6.9 days C
300
280

130

125 125

r = 0.52
130 135 140 145 In situ DOY

40% Leaf color , RMSD = 6.4 days

320 D 310 300 290 280 270 260

r = 0.87

GCC-EOF, DOY

140

130 130

r = 0.38
140 150 160 In situ DOY

80% Leaf color , RMSD = 6.3 days 340
E

320

300

280 r = 0.90

RCC-EOF, DOY

260

240 240

r = 0.84
260 280 300 In situ DOY

10% Leaf fall , RMSD = 3.7 days 320
F
310

300

290

280 r = 0.92

260 280 300 320 In situ DOY

280 300 320 340 In situ DOY

280 290 300 310 320 In situ DOY

GCC-MOF, DOY

FigurFei3g.u(rAe –3F. )(Ale–aFf)lilefeafcylicfeleceyvcelentesvecnotrsrecsoprorensdpionngditnogdtroonderoimneaigmeamgeetmricestrficosr f3o0r t3r0eetsr.eeSsy.mSybmolbsodlsenote
transidtieonnotdeattreasn,switiiotnh dcoatnefis,dwenitchecionntfeidrveanlcsesihnotewrvnalassshcolwonreads lcionleosr.edYellilnoews. Y=erlelodwo=akre(dexocaekp(tetxrceeepst with red sptrreinesgwleitahvreesd, sshporiwngnlienavbelsa,cskhoinw(nAi,nBb))la, crkedin=(Are,Bd)m), raepdle=,rleigdhmt abplulee, l=igAhtmbelureic=anAmbeeercicha,ndbaerekchb,lue = yellowdabrkircbhlu,ean=dyeollroawngbeir=chb, laancdk ooraakn.gDe =ecbillaecskoofaikn. sDiteuciolebssoefrvinatsiiotunsob(is.eer.v, a4t0io%nsle(ai.fe.s,i4z0e%) wleearfessizeel)ected as havwinerge tsheleecmteidniams huamviRngMthSeDminincimomumpaRrMisoSDn tino cdormonpeariimsoangteo mdreotnreicism. aJigtetemr ewtraicssa. dJitdteerdwtoasdaadtdaepdoints that latoydoantathpeoisnatsmtheadt alatye oanndthtehseamoneed-aotneeanlidnethiesosnheo-wonne liinneblisacshko. wCnorirneblalatciokn. Ccoorerfefilactiieonntcsoaerfeficsihenotws n on each paraensehlo(walnl pon<e0a.c0h5)p.anel (all p < 0.05).
3.3. Linking Leaf Area to Color Indices To exWame dinereivheodwtrarnepsirtieosnendatatetsiv(SeOcSo,lMorOinS,deitcce.)sfraorme oPfAtIhteimpersoegrrieess,saisowneollfalseafrfolmifedrcoynceleims aatgethrye tree level,uwsiengcaRlOcuIslactoerdreHspAonCdrineggrteostshioentospb-eotfw-caeneonppyearrceeanatnaagleyszeodf lfiofrePcAycI l(eexcaommpplelemtiiocnroasintedsihnodwenx vinalues acrossFiogubrseer4v).atTiohins dalalotewse.dLuesaftosiezxeamobinseertvhaetiroenlastiwonerbeetmweoerne chaingohplyy cstorrurcetluartaelddwevietlhopGmCeCntvaanludes in springca(naovpeyracgoelorr =de0v.9el5o,p2m5eonft2i8n tdrreoensepim< a0g.0e1s., TThabelteimS1in)gthoafnthbeusdpbriunrgstinocbreseasrevaintioPnAsI (aanvderGaCgCewra=s 0.86, similar across the eight microsites, particularly for later spring transitions. At the first spring transition, SOS, the average difference between PAI and GCC derived dates was −4.1 (negative sign indicates GCC-MOS was earlier) ±3.2 days standard deviation, and was significantly different from zero according to a paired t-test (p < 0.05). However, the MOS and EOS transitions were closer in time, with −0.6 ± 2.5 days and 2.6 ± 4.6 days differences, respectively. Both of these differences were not significantly different from zero. In autumn, declines in GCC substantially preceded declines in

Sensors 2017, 17, 2852

9 of 17

21 of 28 trees p < 0.01). We also note that, pooling across all trees with green spring leaves, GCC values were highly correlated with leaf size observations (r = 0.88, p < 0.001), and to a similar but lesser extent, with budburst observations (r = 0.80, p < 0.001). However, we note that there were significant differences in regression slopes and intercepts between individual trees, even within species, according to t-tests (p < 0.05). This indicates that although there is a clear, statistically significant relationship between color indices and leaf life cycle progressions when all trees are pooled, there is also significant tree-to-tree variability in these relations when regression parameters are allowed to vary by tree.
For the two trees with red spring leaves, RCC increased, and then declined in spring as leaves continued to expand and change in color from red to green (Figure 2F,G). Therefore, we considered RCC values up to and including the spring maximum when examining correlations. Similar to GCC for trees with green spring leaves, we found a higher correlation of RCC with leaf size (Table S2, average r = 0.95, both trees p < 0.001) than budburst (average r = 0.76, neither tree p < 0.001). These results indicate that GCC and RCC values are more representative of the progression of leaf expansion than budburst, as these color indices are still increasing after the completion of budburst (e.g., Figures S2, S3 and S5).
The timing of leaf coloration events in autumn varies more among trees than either budburst or leaf expansion (Figure 3C–E) with a range of roughly 50 days in leaf color deciles across the 30 trees under observation. Species differences become more apparent in fall than in spring, with red oaks senescing after red maples, and American beeches and yellow birches in between these two. The SOF, MOF, and EOF transitions calculated from GCC are fairly evenly spaced throughout the period of leaf color change, corresponding to 20% (RMSD 6.9 days), 40% (RMSD 6.4 days), and 80% (RMSD 6.3 days) leaf color, respectively. RCC-EOF, representing the date of the maximum RCC value in autumn, was closest to 10% leaf fall (RMSD 3.7 days), the earliest abscission metric from in situ observations (Figure 3F).
GCC values in autumn were negatively correlated with leaf color percentages (average r = −0.92, 26 of 30 trees p < 0.01, Table S1). While GCC for deciduous trees generally declines as leaves senesce in autumn, RCC increases to a maximum before declining to the dormant season level. We distinguished between these increasing and decreasing phases when examining correlations between RCC and autumn phenological processes. RCC values up to and including the fall maximum correlated with leaf color percentages (average r = 0.94, all trees p < 0.01, Table S3), while fall RCC values including the maximum and later values were correlated with leaf fall percentages (tree average r = −0.94, 9 of 16 trees with sufficient observations p < 0.01). Similar to the springtime analysis, we found significant correlation when pooling all trees: leaf color was correlated to GCC, which decreased in autumn (r = −0.89, p < 0.001), as well as RCC as it increased to the fall maximum (r = 0.84, p < 0.001). Leaf fall was more highly correlated to RCC as it decreased after the fall maximum (r = −0.89, p < 0.001) than GCC (r = −0.56, p < 0.001), since GCC decreased due to color change, before leaves began falling. As in spring, we noted significant differences in regression slopes and intercepts between trees of the same species in these comparisons.
3.3. Linking Leaf Area to Color Indices
We derived transition dates (SOS, MOS, etc.) from PAI time series, as well as from drone imagery using ROIs corresponding to the top-of-canopy area analyzed for PAI (example microsite shown in Figure 4). This allowed us to examine the relation between canopy structural development and canopy color development in drone images. The timing of the spring increase in PAI and GCC was similar across the eight microsites, particularly for later spring transitions. At the first spring transition, SOS, the average difference between PAI and GCC derived dates was −4.1 (negative sign indicates GCC-MOS was earlier) ±3.2 days standard deviation, and was significantly different from zero according to a paired t-test (p < 0.05). However, the MOS and EOS transitions were closer in time, with −0.6 ± 2.5 days and 2.6 ± 4.6 days differences, respectively. Both of these differences were not significantly different from zero. In autumn, declines in GCC substantially preceded declines in PAI, as leaves changed color, becoming less green, before they were shed, decreasing the canopy PAI. As a

Sensors 2017, 17, 2852

10 of 17

result, PAI and SGensCorCs 20t1r7,a1n7, s28i5t2ions were all significantly different in autumn, with ave10roaf 1g6e differences of

12–15 daSyesns.ors 2017, 17, 2852 1.5

1.5 10 of 16

normalized PAI normalized PAI normalized GCC
normalized GCC

1.5 1

1 1.5

0.5 1

0.5 1

00

0.5 0.5

-0.5 -0.5

0 100

150

200 DOY

250

300 0

-0F.i5gure 4. Smoothed (curve) and raw (× symbols) PAI and GCC curves for an example microsite, a-n0d.5 th1e0d0ates derived from1e5a0ch data source (ve2r0t0ical lines show the2S5O0 S, MOS, EOS, SO3F0,0MOF, and EOF dates). Values are normalized to have minimuDmO0Yand maximum 1 according to the curve fits, to

FtEhiOgeFudrdaeat4et.seFtdfhaisaScgdet)imeu.lcadeisarstarV)oeal.wticteoavu4VesifTl.tlteaatauhdhoShcdltiieumeeeenlfridesdtuifcomasvrortotaeetho(miathcrdchertmreehpeueofcreardsornnocriievrotoin(morescamreevosrumcle)epnma(rhsaaTatvatoicrliaeoaniihgfdsbz)nlaodPlediaatneozAdenatrfSoetadIataf4tPdoahP)wrAss.neAaohotSIdwuIoapc(utavroG×ornc(hrien×CnedGcnaCgmses(GCevvtyyCcGriCemn(etamCairCvijonttCebmmrienbdccaowrtuajoleioRbstalmnclr)iseCtlsicoiitCPehn)wmr0asAviiee.PlgeaassIuheAln.lsunalimhdynneIoscdcemaw,oo0snplGraotroaxdhCsreoinChelmlaiGecdnoStuduedOCwdmrmiavSCccte,eora1tsMscohxPuasfaOeAisocnrmcrSmIdvSo,a(ureiOEpnctdshOmr<iSeoenfSx0s,og,.aic10tMSramet0Onoas1apoOc)F(ntlpFhcwe,SiyMeogemh,urxeclOieEdrnuaceaFrrOimpf,ovn5oas)eaSgp,ionrt,faledeliitsea,tnSoEs,agwO,mOnwtatedhFoFlliell,cerMcousOirtvFe,e,

and and fits,

to facilitate tthheedactoa,maspwaerlilsaosnatoinfdPivAidIuaalnmdicGroCsiCtestr(ap j<e0c.t0o01rifeosr .six of eight microsites). However, there were To sfiugnrtifhicearntindvifefesrteignacetes itnhreegcroesnsnioenctsiloonpesbeatnwd einetnercceoplotsrbientwdeiceensmaincrdositthese (tc-atensot pwyithlepaf< 0a.r0e5a)., we

calculateTdhitshies cpoartrieallaytidouneotfo PthAeIftaoctGthCaCt amnadxiRmCuCmvaPlAuIesv,apluoeoslveadryacmrosres tmhaicnromsaixteims u(FmigGuCrCev5a)l,uaes, well

To fuasrtwhiethrinisnpmavtiieaclsrloytistghitraeostue(gThtahobueltetchSoe4nf)o.nrSeepsctr.tiTinohgenGmbCaCexitwmwausemehsniugmhclmoyleocrtorimrirneelPdaAtiecIdveastlouaePnsAwdIet(ohpbe<se0crv.a0en0d1oa)tpwthyhelemeniacpfroasoirtleeinsag, awlle calculated

tmahsiewcrceoolsrlirateestsTslhsipahageait(ntstdiTiiioafainaislntclbdaayp,vowelnoiavtaenavhthrfirswdiiitrSdcewidPeaoIhadi4nrnuAfuledcf)lfalgfeyre,.allaIoshrfflalSoomedolots,lnmropulumwo3acletir.wteaienG4ticrsiteihntnkenoraoCxiesdgogrnaC6tfictsmoh.ovho4rGiueaerie,ritdnrrewgeCnebfesouraslhedCtdbaeac.gits(sllttTiieepwshnoRmrhtimennvoh<eCaiiacinaancntmCosxtrg0iroisopmr.mhalnovos0eoxfuioliat0aptaighlmtemxte1lehaishdiuosetmuGf(nldtepamodCyhuabCnseer<tmec,davscptr0suawpoaee.iriPcm0aaxnierorse0kAeortemdo1neoseoIrsiflffeflnPrcoevmaroaeAeratdtmutpiihiIlcsemtguteruia0asdoxehnme.cs4csbdaotiP4ntrtnefevoAmoRttoRseoawCp,sIiCPrC0iagyCseyvc.lAhv4teiarlhnm8atmenloIolduam(uuTsfoimi(eecgciairpsaascthbierrtcewrfeel<otroertsashooehss.S)etm0siah.i4oArtiten.)eebtH.s0teeswhsbm0s)sieeeon.ar1g(Hs(afwvtisa)Fxan-peoltilnewilrdewmogiiRvsnnwaetChuuggetCevm,werrrotem,whenfsirttteaG,leaheh5xtpmtaChfii)emopfsCoe,iturfciuro<avaecnremlaadolls0lw,wils.unw0ietee5gersr)ees,.ealllslaigtshnweifiidtchaatinant, differencveasriendrfersogigmrneif3sic.s4ainotcone6(.sp4l<,ow0p.0he5is)l,elaimkneadlyxiidmnuuteeminrcGpeaCprCtttvsoabariesmdtwaflrleeoremsna0mm.4p4liectsroioz0es..4iSt8iem(sTila(atbr-llytee, Sssi4tte).w-levitehl cporr<ela0t.io0n5s)o. fThis is partially

due to the facInt tfRahClaCl,tawnmde aPeAxxaIimmwieunreemdlowtPhAeer Icthovarnarelilunaetsipsornivngbaer(taywvemereangoerPerA=tIh0a.a9n0ndvmeRrsCauCsxvi0am.9lu6uefsomrfrsoGpmrCinCgthtivemaefalGlulCeRCsCa,CnsdmpPaaAxtIii)ma, lulym throughout the forestow.nhTwichahredam,aslsilfgtooahnloxilmfouiicwgmaahrnikunhtsagmlliftnohoeusefarruthbocmeeobmgrsmrieiencrlrneavotiraisnotittgniieomsano(ttefwtphtPoh<aAoet0ft.Idf0ho5eveuc(raprb)elaewuasaeisketdehsosoifuwnnfafoteuihcrtideoeuinnbmctaasonrneybosRrlepveCrayvCesadtiltniesodqaanufitscaaaftrorehtreesesaarn.etmahgAlreyiesscsbisisreinohognasisdnsip,tnareaiissnnstaggtvht,ieosawtfiHrcleiaAeelalfCdyof uffanrlodl,m 3.4 to 6.4, while maexviidmenucmeadfojGursCtlmiCneenvatrafocrorirserederilaafltraioounmtoicno0rpr.4eol4oatlietoodn dc0oa.u4tal8dan(coTrtoabsbes lumeseicSdr4od)su.ietetos,tahlethsmoualgl hsatmhperleeswizae)s. a lower statistical

significance (p < 0.05), likely due in part to a smaller sample size. Similarly, site-level correlations of

RCC and PAI were low0.e4r4thaAn although half of the microsites

in spring (average r = 0.900.5versus 0.96 for (two of four) with sufficient obseBrvations for

springtime GCC and PAI), analysis had a statistically

significant linear corr0e.l4a2tion at p < 0.05 (based on ord0in.4a5ry least squares regressions, as the HAC

adjustment for serial au0t.o4correlation could not be used due to the small sample size).

RCC

GCC

0.38
0.44 0.A36 0.42 0.34
0
0.4

24 PAI

0.4
0.5 0.35 B
0.45
60

24 PAI

6

RCC

GCC

Figure 5. ScattFGeiCgrCu0vrpae.l3lu5o8.esStscfraotatmenrddprloonltseiniamnedasgleionryefsinbofsepbsreitsntgfif(ittNbb=ete4w8t,eweeingehPetAn0mI.i4vcParAoluseiItsesvf)r;oa(mBlu)uRepCswC vafarrdloupemshaofttuoegrprtahwpehamyraadxnidmpu(Ahm)otography and

(mAa)xGimCCumvaRluCeCsRopboiCfnsCoreli00oernfvd..mafa33-adtll64illao,,dtnfaofrr(oroNemrgnir=cmeers1os6isii,moictnferosoauewgqsruiiemttahretiiocysornnoewissnoitareirtsmse)h.poGDrroCeiiCfnonf=beegrs0ee.onr3(vt4Nramt+mioi=0cn0r.o0so.14sr3b9iee8t5*etP,wsoAeeabIerines(gretihh=nredtv0fi.cama8al4tlt,eimiodcparn<xboiysm0s.d0iub0timsee1t)tisonwa)fcn;RtedC(seCyBRanmC)nCbdtRo=htlhsC0e.e.C3Tl4fahasv+etlal lmueasxiamftuermthoef

RCC

and

the

las0t.01o7b*PsAeIr0(vr a= t0i.6o0n, p(<N20.0=5).1M6ic,rfo4osiuter-lemveliccor6roresliattieosn)s.reDp0oirftfeedriennTat2bmle Sic4.ros4ites

are indicated 6

by

distinct

symbols. The pooled-data regresPsiAonI equations are GCC = 0.34 + 0.P0A19I*PAI (r = 0.84, p < 0.001) and

RCC = 0.34 + 0.017*PAI (r = 0.60, p < 0.05). Microsite-level correlations reported in Table S4.
Figure 5. Scatter plots and lines of best fit between PAI values from upward photography and (A)

GCC values from drone imagery in spring (N = 48, eight microsites); (B) RCC values after the maximum

In fall, wRCeC ienxfaallm, foirnmeidcrotshitees cwoitrhroenleaotriomnorbe oebtswerevaetnionPsAbeItwaenendthRe CfaCll mvaaxliumeusmforfoRmCC atnhdethfealalsltRCC maximum onward, folloopbowsoeliernvda-gdtiaootnau(rrNegor=ebs1ss6ie,ornfovueaqrtumiaotiicnornotsshitaearset).GtDhCiCeff=epr0ee.n3at4km+ io0cr.f0o1sa9iut*ePstAuaIrme(rni=ndR0i.c8Ca4tC,edpin<byd0.d0ic0is1at)itneacnstdstyhRmCeCbob=les0.g.3T4ihn+ening of leaf fall, which also m0.0a1r7k*PsAtIh(re= 0b.6e0g, pin<n0.i0n5)g. Moicfrotshitee-ldeveecl croerareslaetioinns trehpeorcteadninoTpaybleleS4a.f area. As in spring, we found
evidence for linear correlation in pooled data across microsites, although there was a lower statistical

significance (p < 0.05), likely due in part to a smaller sample size. Similarly, site-level correlations of

RCC and PAI were lower than in spring (average r = 0.90 versus 0.96 for springtime GCC and PAI), although half of the microsites (two of four) with sufficient observations for analysis had a statistically

significant linear correlation at p < 0.05 (based on ordinary least squares regressions, as the HAC

adjustment for serial autocorrelation could not be used due to the small sample size).

Sensors 2017, 17, 2852

11 of 17

3.4. Microclimate Effect on Phenology
To examine landscape variability in phenology, we calculated MOS and MOF dates of 10 m grid cells (Figure S6). We used GCC-MOS, except for six grid cells with a spring RCC amplitude greater than 40% of the spring GCC amplitude, based on the results of Section 3.1, and where we observed red spring leaves. We also filtered these grid cells to require the spring amplitude in RCC to be greater than 0.02, to ensure that the color metric reflected a substantial change in redness that represented leaf growth as opposed to background noise. For these six grid cells, we used RCC-MOS. Variability across grid cells (Figure S6) was similar to that observed for these dates in in situ observations of trees (14 days range in MOS, 42 days in MOF, Figure 3), with a localized area of later spring and earlier fall phenology along a forest edge bordering a wetland, in the northwest of the study area (Figure S6).
We found that 29% of variance in MOS dates and 46% in MOF dates of 10 m grid cells was explained by spatial differences in the woody plant species assemblage (regression results in Table S5). Regression coefficients indicated a similar phenological order of species to that observed in situ: American beech was among the earliest species in spring, and red maple senescence preceded that of red oak by 23 days. In a previous analysis, we found evidence of spatial autocorrelation on length scales up to at least 30 m, in residuals from a similar regression. This was presumably due in part to fine-scale temperature variation [23]. Consequently, we examined residuals from the regression we performed here, in the 30 m by 30 m neighborhoods of the microsites where temperature measurements were made. Due to the low spatial variation in spring temperatures and a limited number of microsites, we did not detect significant relationships between microsite temperature and MOS residuals. In autumn, we found that the September mean of daily minimum temperature was correlated to the residuals of MOF dates (Bonferroni-adjusted p = 0.08, r = 0.77). We found an effect size of a five day delay (1.7 day SE) in MOF date per ◦C warmer (Figure S7).
4. Discussion
4.1. Leaf Life Cycles According to Image Metrics
4.1.1. Exploration of Budburst, Leaf Growth, Color Change, and Fall
By making in situ observations of leaf life cycle events for individual trees, we were able to determine which phenology events corresponded to transitions identified from the color index (GCC or RCC) analysis of images. Similar to Berra et al. [22], who also compared ROIs of individual trees in drone imagery to ground observations, and Ahrends et al. [9], who used individual tree ROIs from a tower-mounted camera, we found that observed budburst dates corresponded closely to the beginning of the increase in color indices in spring time, i.e., the SOS transition. But, while previous studies only focused on budburst, we also examined the relation of image metrics to in situ observations of leaf expansion, finding statistically significant correlations with this process as well. MOS and EOS were closest in time to 40% and 70% leaf expansion, respectively. Similarly, we found a higher correlation between springtime color indices and the leaf expansion process, than the progression of budburst. This indicates that while the initial increase in greenness (or redness) is a good indicator of the start of budburst, subsequent increases in index values are more representative of the leaf expansion process, as opposed to the completion of budburst.
We found that both GCC and RCC metrics were significantly correlated with leaf color change in autumn. These results varied somewhat by species as suggested by others [36], although we found that both GCC and RCC work reasonably well for all the species we observed. Red maples, which may display a variety of autumn colors from pale yellow to intense scarlet, had leaf color progressions that were more highly correlated to an increase in fall RCC (species average r = 0.93) than a decrease in GCC (−0.84). However, leaf color change in red oaks, which turn a yellow-orange color (Figure 2K), was similarly correlated to both GCC and RCC (species averages r = −0.96, 0.95, respectively). In terms of

Sensors 2017, 17, 2852

12 of 17

phenophase transition dates, we did not observe species-specific biases between in situ observations and dates calculated from image analysis (Figure 3).
Our results suggest that color indices may be used to describe the trajectories of both leaf color change and leaf fall in autumn, as well as key transition dates associated with these processes. Phenophase transitions calculated from GCC (SOF, MOF, EOF) were highly correlated to a range of points in the progress of leaf coloration (20%, 40%, and 80% leaf coloration, respectively) across individual trees. The EOF metric determined from maximum fall RCC was closest to the beginning of leaf fall, indicating that the greatest “redness” of trees occurs just as leaves are starting to abscise, although the redness of individual leaves may continue to increase after this [37]. Finally, a decline in RCC after the autumn maximum was significantly correlated with the trajectory of leaf fall.
4.1.2. Plant Area Index
Measures of the plant area index from an analysis of upward photos reinforce the results of the comparison of in situ observations with drone image metrics. We found significant correlations of PAI values with drone-derived GCC in spring. This agrees with our result that in springtime, GCC from drone images is associated with leaf expansion, which increases PAI. In general, the PAI time series and GCC time series have similar trajectories in springtime (Figures 4 and 5). Previous research using GCC from tower-mounted cameras showed that the oblique angle of the camera field of view caused an early bias in late spring phenophase transitions relative to leaf area index (LAI, closely related to PAI) measured by the LAI-2000 instrument [19]. However, the nadir perspective of drone orthophotos appears to alleviate this effect; we found evidence for a strong linear relationship between springtime GCC and PAI (Figure 5A).
We also found that the decrease in RCC after the autumn maximum was significantly linearly correlated to the decrease in autumn PAI. This agrees with the earlier result that leaf fall, which decreases PAI, is correlated with RCC at the individual tree level. Overall, our comparisons with in situ observations demonstrate the utility of drone image analysis for a thorough description of spring and autumn leaf phenology, including budburst, leaf expansion, senescence, and abscission.
4.2. Leaf Color and Color Indices: Red Spring Leaves
We observed a subpopulation of red oak trees that expressed reddish colors in young leaves. When isolating these trees for analysis, we found that a greenness index typically used in spring phenology studies of digital camera images (GCC) was not an effective characterization of leaf expansion phenology; however, a redness index was (RCC, Figure 2). Subsequently, we chose to use a combination of RCC to describe the few trees or grid cells with spring amplitude in RCC greater than 40% of spring amplitude in GCC, and GCC for all other ROIs. While this criterion was useful in our identification of trees and grid cells with red spring leaves, we do not recommend it as a general rule since digital image color index values may vary by camera according to different image sensor characteristics [25], and different trees may have varying degrees of springtime redness. Researchers who observe red spring leaves should identify an appropriate criterion for their analysis.
We examined other options for color indices, including GRVI [38], ExG [17], and Hue [39]; however, these all had similar trajectories to GCC, which failed to describe leaf expansion phenology for trees with red spring leaves (Figure S2). We also calculated an index as the sum of GCC and RCC (equivalent to 1 minus blue chromatic coordinate, defined similarly to RCC and GCC [11]), in an attempt to combine the relevant information for trees with green and red spring leaves into a unified approach. However, because of the relatively large amount of noise in this index, it was less correlated to in situ observations of leaf expansion, and was less reliable in the 10 m grid cell analysis (caption to Figure S2).

Sensors 2017, 17, 2852

13 of 17

While we found utility in using RCC to quantify the phenology of trees with red spring leaves, to our knowledge, red leaves have not presented an issue for the many recent phenology studies using greenness indices to study spring phenology at Harvard Forest and elsewhere, e.g., [9,10,19,25]. This is presumably because of the relatively small number of trees with red spring leaves in relation to those with green spring leaves (i.e., only two of the 30 trees we observed in situ), and the fact that most phenology studies using digital images typically use regions of interest that integrate several tree crowns over the canopy. Indeed, we did not experience any redness-related issues with springtime GCC for the larger ROIs we used to integrate over several trees for comparison to PAI measurements (533 m2 compared to individual tree ROIs, which were 31 m2 on average).
However, individual tree ROIs have recently been increasingly examined due to the advent of new methods in photographic studies of phenology. In crowd-sourcing [36], the time-intensive task of drawing regions of interest around many individual trees is alleviated by distributing it among many citizen scientists. In photogrammetry studies using drone imagery [22], larger study areas and an undistorted nadir view of the canopy facilitate the delineation of more individual tree crowns in imagery. As greater numbers of individual trees are analyzed, our results indicate that alternatives to greenness indices, such as the redness index we used here, should prove useful in identifying and analyzing trees with red spring leaves. We also note that besides color indices, texture analysis may yield useful approaches to phenological studies of forest canopies using digital images, as suggested in a recent exploratory study [40]. Texture analysis presents the possibility of color-independent image metrics, which may be applicable across trees with different leaf colors.
Red leaf color is more typically associated with autumn than spring in temperate forests, and is known to be associated with increased concentrations of anthocyanin pigments. However, a study on tropical trees indicated that young, immature leaves were more likely than senescent leaves to express a red color and increased anthocyanin concentrations, across different species [41]. Recent research indicated that, across diverse climates including temperate forests, young leaves of many tree and shrub species have a reddish appearance [42], including red oaks [43]. The physiological function of anthocyanins has not been settled; they have been hypothesized to serve a photoprotective role, as senescent leaves resorb nutrients [44], and may also benefit young, developing leaves in this way [41]. However, we observed differences in spring leaf color in conspecific individuals in close physical proximity and with a similar canopy position (Figure 2), suggesting that light microclimate may not be a determining factor. Our results indicate that drone photography would be a useful method for the further study of intraspecific differences in springtime leaf color.
4.3. Spatial Variance in Phenology and Relation to Microclimate
The effect of spatial temperature variation on phenology is most commonly explored over altitudinal gradients in mountainous terrain, in “space for time” studies that seek to use these gradients as a proxy for global change [45], or across large geographic distances [18,46]. However, temperature variation may also act on phenology over smaller spatial scales, due to microclimates induced by topography or canopy cover within forests [47,48]. We explored the effect of microclimate by first accounting for species composition effects, as variability in phenology among tree species is well known [26,49], and examining the relation of the residual phenological variance with air temperature.
Although we did not detect an effect of microsite temperature on spring phenology, we found evidence that warmer daily minimum temperatures in September delay senescence by five days per ◦C. These results agree with the chilling degree days hypothesis, which states that senescence is triggered when accumulated cold temperatures reach a critical threshold in autumn: colder temperatures advance senescence locally [31,46]. We note that microclimates with lower daily minimum temperatures were located along a forest edge, near a wetland to the northwest of the study area (Figure 1). This result highlights the importance of considering forest structure effects on microclimate and phenology [48], in addition to factors such as canopy position and tree developmental stage [26] when assessing the causes of phenological differences between trees.

Sensors 2017, 17, 2852

14 of 17

5. Conclusions
We presented the first synthesis of in situ observations of tree phenology over a complete growing season with an analysis of drone images of the same trees. Using digital image processing techniques, we found that greenness metrics integrating several trees were correlated with canopy structural development (PAI) in springtime, while redness metrics were correlated with PAI decrease in autumn. At the individual tree level, we determined that the onset of rising greenness in spring corresponded with the start of budburst, while continued increases in greenness were driven by leaf expansion. However, we documented intraspecific variation in spring leaf color, noting that some oak trees displayed reddish (orange, pink) leaves during leaf-out while other trees had green leaves. We found that this affected the color index (greenness versus redness) that could be used to describe the leaf expansion phenology of individual trees. In autumn, we found that the decrease of canopy greenness, as well as the increase of redness, were correlated with the percentage of leaves that had changed color on trees. The time of greatest canopy redness corresponded to the beginning of abscission, and the trajectory of leaf fall correlated with decreasing redness. These results leverage the novel method of drone photogrammetry to advance our understanding of how digital image metrics relate to foliar phenology and canopy structural development throughout the leaf life cycle, from budburst to abscission.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-8220/17/12/2852/s1.
Acknowledgments: S.K. was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program—Grant 14-EARTH14R-23. Additional support was received from the National Science Foundation under the Harvard Forest Long-Term Ecological Research Program (NSF DEB-1237491) and the Harvard Center for Geographic Analysis. We acknowledge John O’Keefe of Harvard Forest for helpful discussions on this research, and three anonymous reviewers for their comments, which resulted in many improvements to the manuscript.
Author Contributions: S.K. and A.D.R. conceived and designed the experiments. S.K. performed the experiments and analyzed the data. S.K. wrote the paper with input from A.D.R.
Conflicts of Interest: The authors declare no conflict of interest.
References
1. Lieth, H. Phenology in productivity studies. In Analysis of Temperate Forest Ecosystems; Reichle, D.E., Ed.; Springer: Berlin/Heidelberg, Germany, 1973; pp. 29–46.
2. Fu, Y.H.; Zhao, H.; Piao, S.; Peaucelle, M.; Peng, S.; Zhou, G.; Ciais, P.; Huang, M.; Menzel, A.; Peñuelas, J.; et al. Declining global warming effects on the phenology of spring leaf unfolding. Nature 2015, 526, 104–107. [CrossRef]
3. Vitasse, Y.; Bresson, C.C.; Kremer, A.; Michalet, R.; Delzon, S. Quantifying phenological plasticity to temperature in two temperate tree species. Funct. Ecol. 2010, 24, 1211–1218. [CrossRef]
4. Keenan, T.F.; Gray, J.; Friedl, M.A.; Toomey, M.; Bohrer, G.; Hollinger, D.Y.; Munger, J.W.; O’Keefe, J.; Schmid, H.P.; Wing, I.S.; et al. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. Nat. Clim. Chang. 2014, 4, 598–604. [CrossRef]
5. Westergaard-Nielsen, A.; Lund, M.; Hansen, B.U.; Tamstorf, M.P. Camera derived vegetation greenness index as proxy for gross primary production in a low Arctic wetland area. ISPRS J. Photogramm. Remote Sens. 2013, 86, 89–99. [CrossRef]
6. Toomey, M.; Friedl, M.A.; Frolking, S.; Hufkens, K.; Klosterman, S.; Sonnentag, O.; Baldocchi, D.D.; Bernacchi, C.J.; Biraud, S.C.; Bohrer, G.; et al. Greenness indices from digital cameras predict the timing and seasonal dynamics of canopy-scale photosynthesis. Ecol. Appl. 2015, 25, 99–115. [CrossRef]
7. Fitzjarrald, D.; Acevedo, O. Climatic consequences of leaf presence in the eastern United States. J. Clim. 2001, 14, 598–614. [CrossRef]

Sensors 2017, 17, 2852

15 of 17

8. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol. 2013, 169, 156–173. [CrossRef]
9. Ahrends, H.E.; Bruegger, R.; Stoeckli, R.; Schenk, J.; Michna, P.; Jeanneret, F.; Wanner, H.; Eugster, W. Quantitative phenological observations of a mixed beech forest in northern Switzerland with digital photography. J. Geophys. Res. 2008, 113. [CrossRef]
10. Richardson, A.D.; Jenkins, J.P.; Braswell, B.H.; Hollinger, D.Y.; Ollinger, S.V.; Smith, M.-L. Use of digital webcam images to track spring green-up in a deciduous broadleaf forest. Oecologia 2007, 152, 323–334. [CrossRef]
11. Richardson, A.D.; Braswell, B.H.; Hollinger, D.Y.; Jenkins, J.P.; Ollinger, S.V. Near-surface remote sensing of spatial and temporal variation in canopy phenology. Ecol. Appl. 2009, 19, 1417–1428. [CrossRef]
12. Reed, B.C.; Schwartz, M.D.; Xiao, X. Remote sensing phenology. In Phenology of Ecosystem Processes; Springer: New York, NY, USA, 2009; pp. 231–246.
13. Garrity, S.R.; Bohrer, G.; Maurer, K.D.; Mueller, K.L.; Vogel, C.S.; Curtis, P.S. A comparison of multiple phenology data sources for estimating seasonal transitions in deciduous forest carbon exchange. Agric. For. Meteorol. 2011, 151, 1741–1752. [CrossRef]
14. Lange, M.; Dechant, B.; Rebmann, C.; Vohland, M.; Cuntz, M.; Doktor, D. Validating MODIS and sentinel-2 NDVI products at a temperate deciduous forest site using two independent ground-based sensors. Sensors 2017, 17, 1855. [CrossRef]
15. Zhang, X.; Friedl, M.A.; Schaaf, C.B.; Strahler, A.H.; Hodges, J.C.F.; Gao, F.; Reed, B.C.; Huete, A. Monitoring vegetation phenology using MODIS. Remote Sens. Environ. 2003, 84, 471–475. [CrossRef]
16. Elmore, A.J.; Guinn, S.M.; Minsley, B.J.; Richardson, A.D. Landscape controls on the timing of spring, autumn, and growing season length in mid-Atlantic forests. Glob. Chang. Biol. 2012, 18, 656–674. [CrossRef]
17. Hufkens, K.; Friedl, M.; Sonnentag, O.; Braswell, B.H.; Milliman, T.; Richardson, A.D. Linking near-surface and satellite remote sensing measurements of deciduous broadleaf forest phenology. Remote Sens. Environ. 2012, 117, 307–321. [CrossRef]
18. Klosterman, S.T.; Hufkens, K.; Gray, J.M.; Melaas, E.; Sonnentag, O.; Lavine, I.; Mitchell, L.; Norman, R.; Friedl, M.A.; Richardson, A.D. Evaluating remote sensing of deciduous forest phenology at multiple spatial scales using PhenoCam imagery. Biogeosciences 2014, 11, 4305–4320. [CrossRef]
19. Keenan, T.F.; Darby, B.; Felts, E.; Sonnentag, O.; Friedl, M.A.; Hufkens, K.; O’Keefe, J.; Klosterman, S.; Munger, J.W.; Toomey, M.; et al. Tracking forest phenology and seasonal physiology using digital repeat photography: A critical assessment. Ecol. Appl. 2014, 24, 1478–1489. [CrossRef]
20. Dandois, J.P.; Ellis, E.C. High spatial resolution three-dimensional mapping of vegetation spectral dynamics using computer vision. Remote Sens. Environ. 2013, 136, 259–276. [CrossRef]
21. Lisein, J.; Michez, A.; Claessens, H.; Lejeune, P. Discrimination of deciduous tree species from time series of unmanned aerial system imagery. PLoS ONE 2015, 10, e0141006. [CrossRef]
22. Berra, E.F.; Gaulton, R.; Barr, S. Use of a digital camera onboard a UAV to monitor spring phenology at individual tree level. In Proceedings of the 2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS), Beijing, China, 10–15 July 2016; pp. 3496–3499.
23. Klosterman, S.; Melaas, E.; Wang, J.; Martinez, A.; Frederick, S.; O’Keefe, J.; Orwig, D.A.; Wang, Z.; Sun, Q.; Schaaf, C.; et al. Fine-scale perspectives on landscape phenology from unmanned aerial vehicle (UAV) photography. Agric. For. Meteorol. 2018, 248, 397–407. [CrossRef]
24. Klosterman, S.; Richardson, A.D. Landscape Phenology from Unmanned Aerial Vehicle Photography at Harvard Forest since 2013. Harvard Forest Data Archive: HF294. Available online: http://harvardforest.fas. harvard.edu:8080/exist/apps/datasets/showData.html?id=hf294 (accessed on 1 December 2017).
25. Sonnentag, O.; Hufkens, K.; Teshera-Sterne, C.; Young, A.M.; Friedl, M.; Braswell, B.H.; Milliman, T.; O’Keefe, J.; Richardson, A.D. Digital repeat photography for phenological research in forest ecosystems. Agric. For. Meteorol. 2012, 152, 159–177. [CrossRef]
26. Richardson, A.D.; O’Keefe, J. Phenological differences between understory and overstory. In Phenology of Ecosystem Processes; Noormets, A., Ed.; Springer: New York, NY, USA, 2009; pp. 87–117.

Sensors 2017, 17, 2852

16 of 17

27. Andrews, D.W.K. Heteroskedasticity and autocorrelation consistent covariance matrix estimation. Econometrica 1991, 59, 817–858. [CrossRef]
28. Andrews, D.W.K.; Monahan, J.C. An improved heteroskedasticity and autocorrelation consistent covariance matrix estimator. Econometrica 1992, 60, 953–966. [CrossRef]
29. Macfarlane, C.; Hoffman, M.; Eamus, D.; Kerp, N.; Higginson, S.; Mcmurtrie, R.; Adams, M. Estimation of leaf area index in eucalypt forest using digital photography. Agric. For. Meteorol. 2007, 143, 176–188. [CrossRef]
30. Macfarlane, C. Classification method of mixed pixels does not affect canopy metrics from digital images of forest overstorey. Agric. For. Meteorol. 2011, 151, 833–840. [CrossRef]
31. Richardson, A.D.; Bailey, A.S.; Denny, E.G.; Martin, C.W.; O’Keefe, J. Phenology of a northern hardwood forest canopy. Glob. Chang. Biol. 2006, 12, 1174–1188. [CrossRef]
32. Vitasse, Y.; François, C.; Delpierre, N.; Dufrêne, E.; Kremer, A.; Chuine, I.; Delzon, S. Assessing the effects of climate change on the phenology of European temperate trees. Agric. For. Meteorol. 2011, 151, 969–980. [CrossRef]
33. Friedl, M.A.; Gray, J.M.; Melaas, E.K.; Richardson, A.D.; Hufkens, K.; Keenan, T.F.; Bailey, A.; O’Keefe, J. A tale of two springs: Using recent climate anomalies to characterize the sensitivity of temperate forest phenology to climate change. Environ. Res. Lett. 2014, 9. [CrossRef]
34. Orwig, D.; Foster, D.; Ellison, A. Harvard Forest CTFS-ForestGEO Mapped Forest Plot since 2014. Available online: http://harvardforest.fas.harvard.edu:8080/exist/apps/datasets/showData.html?id=hf253 (accessed on 25 August 2017).
35. Draper, N.R.; Smith, H. Applied Regression Analysis, 3rd ed.; John Wiley & Sons: New York, NY, USA, 1998. 36. Kosmala, M.; Crall, A.; Cheng, R.; Hufkens, K.; Henderson, S.; Richardson, A. Season Spotter: Using citizen
science to validate and scale plant phenology from near-surface remote sensing. Remote Sens. 2016, 8, 726. [CrossRef] 37. Archetti, M.; Richardson, A.D.; O’Keefe, J.; Delpierre, N. Predicting climate change impacts on the amount and duration of autumn colors in a new england forest. PLoS ONE 2013, 8, e57373. [CrossRef] 38. Motohka, T.; Nasahara, K.N.; Oguma, H.; Tsuchida, S. Applicability of green-red vegetation index for remote sensing of vegetation phenology. Remote Sens. 2010, 2, 2369–2387. [CrossRef] 39. Nagai, S.; Maeda, T.; Gamo, M.; Muraoka, H.; Suzuki, R.; Nasahara, K.N. Using digital camera images to detect canopy condition of deciduous broad-leaved trees. Plant Ecol. Divers. 2011, 4, 79–89. [CrossRef] 40. Almeida, J.; dos Santos, J.A.; Alberton, B.; Torres, R.S.; Morellato, L.P.C. Applying machine learning based on multiscale classifiers to detect remote phenology patterns in Cerrado savanna trees. Ecol. Inform. 2014, 23, 49–61. [CrossRef] 41. Lee, D.W.; Collins, T.M. Phylogenetic and ontogenetic influences on the distribution of anthocyanins and betacyanins in leaves of tropical plants. Int. J. Plant Sci. 2001, 162, 1141–1153. [CrossRef] 42. Lev-Yadun, S.; Yamazaki, K.; Holopainen, J.K.; Sinkkonen, A. Spring versus autumn leaf colours: Evidence for different selective agents and evolution in various species and floras. Flora-Morphol. Distrib. Funct. Ecol. Plants 2012, 207, 80–85. [CrossRef] 43. Key, T.; Warner, T.; McGraw, J.; Fajvan, M. A comparison of multispectral and multitemporal information in high spatial resolution imagery for classification of individual tree species in a temperate hardwood forest. Remote Sens. Environ. 2001, 75, 100–112. [CrossRef] 44. Lee, D.W.; O’Keefe, J.; Holbrook, N.M.; Feild, T.S. Pigment dynamics and autumn leaf senescence in a New England deciduous forest, eastern USA. Ecol. Res. 2003, 18, 677–694. [CrossRef] 45. Vitasse, Y.; Porté, A.J.; Kremer, A.; Michalet, R.; Delzon, S. Responses of canopy duration to temperature changes in four temperate tree species: Relative contributions of spring and autumn leaf phenology. Oecologia 2009, 161, 187–198. [CrossRef] 46. Dragoni, D.; Rahman, A.F. Trends in fall phenology across the deciduous forests of the Eastern USA. Agric. For. Meteorol. 2012, 157, 96–105. [CrossRef] 47. Fisher, J.; Mustard, J.; Vadeboncoeur, M. Green leaf phenology at Landsat resolution: Scaling from the field to the satellite. Remote Sens. Environ. 2006, 100, 265–279. [CrossRef]

Sensors 2017, 17, 2852

17 of 17

48. Motzkin, G.; Ciccarello, S.C.; Foster, D.R. Frost pockets on a level sand plain: Does variation in microclimate help maintain persistent vegetation patterns? J. Torrey Bot. Soc. 2002, 129, 154–163. [CrossRef]
49. Lechowicz, M.J. Why do temperate deciduous trees leaf out at different times? Adaptation and ecology of forest communities. Am. Nat. 1984, 124, 821–842. [CrossRef]
© 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).