Phenological Dynamics Characterization of Alignment Trees with Sentinel-2 Imagery: A Vegetation Indices Time Series Reconstruction Methodology Adapted to Urban Areas

: This article presents a novel methodology for the characterization of tree vegetation phenology, based on vegetation indices time series reconstruction and adapted to urban areas. The methodology is based on a pixel by pixel curve ﬁtting classiﬁcation, together with a subsequent Savitzky–Golay ﬁltering of raw phenological curves from pixels classiﬁed as vegetation. Moreover, the new method is conceived to face speciﬁcities of urban environments such as: the high heterogeneity of impervious/natural elements, the 3D structure of the city inducing shadows, the restricted spatial extent of individual tree crowns and the strong biodiversity of urban vegetation. Three vegetation indices have been studied: Normalized Difference Vegetation Index (NDVI) and Normalized Difference Red Edge Index 1 (NDRE1), which are mainly linked to chlorophyll content and leaf density and Normalized Burn Ratio (NBR) mostly correlated to water content and leaf density. The methodology has been designed to allow the analysis of annual and intra-annual vegetation phenological dynamics. Then, different annual and intra-annual criteria for phenology characterization are proposed and criticized. To show the applicability of the methodology, this article focuses on Sentinel-2 (S-2) imagery covering 2018 and the study of groups of London planes in an alignment structure in the French city of Toulouse. Results showed that the new method allows the ability to (1) describe the heterogeneity of phenologies from London planes exposed to different environmental conditions (urban canyons, proximity with a source of water) and (2) to detect intra-annual phenological dynamics linked to changes in meteorological conditions.


Introduction
The presence of vegetation improves the responsible and sustainable development of urban areas by offering several ecosystem services [1][2][3] such as: vegetation-climate interactions (air quality increase, urban cool island, rain water management and thermal comfort) [4,5], energetic consumption and CO2 track reductions [6], biodiversity conservation [6], human well-being [7], ambiance and socio-cultural benefits [7]. Among this urban vegetation, alignment trees create shade contributing to reduce the urban heat islands and fit the geometrical urban planning [8]. As an example in Paris, about 100,000 street trees cover around 700 km of roads and concern approximately 1600 roads out of 6000, inducing a shade covering about 3% of the city surface [9]. However, urban trees are exposed to increasing anthropic and climatic pressures [10,11]. Hence, monitoring their health is essential to protect green areas in cities, anticipate changes in climatic forcing and guide environmental policies [12]. Ground surveys performed by members of the city hall departments managing the green areas of the city, or by external experts, have been traditionally used [13]. However, they are costly, require human resources and suffer a lack of high temporal repeatability to give an accurate report of health condition and its evolution. Remote sensing data is a good alternative and has been already largely used in urban environments for vegetation classification [14][15][16][17][18][19] and health assessment [20][21][22][23], with acquisitions from unmanned aerial vehicle sensors [15,24], airborne [14,20,[25][26][27] and spaceborne sensors [16][17][18]28]. However, most of these previous studies work on a limited number of acquisitions along the year (often only one date) and as a consequence they cannot properly differentiate normal vegetation dynamics from a change in health condition.
Actually, vegetation dynamics are essentially controlled by the annual phenology, which depends on the climate seasonal variations determining key phases like greenup, maturity, senescence and dormancy periods and also depends on the biome type [29,30]. In addition, each vegetation species and individual has a specific phenology. Then, stress events can be detected by its intra-annual and/or inter-annual fluctuations. Studying vegetation phenology from remote sensing data requires a high temporal revisit, thus being mainly restricted to the use of satellite imagery. Thereafter, the term "phenology" frequently refers to the analysis of time series of spectral Vegetation Indices (VIs) computed from satellite images. The Normalized Difference Vegetation Index (NDVI) is the most used VI correlated to green leaf density (LAI), green biomass and green vegetation cover and has been historically obtained from AVHRR and MODIS sensors [29,[31][32][33][34]. Although both of them take daily images, their spatial resolutions are coarse (NDVI is obtained at 500 m and 250 m, respectively), and only allow capturing the global phenology of vegetation communities from regional to national scale, but not at the local or individual scale needed in urban environments. On the other hand, other satellite sensors have a better spatial resolution such as the Landsat series (ground sample distance of 30 m), SPOT (10 m) and RapidEye (5 m) constellations, and with them more heterogeneous ecosystems can be studied [35,36]. However for Landsat sensors, the spatial resolution is still not adapted to urban tree crown dimensions and their lower temporal revisit of around 16 days considerably reduces the number of cloud free acquisitions to integrate for building the phenological curve [35]. The use of SPOT and RapidEye appears to be more adapted to study urban tree phenology with their higher temporal revisit between 1 and 5 days. However, their images are not freely available and this can represent a significant cost for accurate monitoring. Recently, new multi-spectral satellite sensors such as the European Sentinel-2 and the Israeli-French Venµs with a spatial resolution of 10/20 m and 10 m, and a temporal revisit of 5 and 2 days respectively, have the advantage to acquire information in a larger number of spectral bands (more than the traditional Red Green Blue (RGB) + Near InfraRed (NIR)). This opens the way to compute more VIs related to other vegetation bio-physico-chemical and structural variables of interest for tree health monitoring [37][38][39][40][41], for instance with leaf pigments content from the visible (VIS) and NIR spectral domains [42], and leaf water and dry matter content in the Short-Wave InfraRed (SWIR) domain [42].
In addition, VI time series from satellite imagery contain noise induced by atmospheric variability, cloud cover variability and illumination and observation geometries variability [43]. Consequently, dedicated techniques for noise reduction were developed and can be split into three different approaches: function fitting methods [44,45], harmonic analysis methods [46,47] and local filtering methods [48][49][50]. The first approach consists of fitting the measured VI time series with a fixed mathematical function, assuming that healthy vegetation phenology behaves as a given function of time. The main disadvantage is that this approach prevents us from observing possible intra-annual phenological dynamics. Several fixed functions have been tested: double logistic (D-L) function [51,52], double hyperbolic tangent (DHT) function [45] or assymetric Gaussian function [51]. Then, the harmonic analysis approach, based on modeling the VI curve by series of sines and cosines, hypothesizes a symmetrical behavior of phenology, and so is not able to capture asymmetrical phenological dynamics. Some methods belonging to this approach are: HANTS [53], Fourier Filtered Cycle Similarity [47] or the Sellers et al. 1994 Fast Fourier Transform technique [54]. Finally, local filtering approach consists of time localized smoothing of the phenological curves constrained to some given hypotheses. Among the most used filters are the Savitzky-Golay (S-G) filter [50] and Best Index Slope Extraction (BISE) methods [48,49]. This last approach family allows us to observe asymmetries and intra-annual events, but the smoothing hypotheses should be well chosen. Finally, Vrieling et al. 2018 showed the importance of the number of available images to correctly characterize phenological dynamics with reconstruction methods (in their case D-L fitting) [45]. For instance, the performance of the phenology reconstruction methodology decreases when the number of available images decreases (for them from 27 to 15). In addition, they demonstrated that the distribution of the images across the year is also important, being necessary to sample greenup and senescence periods as well as minimal and maximal values (dormancy and maturity periods).
Previous studies characterizing the phenological dynamics of vegetation with the use of VIs derived from satellite imagery focus on natural (woodlands and wild areas) or agricultural environments [35,45,[50][51][52], with urban environments remaining scarcely investigated [55][56][57]. In addition, they are mainly based on NDVI phenology curves [35,[50][51][52] leaving the rest of VIs only slightly analyzed [45], and so, the spectral richness of new satellites is not fully taken into account. The lack of studies focused on the monitoring of urban tree phenological dynamics can be explained by the additional issues due to the specificities of urban environments such as: the high heterogeneity of impervious/natural elements and the high geometry variability of the city, as well as its 3D structure heterogeneity inducing shadow effects, the restricted spatial extent of individual tree crowns (around ten of meters) and the strong biodiversity of urban vegetation. Thus, studying phenology curves of urban vegetation, not only needs satellites with a high temporal resolution, but also satellites with high spatial resolutions (a few tens of meters). Nowadays, missions such as the European Sentinel-2 or the Israeli-French Venµs [58] seem to be appropriate for urban areas.
The objective of this article is to propose a new methodology adapted to urban areas to characterize tree phenology by using VI time series reconstruction from Sentinel-2 imagery. Globally, it is based on (1) a classification of vegetation pixels relying on phenological curve fitting, considering healthy tree phenology to approximately behave as a given fixed function and (2) the application of a filter on raw VI time series for pixels classified as vegetation. This methodology has been developed to successfully: face high heterogeneity of urban environments and shadows (challenge of mixed pixels), face atmospheric, cloud cover, illumination-observation geometry and registration variability between dates and characterize annual and intra-annual phenological dynamics. A first good study case to test the new method has been selected based on criteria such as the choice of mono-species trees (to avoid inter-species phenology variability which can hide or generate mixed phenological dynamics) with a large tree crown size (in the order of magnitude of the satellite spatial resolution) and in alignment structure (to facilitate the number of vegetation pixels). The London plane is a common deciduous broadleaf tree species in European cities like Brussels, London, Madrid and Toulouse, and its monitoring is extremely important since some trees have been sensitive to canker stain disease and dying for the last 70 years in Italy, France, United kingdom, Spain and Greece [59,60]. So, the study case for this article is focused on alignment London plane trees in the city of Toulouse, France, with a monitoring of their phenology over the year 2018 with Sentinel-2 imagery (2018 is the first complete year for which S-2 has a revisit of 5 days). For that purpose, different environmental conditions (urban canyons, the proximity with a source of water) for urban and peri-urban London planes are investigated to observe variations in phenological curves as described by three distinct VIs. The proposed methodology finally allows us to study both annual and intra-annual phenological characteristics. The second is defined as anomalous events in the global annual dynamic [35,61,62] possibly corresponding to stress/disturbance events when analyzed with meteorological information.
Section 2 describes the studied area, the satellite and climatological data used in the study. In Section 3, the proposed methodology is presented together with the annual and intra-annual criteria for phenology monitoring. Results are showed in Section 4 for the study case and the different environment scenarios for the London plane trees. Section 5 discusses the appropriateness of the chosen hypotheses of the methodology formulation, the advantages and disadvantages of the criteria to characterize the phenology and the influences of meteorological factors on phenological curves. Finally, Section 6 presents conclusions and perspectives of this work.

Case Study
Toulouse is located in the south-west of France, close to the Pyrenees mountain range and between Atlantic ocean and Mediterranean sea. With around 500,000 inhabitants, it is the fourth largest city of the country. It is characterized by a humid temperate climate favoring vegetation development. Toulouse is crossed by the Garonne river and by the historical "Canal de Midi" and "Canal de Brienne" artificial canals. London planes are found in large avenues in the city center, but also along the river and the canals.
In this work, five different areas of Toulouse, where London planes can be found exposed to different environmental conditions, have been studied (see Figure 1): (1) Canal de Brienne (Brienne): These old London planes are located near the city center and grow along a canal. The orientation of this alignment is north-west, with two rows of trees (one at each side of the canal) and crown widths between 10 and 20 m. London planes from the city center (Brienne, Midi Int, JGuesde and FVerdier) are supposed to be exposed to pollution and high temperatures during summer due to the urban heat island effect [63], while London planes from Midi Ext do not. In addition, London planes from city center avenues (JGuesde and FVerdier) may be exposed to lack of water during droughts. This is not the case for London planes from the canals (Brienne, Midi Int and Midi Ext), which do not suffer from lack of water since water from canals moisten the soil under them, and in some cases roots reach directly the canal water. Furthermore, London planes from the canals (Brienne, Midi Int and Midi Ext) are rooted on wide soil surfaces, while planes from the city center avenues (JGuesde and FVerdier) are rooted on small soil squares surrounded by impervious materials such as asphalt or slabs (see Figure 1). At last, London planes from JGuesde also suffer from the impacts of their root system alteration and a soil replacement due to 2011-2013 road works in order to extend a tramway line (information from Town-Hall).   Table A1. In addition, the viewing angle of both Sentinel-2 satellites is almost constant across the year and very close to nadir (viewing zenithal angle θ v z ≈ 3 • and viewing azimuthal angle φ v a ≈ 202 • ). In Toulouse, the illumination angles (zenithal θ i z and azimuthal φ i a angles) at S-2 acquisition time (≈11:00 local time) varies along the year going from θ i z ≈ 165 • and φ i a ≈ 65 • during winter to θ i z ≈ 145 • and φ i a ≈ 25 • during summer, see Figure 2.  THEIA platform (https:www.theia-land.fr) from CNES (Centre National d'Etudes Spatiaux (France)) gives free access to Sentinel-2 Bottom Of Atmosphere reflectance (level 2A) processed images of Toulouse. Moreover, THEIA provides co-registered images, with cloud masks, and quality and atmospheric information. For Sentinel-2A and Sentinel-2B combined products, the multi-temporal spatial registration performance varies between pixels, with approximately a third of pixels with co-registration errors of 0-50% of the pixel size, another third with errors around 100-150% of the pixel size, and the last third with errors in between [64]. For 2018, THEIA provides up to 37 S-2 images of Toulouse, among which there are 20 exploitable ones (images with strong cloud coverage, higher than 50% or hidden the studied areas, are not considered). These 20 exploitable images are distributed across the year in such a way that they sample both dormancy and maturity period, as well as greenup and senescence, thus allowing us to correctly apply phenology reconstruction methodologies.
Three VIs accessible with Sentinel-2 are studied in this work: NDVI and Normalized Difference Red Edge Index 1 (NDRE1) were considered to give information on leaf density and chlorophyll content [65] which are measured at 10 m and 20 m resolution respectively, and Normalized Burn Ratio (NBR), which is supposed to provide information on leaf density and water content [65] at 20 m resolution, see Table 1. Normalized indices are selected to facilitate quantification, comparison and description of phenological dynamics. The meteorological data used over the year 2018 come from the Meteopole-Flux station installed on the Meteo France site 6.5 km south-west of the city center, see Figure 1. This permanent station has been operational since 2012 and allows the long-term monitoring of radiation and energy exchanges, meteorological variables and soil water and thermal status, for a grassland site on the outskirts of the Toulouse urban area. For this study, the data analyzed are the incoming short-wave radiation (in W·m −2 ) and the air temperature (in • C) recorded at 2 m above the ground, the soil water content (in m 3 ·m −3 ) recorded by 16 sensors installed at different depths (from 10 to 220 cm depth) and precipitation (mm). Data are hourly available and are here post-processed to provide daily information of mean global incoming radiation, maximum temperature and cumulative precipitation. For water content, mean daily values are calculated by soil layers i.e., 10-50 cm, 50-100 cm and 100-220 cm, because the temporal dynamics of soil water status (especially its response to daily and subdaily precipitation) varies between near-surface layers and deep soil.

Phenology Time Series Reconstruction
The proposed methodology presents four steps, see Figure 3. (1) A manual masking aims at discriminating the vegetation study areas (in this study five masks are created, see Section 2.1). This mask is visually delineated by using S-2 images but also georeferenced open access airborne images of the city (Google Earth images). Only one mask per green studied area was used for the whole studied year. So, to entirely contain the vegetated area, the mask should be set during the period of the year for which vegetation occupies the largest area, both in terms of quantity of pixels and quantity of vegetation per pixel (this period usually corresponds to the maturity period). For this work, QGIS software was used to create the masks.
The chosen VI was calculated for each pixel of the masked areas and for all the available dates of the year. Hence for each pixel, a raw VI time series, V I 0 t , with a number of samples equal to the number of available images, was obtained.
Then, an unsupervised classification based on pixel by pixel weighted iterative fitting of phenological curves was applied. Weighted iterative curve fitting has been traditionally used to study vegetation phenology, by considering that healthy vegetation presents phenological curves that can be described by fixed functions [45,51,52]. The proposed methodology takes advantage of this hypothesis and considers that pixels where the error fit is large are mixed or non-vegetated pixels. On the other hand, pixels where the error fit is considerably small, i.e., behaving as the imposed function, are considered vegetation pixels. In this work, double logistic function was chosen to describe phenology time series of VIs, since it presents less free parameters than the double hyperbolic tangent function and it has been shown to perform better than the asymmetric Gaussian function [51]. The double logistic function is expressed as: where V I min , V I max , S, A, m S and m A are the free parameters of the fit. With V I min the minimum yearly value of the VI, V I max its maximum yearly value, S the Day of the Year (DoY) indicating the inflection point when the curve raises (greenup period), A the DoY indicating the inflection point when the curve drops (senescence period) and m S and m A the slopes of the curve at days S and A respectively.
The weighted iterative process relies on the hypothesis that in urban environments, small co-registration errors between dates are the main source of noise induced on V I 0 t time series, and thus, considers that sudden increases or drops of V I 0 t values between dates are mainly due to changes in the fraction of vegetation within a pixel. To correctly characterize the vegetation phenology and reduce mixed pixels effect, more weight should be given to high VI values.
Then from a first fit (V I 1 t ) of the pixel raw time series (V I 0 t ), a set of weights are computed to give more prominence to high VIs values than to low ones (upper envelope of the curve): where W t i is the weight of sample at date t i , V I 0 t is the raw VI time series, V I 1 t is the first fitted VI time series, Once the weights have been fixed from the first fitting, an iterative process starts until minimizing the fitting error defined as: where F k is the error at iteration k and N is the number of samples in the time series. Fits are performed with the Levenberg-Marquardt least squares method [72,73] of the SciPy library [74] for Python 3.6.
Finally, pixel classification is based on F k errors and two thresholds, the pixel is considered as mixed and if F k > B m the pixel is considered as not (or poorly) vegetated. B v and B m were empirically chosen to be respectively 5% and 10% of the annual maximal value VI of the pixel. For our study case, these bounds appear strict enough to importantly reduce non-vegetated pixels in those classified as vegetation. An additional condition was applied: among the pixels initially classified as vegetation on each masked area, only pixels whose mean VI value during the maturity period is comprised within one standard deviation (std) of the ensemble of vegetation pixels over the same period, were retained. For London planes in Toulouse, we defined this period between 1st of May and 1st of October, since for this species in Toulouse, greenup and senescence periods correspond to the beginning of spring and the beginning of autumn respectively. Background effects are major during winter (when there are no leaves) and during the start of the greenup period and the end of the senescence period. So, the choice of the maturity period bounds were done considering that from 1st of May to 1st of October the vegetation cover should be enough to minimize the influences of ground in the statistical measures. Adding this condition on the mean VI values allows us to avoid, as much as possible, false-positives in the classification. (4) A weighted iterative Savitzky-Golay filtering was applied on vegetation pixels as classified in step (3). It allows us to reduce the noise that registration variability between dates induces on time series, and as it does not fix the phenological behavior, it allows us to observe intra-annual periods of disturbance and anomalies in vegetation phenological curves. The weighted iterative Savitzky-Golay filter used in this work is our own implementation (in Python 3.6) of the filter presented in Chen et al. 2004 [50], but it was adapted to process not uniformly sampled time series. Savitzky-Golay filter presents two free parameters: the width of the smoothing window (m) and the degree of the fitting polynomial (d).
(a) First, an automatic search of the best parameter values within a given range, m ∈ [6, 10] and d ∈ [2,4], was done [50]. These ranges were fixed following [50] and considering the number of available dates. On one hand, too small m values may lead to an over-fit of the raw VI time series, and so, difficulties can raise in capturing long-term trends. However, too large m values may lead to neglecting some important variations in the phenological curve. On the other hand, small values of d may lead to smoother results, and large values of d may over-fit the raw VI time series and generate noisy results. In addition, to perform a polynomial fit, the degree d of the polynomial must be lower than the size m of the smoothing window.
(b) Then, once the first filtering was performed providing V I 1,filt From step (3), another characterization of vegetation health status can be directly obtained from the D-L fitted curves of pixels classified as vegetation. Thus in Section 4.3, this alternative methodology is compared to the proposed one. However, D-L fitting characterization does not allow us to describe intra-annual dynamics and is then not advised. On the contrary, the proposed methodology by using S-G filtering does not fix a given behavior for the phenological curve and then it allows us to reduce VI time series noise without masking intra-annual anomalies in the phenological dynamics. So, annual and intra-annual dynamics can be studied, see Sections 3.2 and 3.3.

Annual Phenology Characterization
The final VI time series obtained from the previous reconstruction methodology should be analyzed by looking at different phenological information that can be compared from one year to another [45]:

Intra-Annual Phenology Characterization
The S-G reconstruction methodology allows us to observe intra-annual dynamics, that might be related to phenological disturbance periods during the year. Here, these periods are defined as periods with anomalous dynamics (that are unpredictable with reconstruction methods imposing a function behavior): periods where the concavity/convexity of the phenological curve is suddenly broken [35,61,62]. Their quantification can be assessed with: (e) SlVI t i : Negative slope (VI/day) of the VI curve at the beginning of the phenological disturbance period, see Figure 5. This indicator provides information on the strength of the perturbation. The more negative the slope is, the more important the phenological disturbance is, i.e., the more important is the concavity break. On a global overview of the chosen annual criteria, vegetation phenology is analyzed through the maximum values of the VI (V I max ), but also, by studying how long vegetation maintains high VI values (Green period , PS90 S , PS90 E ), how fast is greenup (SOS20,SOS50) or when senescence arrives (EOS20, EOS50). Two additional criteria are also proposed to quantify intra-annual anomalous phenology events (slVI t i and DiffA).

Results
This section studies the performances of the above presented methodology when it was applied to characterize urban and peri-urban London plane phenological curves as described by NDVI, NBR and NDRE1, see Table 1. First, the performances of the pixel classification and the S-G filtering steps are shown. Then, the reconstructed time series are used to characterize the vegetation phenological dynamics. The results obtained with the S-G filtered time series are compared to those obtained with a fitting function (double logistic) time series reconstruction method applied on the pixels classified as vegetation. Only pixels in the masked areas shown in Figure 1 are processed. Figure 6 shows the classification map obtained with curve fitting on the Canal de Brienne area, with NDVI (10 m resolution), NBR (20 m resolution) and NDRE1 (20 m resolution) respectively. Based on step (3) of the methodology described in Section 3.1, vegetation pixels are colored in fuchsia, mixed pixels in cyan and non-vegetation pixels in blue. As expected, the number of vegetation pixels is higher with NDVI than with NBR or NDRE1 due to the difference in spatial resolution. Further, the location of the vegetation pixels at 20 m overlaps that of the vegetation pixels at 10 m. However, vegetation pixel locations from NBR and NDRE1 do not fully coincide. Figure 7 shows the raw VI curves for pixels classified as vegetation on the five studied areas for NDVI, NBR and NDRE1. After classification, pixels classified as vegetation present a common phenology behavior associated to vegetation in any of the studied sites. In addition, whatever the indices, the differences of the VI values between time series are smaller on Brienne, Midi Int and Midi Ext compared to JGuesde and FVerdier sites. A larger range of VIs value variations during the winter season is also noticed, whatever the location, but being larger in FVerdier. This effect is due to the background heterogeneity, since during the dormancy period, the lack of leaves increases the background influences on the phenological curve.  Table 2 shows the number of pixels classified as vegetation together with the percentage they represent in each studied area. It also shows, for each studied site, the standard deviation of the VI values of the vegetation pixels during the maturity period. For every VI, Brienne and Midi Ext present the lowest standard deviations (between 0.02 and 0.03 depending on the VI), followed by Midi Int (0.03, 0.04 and 0.05 for NDVI, NBR and NDRE1 respectively). For NDVI, the percentage of vegetation pixels is higher for Brienne and Midi Ext than over the other sites. On the other hand, for both NBR and NDRE1 the highest percentages of vegetation pixels are found for Midi Ext and JGuesde. Moreover, the percentage of pixels classified as vegetation obtained for the studied areas is more regular with NDRE1 (between 29% and 40% of pixels), while for NDVI and NBR larger variations are found (between 22% and 58% for NDVI and 13% and 46% for NBR). This can indicate that D-L fitting classification with NDRE1 is less influenced by the type of background.   Figure 8 allows to compare the NDVI (a), NBR (b) and NDRE1 (c) Savitzky-Golay filtered time series for pixels from Canal de Brienne classified as vegetation, mixed and non-vegetation pixels. Green lines represent pixels that are completely occupied by London planes. They show the highest VI values and after application of the proposed methodology, no sudden variability between dates remains. They also exhibit the smallest standard deviation. Black lines represent non-vegetated (or very few vegetated, or very unhealthy vegetation) pixels. They show the lowest VI values and, after applying the proposed methodology, noise remains important. This can be due to misregistration between dates, i.e., the composition of the pixels varies from one date to another, sometimes partially containing vegetation. In between, yellow line represents mixed pixels, with variable vegetation proportion. It can also represent unhealthy vegetation pixels, i.e., pixels for which the phenology curve can not be fitted by a given fixed function such as "double logistic" or "double hyperbolic tangent". Both, mixed and non-vegetated pixels present large standard deviations showing the existence of very different phenology behaviors in these groups. However, Figure 8 shows that the yellow means are very close to the green ones with a very similar behavior. This may indicate that the classification procedure has been restrictive. This strict classification is used to avoid false-positive vegetation pixels in the analysis. Identical qualitative results are found for the other studied areas (data not shown). Figure 8 also shows, in red, the mean of the D-L fitted curve for vegetation pixels. While for NBR green and red curves are very similar, for NDVI and NDRE1 the red line does not describe the drop appearing for the green one at DoY ≈ 205. In addition, D-L fitting, imposing a fixed function, provides smoother VIs time series, as expected.

Annual Dynamics
Two methodologies to reconstruct vegetation phenological curves are compared in this section: Savitzky-Golay filtering and double logistic fitting, both applied on pixels previously classified as vegetation. Figure 9 shows the phenology curves of London planes as characterized by NDVI (a, b), NBR (c, d) and NDRE1 (e, f) for the five studied areas described in Section 2.1. Curves in left column of Figure 9 are defined as the mean of the Savitzky-Golay filtered phenology curves, while curves in right column of Figure 9 are defined as the mean of the double logistic fitted phenology curves. In Figure 9a-d differences between the five sites can be observed with the NDVI and NBR reconstructed time series, and independently of the reconstruction methodology. These differences are quantified by V I max , and Green period criteria, see Tables 3 and 4     On the other hand, criteria such as SOS20, SOS50, PS90 S PS90 E and EOS50 do not allow us to observe differences between the annual phenological dynamics of the different sites as described by NDVI, NBR and NDRE1, see Tables 3 and 4. For any studied site and with any VIs, the start of greenup, indicated by SOS20, is found around DoY 90 with variations between sites of less than one week and standard deviations between a couple of days and a couple of weeks. SOS50, which indicates the midpoint of the greenup period, appears around DoY 105-115 depending on the site, with standard deviations of the same magnitude as the differences. These results are independent of the reconstruction methodology. The position of the peak of maturity, located between PS90 S and PS90 E , varies between VIs and reconstruction methodologies. Thus, for NDRE1, independently of the methodology, PS90 S PS90 E criteria predict the peak of maturity to be between DoY 140-170 and DoY 195-210. However, for NDVI and NBR, PS90 S and PS90 E criteria present differences between S-G and D-L time series. While S-G time series for both NDVI and NBR predict PS90 S around DoY 140-150 for any studied area, D-L time series predict earlier greenup of about 20 days (DoY ≈ 120-130). This behavior is inverted for PS90 E . In this case, D-L time series estimate the starting of senescence around DoY 220-240 for NDVI and 250-280 for NBR, while S-G ones estimate earlier starting of senescence at DoY 200-210 for NDVI and 230-250 for NBR. In addition, whatever the selected VIs, both PS90 S and PS90 E present high variabilities with high standard deviations of around 20 days, especially when differences between the mean values for the different sites are more important. EOS50, indicating the midpoint of senescence, also depends on the observed VI but is independent of the methodology. Thus, for NDVI EOS50 approximately appears on DoY 280-290, for NBR it is around DoY 300-310 and for NDRE1 around DoY 240-260. In the case of NDVI and NBR, the differences of PS90 S and PS90 E between S-G and D-L time series are explained because D-L does not characterize intra-annual sudden drops of VIs, which appear during the maturity period. Hence, we consider that these criteria on S-G time series provide more confident values. In the case of NDRE1, the shape of the phenology curve as described by this VI reduces the differences of PS90 S , PS90 E and SOS50 between S-G and D-L methodologies.

Intra-Annual Phenological Disturbance Events Analysis
Only the proposed methodology, with S-G filtering, allows to observe intra-annual phenological disturbance events (defined as periods with anomalous phenological dynamics), while D-L curve fitting method does not (see Figure 9). Thus for NDVI and NDRE1, Figure 9a,e shows a marked disturbance period, characterized by a sharp fall of both VIs [35,61,62], occurring from the end of July (DoY ≈ 205) to September (DoY ≈ 270). This period can be observed on every studied area with a decrease of NDVI and NDRE1 going from around −0.01 NDVI (NDRE1) per week, for London planes from Canal de Brienne, to -0.005 NDVI (NDRE1) per week, for London planes from peri-urban Canal de Midi. Another (less important) intra-annual phenological disturbance period, also characterized by a sharp decrease of VIs appears at the end of June (DoY ≈ 170) on every studied area. Table 5 shows that for NDVI and NDRE1 both, the slope and the area criteria, indicate that London planes from peri-urban Canal de Midi are the less disturbed trees (for NDVI and NDRE1 SlVI 205 = −0.005 VI/week, and DiffA = 0.010 and 0.007 respectively), while those of Canal de Brienne are more disturbed (for NDVI and NDRE1 SlVI 205 = −0.008 VI/week and −0.009 VI/week respectively, and DiffA = 0.025 and 0.033 respectively). The rest of areas present phenological disturbance in between.
For NBR, Figure 9c shows a slight phenological disturbance period occurring at the same moment of the one observed in NDVI and NDRE1 (DoY ≈ 205), with a decrease of NBR of −0.004 NBR/week, for London planes from Canal de Brienne, and of −0.002 NBR/week, for London planes from peri-urban Canal de Midi. Thus, NBR seems to be less sensitive to the anomalous phenological dynamics appearing during this period than NDVI or NDRE1. Table 5 shows that for NBR the slope criterion agrees with the results obtained for NDVI and NDRE1. Peri-urban Canal de Midi trees have the less anomalous dynamics, and newly, Canal de Brienne trees show the most anomalous dynamics. Table 5. Mean of the intra-annual disturbance criteria on S-G filtered phenological curves for NDVI, NBR and NDRE1 and for the five studied areas.

Site
SlVI   2018) with a maximum daily temperature greater than 30 • C associated to a high radiation forcing. A short heat-wave event is even recorded from 2 to 6 August. In addition, precipitations are low during the period DoY 170-287 leading to a drying of the soil that is more or less rapid and durable depending on the depth of soil layers analyzed. The water content for the near-surface layer (10-50 cm) is very sensitive to short-term precipitation variation. A rapid decrease from 0.30 to 0.12 m 3 ·m −3 is noted from DoY 170 followed by a slight increase around DoY 197 due to some rainfall events, and then the water content continues to decrease down to 0.09 m 3 ·m −3 until the end of the dry period by October. The water content of soil layer 50-100 cm decreases rapidly at the beginning of the dry period (DoY 287) and then stabilizes at around 0.12 m 3 ·m −3 . At this depth, the water content is not sensitive to short-term rainfall events. Finally the deeper soil layer presents a slower dryness following a linear decrease in soil water content from DoY 168 to DoY 280 down to 0.19 m 3 ·m −3 , and then remains unchanged for the rest of the year.  Figure 9 gives the phenological curves as described by NDVI, NBR and NDRE1 for the five studied areas together with the vertical dashed lines containing the intra-annual periods of anomalous phenological dynamics (same dashed lines as in Figure 10). These intra-annual phenological disturbance events seems to be linked with the meteorological conditions of the periods: increase in temperatures and total received irradiance and decrease in precipitations and soil water content, see Figure 10.

Methodology Hypotheses and Performances
As said in Section 3.1, the proposed method supposes that for VIs time series in urban environments, sudden falls are mainly due to noise induced by registration variability between dates. The magnitude of these sudden drops for the same co-registration error can depend on the spectral contrast between vegetation and background and on the chosen VI used to reconstruct the phenological curve. For these reasons, both the D-L curve fitting and the S-G filtering contain a weighting hypothesis to make reconstructed phenological curves converging to the upper VI envelope and thus correcting this mixed-pixel effect.
The same hypothesis has been proposed for reconstructing VIs in natural environments [45,50,51], where sudden falls in VI time series are considered to be due to atmospheric and cloud cover variability and/or bi-directional effects [51]. Hence, in these cases, the hypothesis grounds on the NIR reflectance decrease when atmospheric transmission is reduced. On the other hand, reflectances in the visible domain increase due to additive path radiance [51]. However, in urban environments, registration variability between dates appears as the main source of noise due to the important heterogeneity of impervious/natural elements, the 3D structure of cities inducing shadows and the restricted spatial extent of vegetation. This is shown in Figure 8a,b, since the amplitude of the noise depends on the pixel, this noise can not be due to atmospheric effects or the latter are strongly masked by the first.
In this work, classification is based on curve fitting and a statistical condition to avoid pixels with low VI values during the maturity period, see Section 3.1. The main advantage of this classification methodology based on curve fitting is that most of the non-vegetated pixels are excluded for next analysis, see Figure 7. Its drawback is that strongly unhealthy vegetation pixels, with phenologies not following the imposed function, may be considered as non-vegetated ones. In order to improve the classification performances, refined definitions of B v and B m depending on VIs, background and vegetation species are further recommended. In addition, to study other species, the bounds of the maturity period may be re-defined following suggestions from Section 3.1.
In the S-G filtering, the weighting hypothesis may hide or reduce possible existing intra-annual phenological disturbance periods, since the filter tends to approach the upper envelope of the time series. However, even if these intra-annual periods are difficult to distinguish from noise, it can be seen in Figure 8 that if these periods are marked enough, they are correctly characterized by the methodology (comparison between green and red curves). The anomalous dynamics of the London planes phenology during a punctual short period may be underestimated, but they are observed.
In addition, the proposed time series reconstruction methodology is semi-automatic since for both steps: curve fitting classification and Savitzky-Golay filtering, the configuration of input parameters is done by the algorithm. Thus, in the minimization needed to time-series fitting, the initial values of the parameters of Equation (1) are directly approximated by the code from the raw time series. For the S-G filtering only two parameters are needed, the size of the smoothing window (m) and the degree of the smoothing polynomial (d). As explained in Section 3.1, the code selects those parameters that lead to smaller errors. Only the range of possible m and d values should be defined, see Section 3.1. Finally, this methodology is generic and appears as applicable on any large-enough group of mono-species urban vegetation, the main limitation being the resolution of the satellite used in the study. The phenology of mono-species grouped trees located in non-urban environments where mixed pixel problems can appear, such as in orchards [75,76] or woodlands/savannas [77][78][79], might be also studied with the proposed methodology.

Vegetation Phenology Characterization
Contrary to other works that focus on NDVI time series to describe vegetation phenology [50][51][52], the proposed methodology enlarges to study other VIs such as: NBR and NDRE1. The use of normalized indices simplifies the quantification of the magnitude of intra-annual phenological disturbance periods, as well as the comparison between VI phenological curves, since the maximum and minimum possible values of the indices are fixed. In addition, analyzing several VIs allows for a more complete characterization of the phenological dynamics since: 1) the different VIs provide information on different vegetation functional and structural variables helping to have a better assessment of vegetation health [37][38][39][40][41] (Table 1) and 2) the sensibility to background, shadows and atmospheric effects varies between VIs.
Among the criteria considered as characterizing the annual phenological dynamics of urban vegetation, V I max , and Green period provide the best performances for describing the heterogeneity of the phenological dynamics of the studied sites. While V I max characterizes the annual dynamics with single-date information, that for some VIs can be interpreted as the maximal capacity of vegetation production along the year, Green period is not based on a single-date but on the whole maturity period, and is thus less dependent on the number of available images. Taking both of them into account allows to well characterize the VI behavior during the green period. Several potential factors can lead to the obtained phenological dynamics heterogeneity between sites, among them, differences in: 1) The magnitude of the influence of background (including self and cast shadow effects) on the phenology characterization, 2) the illumination and viewing conditions, 3) the environmental conditions (root space, water availability, temperature differences, etc.), 4) the age of the trees and 5) the orientation of the alignment, etc. To reduce the impact that background and shadows can introduce on the phenology characterization, the pixel classification was performed to discriminate pixels mainly covered by vegetation (where the effects of background are supposed to be lower) from mixed and non-vegetated pixels. Furthermore, the S-G filter applied on vegetation pixels has been also developed to reduce possible remaining influences of mixed pixel effects. In addition, for the five studied areas, buildings are smaller than London planes (except marginal exceptions), and so, no shadows are expected to be cast on the tree top of canopy. The illumination and viewing conditions can be discarded since illumination evolves across the year equally for the five studied sites, and for S-2 used images the viewing angles does not vary, see Section 2. Environmental differences between the studied sites are clear and explained in Section 2.1, for example, the Toulouse town hall indicated that London planes from J. Guesde avenue present health issues due to the long term impact of road works on their root system development and the quality of the soil that has been replaced for this purpose.
The other studied annual criteria: SOS20, SOS50, PS90 S , PS90 E and EOS50, based on the characterization of phenological asymmetries [80], are not able to capture the heterogeneities between the different sites. These indicators are based on the date of a given VI amplitude and are very affected by the number of available images, the maximum and minimum VI values and the existence of intra-annual periods of phenological disturbance. Hence, their interpretation when characterizing the phenological dynamics heterogeneity of urban mono-species vegetation is difficult, since differences on start of greenup and senescence are expected to be small. Nevertheless, Vrieling et al. 2018 and Zhou 2019 showed that these criteria are adapted when discriminating different species of vegetation [45,80].
For every tested annual criteria, S-G filtering and D-L fitting methodologies, both on pixels previously classified as vegetation, present similar performances. Since S-G filtering takes into account intra-annual phenological disturbance periods, its description seems to be more accurate even for annual characterizations. However, noise can remain in the S-G reconstructed time series, being a source of error.
Two criteria to quantify the influence of intra-annual phenological disturbance periods (characterized by concavity breaks [35,61,62]) on VIs reconstructed time series have been used: SlVI t i and DiffA, see Figure 5. While the slope criterion allows us to quantify the disturbance period, very locally in time, by its initial and sudden effect on VI, the area criterion allows to quantify its whole effect on vegetation phenology. The comparison between the reconstructed phenology curve and the straight line has been chosen since the straight line represents the most disturbed possibility for the phenology curve without breaking the concavity/convexity of the time series. Another option would have been to compare, during the intra-annual phenological disturbance period, the S-G reconstructed phenology area to the D-L reconstructed phenology area, see Figure 8. However, the D-L method is based on curve fitting, and thus, the estimated area during the punctual disturbance period will depend on every date of the fitted raw phenology time series. The area criterion should be used carefully, because the same disturbance period can last more or less depending on the used VI. This is not surprising since different VIs provide different information on the vegetation health status and are differently influenced by background. However, intra-annual disturbance events represent slight modifications of phenological curves and then fine spatial resolutions and band widths are important. NDVI, NBR and NDRE1 are calculated at different spatial resolutions: 10 m for NDVI and 20 m for NBR and NDRE1; and band widths: around 100-150 nm for NDVI and NBR and 15 nm for NDRE1. Thus, both low spatial resolutions and wide band widths can hide intra-annual details on phenology curves. The spatial resolution and band widths of NBR may also explain the difficulty to observe the phenological disturbance period during the first week of August 2018 on NBR time series.

Meteorological Influences on Intra-Annual Phenological Dynamics
The analysis of the evolution of meteorological conditions during 2018 shows their possible influence on the intra-annual phenological disturbance phases of the London planes that have been identified on the monitoring of NDVI, NBR and NDRE1, see Figures 9 and 10. First, air temperature and vapor pressure deficit are the two environmental parameters that govern stomatal opening and plant transpiration. The rapid change to very high heat conditions can also increase rapidly the vapor pressure deficit, that can hasten the decline of vegetation [81,82]. Thus, several works indicated the "strong statistical association between phenology and temperature" [83][84][85][86]. On the other hand, the drop in precipitation and soil water availability limits the capability of trees to uptake water by roots for transpiration [82]. For 2018, there are no measurements available at different points in the city to assess the specific meteorological conditions at the different studied sites. This would allow for a better understanding of the intra-annual phenological dynamic differences between sites, since such as Wang et al. 2016 suggested, presenting VI phenology time series versus temperature cumulative indices, and not versus DoY, allows us to better discriminate phenological dynamics [87]. Despite this lack of multiple meteorological stations, in terms of temperature, data collected over Toulouse in 2004-2005 during the CAPITOUL field campaign [88] allows us to visualize the temperature differences according to urban typologies [89]. During summertime, daytime air temperature differences of 1-2 • C were observed between city center and outskirts. These differences was larger at night up to 3-4 • C. Hence, London planes from Midi Ext are considered to be exposed to lower temperatures than London planes from the city center. On the other hand, precipitations are supposed to not strongly vary from one point to the other of the city (especially taking into account the anticyclonic weather prevailing during the most part of the summer of 2018). Thus, the precipitations measured at Meteo-France can illustrate the precipitations at the studied sites. From this consideration, the soil water content sensors located in a grassland ecosystem at Meteo-France can be interpreted as the maximal soil water contents at the Toulouse city center, where impervious surfaces reduce the soil moisture. This hypothesis is only valid for sites far from canals. Then, when soil water content at Meteo-France is low due to lack of precipitations, it is supposed to also be low in the city center far from water sources.
Finally, for Canal de Brienne, Urban Canal de Midi and Peri-urban Canal de Midi, the trees are planted near a water source in an undisturbed natural soil. Consequently, the roots can develop freely (down to deep soil) and uptake soil water. For J. Guesde Avenue and F. Verdier Avenue, trees are planted in a built-up environment where the ground is disturbed and heterogeneous. The root system has less space to expand, which limits access to water and nutrient supplies [90]. In addition, these trees are less exposed to precipitation due to the limited space of soil in which they are planted. In such conditions, trees may be more sensitive to specific disturbance events.

Conclusions
This work proposes a new methodology that enables the monitoring of urban tree vegetation phenological dynamics from Sentinel-2 VIs time series reconstruction. First, a pixel by pixel curve fitting classification was performed, based on the hypothesis that healthy vegetation phenology follows a given fixed annual behavior. This allows us to discriminate vegetation, mixed and non-vegetation pixels. Then, on vegetation pixels, a Savitzky-Golay filter is applied to reduce induced noise mainly due to registration variability between dates. This methodology does not describe the phenology curve with a fixed function. So that, with the appropriate criteria, it allows us to characterize annual and intra-annual vegetation phenological dynamics for groups of mono-species urban trees.
In this article, London planes located at five sites with different environmental conditions have been studied during a single year (2018). From one side, V I max and Green period criteria show the phenology heterogeneity existing between London planes located at the different sites. However, from another side, criteria based on the date of a given VI amplitude such as SOS20, do not appear to be adapted to illustrate this heterogeneity.
The ability of the proposed methodology to detect the existence of intra-annual phenological disturbance events, that may be linked to environmental conditions, is also shown. Two intra-annual periods of anomalous phenological dynamics were detected during summer. Their influence on the phenological curves has been quantified with two different criteria: SlVI measuring the negative VI time series slope at the beginning of the disturbance period, and DiffA measuring the lose of VI area due to the disturbance event. Furthermore, these phenological disturbance periods have been related to summer meteorological conditions. Future steps will focus on multi-year studies, for which annual criteria are expected to show differences between phenological curves from one year to another. Combined with meteorological data and field campaigns, these differences would be linked to different yearly meteorological conditions, diseases, pests, pollutants, climate change, etc. Furthermore, the proposed methodology is applicable to single-season deciduous vegetation. However, for evergreen and double-season vegetation (as it appears on tropical climates) the curve fitting classification step should be adapted [52,91]. Finally, the methodology should be tested for different satellites. The new mission VENµS with a spatial resolution of 10 m with all its bands (400-900 nm), and with a revisit time of 2 days appears as a good candidate.