Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Multi-Scale Proposal Generation for Ship Detection in SAR Images
Next Article in Special Issue
Impact of Soil Reflectance Variation Correction on Woody Cover Estimation in Kruger National Park Using MODIS Data
Previous Article in Journal
Application of Mathematical Morphological Filtering to Improve the Resolution of Chang’E-3 Lunar Penetrating Radar Data
Previous Article in Special Issue
Sensitivity of Landsat-8 OLI and TIRS Data to Foliar Properties of Early Stage Bark Beetle (Ips typographus, L.) Infestation
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 11
Issue 5
10.3390/rs11050525
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Long-Term Monitoring of Cork and Holm Oak Stands Productivity in Portugal with Landsat Imagery

by
Valentine Aubard
*,
Joana Amaral Paulo
and
João M. N. Silva
Forest Research Centre, School of Agriculture, University of Lisbon, Tapada da Ajuda, 1349-017 Lisbon, Portugal
*
Author to whom correspondence should be addressed.
Remote Sens. 2019, 11(5), 525; https://doi.org/10.3390/rs11050525
Submission received: 31 January 2019 / Revised: 20 February 2019 / Accepted: 26 February 2019 / Published: 4 March 2019
(This article belongs to the Special Issue Forest Health Monitoring)
Browse Figures
Review Reports Versions Notes

Abstract

:
Oak stands are declining in many regions of southern Europe. The goal of this paper is to assess this process and develop an effective monitoring tool for research and management. Long-term trends of the Normalized Difference Vegetation Index (NDVI) were derived and mapped at 30-m spatial resolution for all areas with a stable land cover of cork oak (Quercus suber L.) and holm oak (Quercus ilex L.) forests and agroforestry systems in mainland Portugal. NDVI, a good proxy for forest health and productivity monitoring, was obtained for the 1984–2017 period using Landsat-5 TM and Landsat-7 ETM+ imagery. TM values were adjusted to those of ETM+, after a comparison of site-specific and literature linear equations. The spatiotemporal trend analysis was performed using only July and August NDVI values, in order to minimize the spectral contribution of understory vegetation and its phenological variability, and thus, focus on the tree layer. Signs and significance of trends were obtained for six representative oak stands and the whole country with the Mann Kendall and Contextual Mann-Kendall test, respectively, and their slope was assessed with the Theil-Sen estimator. Long-term forest inventories of six study sites and NDVI time series derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) allowed validating the methodology and results with independent data. NDVI has a good relationship with cork production at the forest stand level. Pettitt tests reveal significant change-points within the trends in the period 1996–2005, when changes in drought patterns occurred. Twelve percent of the area of oak stands in Portugal presents significant decreasing trends, most of them located in mountainous regions with shallow soils. Cork oak agroforestry is the most declining oak forest type, compared to cork oak and holm oak forests. The Google Earth Engine platform proved to be a powerful tool to deal with long-term time series and for the monitoring of forests health and productivity.

Graphical Abstract

1. Introduction

Cork oak (Quercus suber L.) and holm oak (Quercus ilex L.) cover one-third of the woodland areas in continental Portugal. They have been exploited for centuries [1] mostly as an artificial extensive silvopastoral system called montado [2]. Those two species have an important value for Portugal’s economy, society, and environment, providing various forest ecosystem services as landscape and tourism hotspots, reducing fire risk and soil erosion, and increasing carbon sequestration and key-habitats for rare and endemic species [3,4,5]. Under the typical silvopastoral management system, those trees are good shelters for animals, either in winter against frost, or in summer against heat and sun exposure [3,6]. Holm oaks produce characteristic acorns that are used to feed livestock, in particular black pigs which are associated with these stands and produce premium meat, and culinary purposes. The wood is known to be extremely calorific and is traditionally used as firewood, and, before the development of fossil fuel, for charcoal production [3,7]. Cork oaks are also exploited for their acorns, but mainly for their cork, which is one of the most valuable products of non-wood forest exploitation worldwide, having remarkable mechanic properties such as low permeability to liquids and gases, resistance to rot, and high elasticity [8]. It is used mainly for bottle stoppers, but also isolation material, clothing, and other products. Portugal is the main producer of cork in the world, owning 50% of the global market, and the sector contributes to 3% of the gross domestic product of the country [9], providing many jobs, especially in rural areas.
An important decline of oaks has been noticed in the last decades in Portugal due to various possible factors such as land use changes, especially livestock grazing [10], difficult access to groundwater in summer in hilltops and shallow soils [11], climate change [12], and new pests and diseases [13,14]. This decline has been observed in several ecosystem features such as the number of trees per hectare [15,16], tree regeneration, tree canopy cover [17], or cork growth [18]. To the best of our knowledge, the present study is the first to quantify and spatialize the current decline of oak woodlands at the national scale.
The Normalized Difference Vegetation Index (NDVI) has been used for monitoring purposes [19] and is considered an adequate proxy of the global health of holm oak and cork oak stands and correspondingly of cork and acorns productivity. A strong correlation was found between wood, leaves and seed production of North American oak species and NDVI derived from the National Oceanic and Atmospheric Administration Advanced Very High Resolution Radiometer satellite (NOAA/AVHRR, 1-km of spatial resolution) [20]. NDVI from the Moderate Resolution Imaging Spectroradiometer (MODIS, 250-m resolution) was also proved to be an efficient indicator of forest biomass growth and dieback in Mediterranean holm oak forests [21]. Few studies used NDVI or other vegetation indices to estimate land cover changes and vegetation trends in the Iberian Peninsula [22,23], pre- and post-fire vegetation dynamics [24], drought response of Mediterranean evergreen oaks [25], the relationship between their phenology and precipitation [26]. Also, the Enhanced Vegetation Index (EVI) derived from Landsat imagery was used to detect productivity trends of cork and holm oaks in Serra de Grândola, Portugal [27]. This study used a 13-year time series for a 20 km2 area and a Mann-Kendall test to determine the significance of the trends. In our study, the analysis was extended to all cork and holm oak stands in Portugal, using a longer time range of 34 years. EVI could have been an interesting vegetation index for this analysis since it is less sensitive to understory spectral response. However, it has been proven to be significantly affected by differences in the sensor view geometry of Landsat satellites, as opposed to NDVI [28].
Since Landsat 4 satellite was launched in 1982, Landsat imagery allows producing long-term trends and spatial analysis at the fine resolution of 30 meters, facilitating detailed forest monitoring. A comparison of Landsat 30-m and MODIS 500-m-derived NDVI with flux tower data [29] revealed Landsat gives a vegetation phenological signal closer to the flux tower than MODIS. This approach implies to use more than one Landsat sensor. Several studies revealed disparities between the sensors, in particular between Landsat TM 5 and ETM+ 7 NDVI values that can introduce artificial long-term trends, due to differences in the bands wavelength sensibility and atmospheric conditions [28,30,31]. Those values need to be adjusted. As a consequence, our study will rely on inter-calibrated Landsat-derived NDVI time series, and will also use MODIS data to validate the adjustments performed. Phenological monitoring of both Californian Quercus douglassii savannas and pure grasslands areas with NDVI [32] showed phenological dates were consistent across spatial resolution for grasslands but varied in savannas. The conclusion of this study was that 500-m pixels give an average of the wide range of information that can be observed at a finer scale. Thus, we expect to find a difference in NDVI values between Landsat and MODIS time series, but a coherent trend direction, as long as the study-area used to compare them is large enough.
Stand crown cover in montado systems is commonly low, ranging from 10 to 50% [17]. A study estimating tree canopy cover in cork and holm oak stands using vegetation indices pointed out the necessity to use a period when the spectral contrast between the overstory and the understory is maximum to focus on overstory signal [33]. This strategy was also applied for monitoring woody vegetation in Israel [34] and for adjusting the NDVI values of two Landsat sensors [35]. Moreover, a comparison of trees (Q. suber), grass and shrub reflectance over the years [36] has proven that the cork oak NDVI response was close to steady through the year, while the herbaceous vegetation reflectance shows large variations, due to phenology, being the lowest in July and August. Our time series will be based on those months.
To follow the long-term productivity trends of cork and holm oaks for a whole country, the potentialities of the new open access platform Google Earth Engine (GEE, [37]) were exploited. GEE gives free access to thousands of geographic information system features (vector and raster) from all over the world. Its JavaScript code editor allows a user to manipulate them online, without downloading, and to save images or table results on an online storage platform. Therefore, this work will be an opportunity to evaluate to what extent GEE is adapted to the needs and purpose of image processing with the goal of forest monitoring.
The aim of this study was to analyze the spatial and temporal trends of NDVI in cork and holm oaks forests of Portugal and to identify and quantify the areas with increasing and decreasing productivity trends. Here, we extended previous research by using (i) a long time series (34 years), (ii) a higher spatial resolution (30 m), and (iii) a contextual statistical technique to take advantage of the fine spatial resolution and account for neighboring information. Additionally, we selected six case studies located in cork and holm oak permanent forest inventory plots to further explore the relationship between NDVI trends and oak woodland productivity. As a secondary aim, we developed a flexible methodology of long-term monitoring based on remotely sensed data that can be applied to other types of forests and vegetation indices.

2. Datasets and Pre-Processing

The main steps of the methodology, summarized in Figure 1, were: (i) to define a study area including all unburned cork and holm oak areas of Portugal with a constant land cover class from 1984 to 2017; (ii) to identify six representative study sites in Portugal; (iii) to produce Landsat NDVI times series, using both TM and ETM+ imagery, which required an adjustment of NDVI values from one sensor to the other; (iv) to elaborate a MODIS NDVI time series to be compared with the Landsat time series in the study sites; (v) to determine NDVI long-term trends with Mann-Kendall (MK) and Contextual Mann-Kendall (CMK) tests, and their rate of change with the Theil-Sen slope estimator (TS).

2.1. Spatial Distribution of Cork and Holm Oaks in Portugal

The spatial distribution of cork and holm oaks was derived from the Land Use and Land Cover Maps of Portugal (Cartas de Uso e Ocupação do Solo—COS, [38]), which is currently available for 1995, 2007, 2010, and 2015. The COS characterizes homogeneous areas of at least 1 hectare with a five-level hierarchical classification (193 classes for the most detailed level). The most common cork and holm oak stand cultural systems were summarized into five classes, defined following COS class descriptions [39]:
  • Cork oak forests: The tree canopy covers no less than 10% of the area (for 4-meter canopy radius, 20 trees/ha). The understory is not used for agriculture (agroforestry) or recreative activities (urban parks). This class regroups pure and mixed forests (hosting other broadleaf and coniferous species) dominated by cork oaks;
  • Holm oak forests: Forests as defined above for cork oaks, in this case dominated by holm oaks;
  • Cork oak agroforestry systems (AFS): AFS are associations of temporary crops or grasslands with permanent tree species. At least 10% of the area is covered by the forest canopy and dominated by cork oaks;
  • Holm oak agroforestry systems: AFS, as defined above, here dominated by holm oaks;
  • Agroforestry systems with cork and holm oaks: Particular case of AFS, as defined above, in which both species co-occur and none represents more than 75% of the forest canopy cover.
Since our goal was the monitoring of cork and holm oak stands and not the study of transitions from or to oak land cover classes, only the areas with a constant land cover class, i.e., areas where class transitions were not observed in the period 1995–2015, were considered in our analysis. No COS data exists to guarantee the land cover of those areas remained constant between 1984 and 1995, as well as for the period 2015–2017. New plantations were important in the 1990s due to financial programs and could produce an over-estimation of positive NDVI trends. This is anyway reflecting the global evolution of cork and holm oaks in Portugal. Most of all, the stands covered by cork and holm oaks are very stable over time, due to their type of exploitation and potential lifetime and since those protected species cannot be cut without formal permission, according to the current national legislation in Portugal. Using the available COS, it was possible to estimate that the land cover changes from 1995–2015 were mainly new plantations (35%) and conversion between a forest and an agroforestry system of the same species (31%). Many areas corresponded also to conversions from one oak species to another, always remaining in the five classes we defined (24%). The last 10% corresponded to stands replaced by urban areas, agriculture, shrublands or other woody species or to abrupt transitions (absence for 3–8 years), possibly due to fires and replantation. With those limits in mind, it was assumed the stands with a stable land cover class from 1995 to 2015 were suitable for a 1984–2017 trend analysis (Figure 2A).
Areas burned between 1984 and 2017 (Figure 2B, in grey) were used as a mask to exclude from the analysis the burned stands and again to focus on the areas that remained covered by cork or holm oaks. Burned area perimeters have been extracted from the Landsat-based Fire Atlas of Portugal [40], produced by the Forest Research Centre (School of Agriculture, University of Lisbon). This dataset provides burned area polygons over the entire mainland Portugal for each fire season between 1975 and 2013. It was produced with a semi-automatic classification of Landsat satellite imagery (Landsat TM/ETM+ era) and has a minimum mapping unit of 5 hectares. Data from 1984 to 2013 were merged with annual burned areas maps from the Institute for Nature Conservation and Forests [41] for the 2014–2017 period.
Cork and holm oak stands are not situated in the regions of Portugal most affected by fire (Figure 2), mostly located in the northern half of the country. Only the most southern region of Portugal, Algarve, as defined in the Nomenclature of Territorial Units for Statistics level 2 (NUTS II) [42], has large areas of oak stands and presents a high fire incidence. Table 1 summarizes the area covered by cork and holm oaks, before and after removing the burned areas, and the proportion of burned area removed per land cover class. Close to 13% of the stable oak stands were masked, mainly from cork and holm oak forest classes. The total area finally used in our analysis covers 1,902,360 ha.

2.2. Representative Study Sites

Six representative cork oak stands were selected to help to choose the best adjustment between TM and ETM+ NDVI values and to explore the relationship between observed NDVI trends and field productivity. Those sites contrast in their density, understory composition, cork production, debarking history, and climate, following a coast-inland gradient (Figure 2A).
Bêbeda plot (38°0’2.7″ N, 8°46’52.0″ W), situated on Arenosols [43], 6 km from the sea, is a stand of 70 years old with 13.1% canopy cover (according to the 30-m pixel Landsat Tree Cover Continuous Fields layer for 2010, which gives a percentage of the vertically projected area of vegetation of woody plants greater than 5 meters in height [44], although those values seem very low compared to field estimations (personal communication)). Contenda stand (38°4’32.6″ N, 7°6’14.2″ W) is an uneven-aged plot of 100 years old or more with natural regeneration on Leptosols, with only 8.8% canopy cover. It is the farthest study site from the sea (150 km), with a drier climate than the other stands (its mean precipitation per year is 450 mm against around 550 mm) but with similar Köppen Csa characteristics, especially a wet winter (around 55 mm of rain in December and around 75 for the other stands) and dry and hot summer months (less than 5mm of rain in average in July for all plots). Coruche (39°8’21.0″ N, 8°20’8.1″ W) is the youngest stand, about 60 years old, and has 18.8% canopy cover (estimated at 40% with unmanned aircraft vehicle image mosaic; study in preparation), hosting a cow pasture on its Podzols soil 90 km from the sea. Grândola (38°6’8.1″ N, 8°26’18.8″ W), at 30 km from the sea, is an uneven-aged plot of more than 100 years old hosting a sheep pasture on Arenosols and Podzols under only 7.3% canopy cover. Lezírias (38°48’55.7″ N, 8°51’20.1″ W) is a 70-year-old plot on Arenosols 40 km of the sea, with 23.3% canopy cover, the highest proportion of the six study sites. At last, Portel (38°15’14.2″ N, 7°38’3.1″ W) is a 100-year-old plot with 12.8% of canopy cover 100 km from the sea, hosting sheep and cows on Leptosols.
The cork productivity of those plots has been followed since the 1990s by the Forest Ecosystem Management under the Global Change (ForChange) research group of the Forest Research Centre. The area of each plot is around 1 hectare, where it fits 11 Landsat pixels. The plot center was selected in order to get the most homogeneous canopy cover possible.

2.3. Landsat and MODIS Imagery

Landsat imagery was selected because they allow a long-term time series analysis at a high spatial resolution. Landsat-5 Thematic Mapper (TM) and Landsat-7 Enhanced Thematic Mapper Plus (ETM+) data were used to create a 34-year time series from 1 January 1984 to 31 December 2017 with a 30-m spatial resolution and a 16-day revisiting period. The NDVI (Equation (1)) [45] was calculated with the Surface Reflectance (SR) images provided by GEE:
NDVI = (NIR − Red)/(NIR + Red),
where NIR and Red are TM or ETM+ Near Infrared and Red bands, respectively. Only the first category of data (Tier 1: “T1_SR”) produced by the United States Geological Survey (USGS) was used (Table 2). According to the USGS Landsat, scenes with the highest available data quality are placed into Tier 1 and are considered suitable for time-series processing analysis. A Fmask function [46] was implemented in GEE to remove clouds and shadows from Landsat images, using the quality of pixel bands (pixel_qa). Moreover, all images of ETM+, taken after 31 May 2003, contain “SLC-off gaps” due to the failure of the Scan Line Corrector (SLC) of the sensor. On GEE, the values of the missing pixels are classified as ‘NA’ values. No interpolation to fill those gaps was made.
As referred to in the Introduction section, although USGS considers all Tier 1 Landsat data radiometrically calibrated and geolocated consistently across the full collection for all the sensors, several studies revealed differences between the values of NDVI due to differences in the bands wavelength sensibility between TM and ETM+ [28,31] or atmospheric conditions [30]. Some linear relations between NDVI values from ETM+ and TM sensors have already been established for Northern Europe [47] (Equation (2)):
ETM+ = 1.0210 * TM − 0.0010,
and more recently for Canada [35] (Equation (3)):
ETM+ = 1.0370 * TM.
To compare site-specific calibration factors with those retrieved from the literature review, site-specific equations were derived from the six study sites described above, by adapting the methodology of tandem images cross-calibration [31] to the Landsat images available in Portugal. For each site, NDVI values were extracted from 1 January 1999 to 5 May 2012, the period shared by the two sensors, for each pixel of each site. Site-specific relations were assumed to be linear, with an intercept equal to zero, given that the literature equations’ intercepts were very small or null (Equations (2) and (3)). Only image values from different sensors separated by a time-gap of 8 days (minimum gap possible for Portugal) were used, backward and forward. NDVI adjustments derived from each study site were compared with those obtained from Equations (2) and (3). Since the populations do not follow a normal distribution and time series observations are auto-correlated by nature, the percentages of mean and median differences between raw and adjusted NDVI values were compared to determine the best relation to use for Portugal oak areas.
To assess the effect of using coarse spatial resolution satellite data, which would be more suitable for monitoring large study areas, the trends derived from Landsat were compared with those derived from MODIS imagery. We have selected the version 6 of the MODIS Terra Vegetation Indices (Table 2), which provides the NDVI at 250-m spatial resolution every 16 days between 2000 and 2017, already fully corrected and cloud-masked.

2.4. Trend Analysis

Under the hypothesis that the trends were monotonic, both Mann-Kendall (MK) [48,49] and Contextual Mann-Kendall tests (CMK) [50] were run in IDRISI 17.0 Selva Edition to obtain the strength and direction of the trends. The original pixel-based non-parametric test of MK, robust to outliers, determines if a time series has no monotonic trend, either uneven or no clear trend (null hypothesis). The alternative hypothesis is the existence of a significant monotonic trend (p-value < 0.05). The sign of the test’s statistics Tau, positive or negative, indicates respectively increasing or decreasing trends. The CMK test is a geostatistical object-based approach of the MK test, which takes into account the behavior of the first order neighbors to judge the significance of a pixel monotonic trend. It was assumed that the serial autocorrelation was of an order of one in the 34-year time series so obtained, most of which can be removed using the Durbin-Watson statistic [50]. The time series was thus first pre-whitened using Durbin-Watson residuals series pre-process of the IDRISI Earth Trends Modeler. A value of the trend rate of change (slope) was assessed by the Theil-Sen robust linear estimator (TS) [51,52]. The expected results were maps of monotonic trend significance obtained with both MK and CMK tests and a map showing the TS slope estimate of those trends.
Based on previous works exposed in the Introduction section regarding intra-annual variations of NDVI in woodlands, only July and August NDVI values were used for the trend analysis. This allows to minimize the spectral contribution of understory vegetation and its phenological variability, and in this way, focus on the tree layer. Moreover, Ju et al. already used exclusively those 2 months to obtain better linear relations between TM and ETM+ NDVI values [35].
A similar analysis was made for the six study sites, to compare Landsat and MODIS NDVI values and trends, validate Portugal trend map results and present an overview of stand monitoring. Trends of 1-ha average summer NDVI were estimated with the MK test and the TS estimator for MODIS and Landsat imagery. Moreover, the non-parametric Pettitt test [53] was used to approximate the year of the significant change-point within the site-study trends. Those tests were run in R software using packages stats version 3.4.2 [54], zyp version 0.10-1 [55], and trend version 1.1.0 [56].

2.5. Exploring Relationships with Field Productivity

The study sites were an opportunity to verify if the NDVI trends were coherent with the field data. The current legislation in Portugal allows stripping trees with a minimum gap of 9 years. After the not-usable first two strippings, the workable cork is extracted from 50 to 150-year-old trees or more. The cork production of the study sites had been measured for two consecutive debarking periods (a total of 17 to 22 years depending on the study site), within the 1985–2012 period. Annual rings thickness was measured after boiling cork samples [18,57]. The first idea was to compare cork production and NDVI trends. However, the cork production is higher on the year following the debarking, due to the induced hydric stress, and decreases then regularly until the next cork harvest [58]. The response of the trees depends on their own phloem tissues composition as well as on the intensity of the stripping [59], with the last varying with site conditions [60]. Only cork samples from living trees had been followed in each plot, introducing a bias in the plot’s actual health. Thus, to validate NDVI trend results with field data, it was chosen to look only at the global productive state of each stand, namely the average production of cork per tree per year. The average was calculated for each stand from samples of 23 to 36 trees and compared with the NDVI TS slope, which is judged a better indicator of the productivity state of the stands and its evolution than NDVI means, which depends mainly on the canopy cover [61], and thus, on the density of trees more than their health state.

3. Results

3.1. TM and ETM Adjustment

Between 1165 and 6445 values were used for each site, depending on the cloud cover, which was larger on the sites closer to the sea. A strong linear relationship between TM and ETM+ NDVI values was found for each site, with the coefficients of determination (R2) ranging from 0.9754 to 0.9922 (Table 3).
Table 4 presents the proportions of difference between raw or adjusted TM and ETM+ NDVI values. The adjustments give similar results. Given that using a unique equation is much more efficient in terms of processing, especially for a whole country, we selected the one presenting the minimal average error for the six study sites: Equation (3).

3.2. Country Level Cork and Holm Oak NDVI Trends

Figure 3 presents the maps of NDVI evolution for all significant trend areas following the CMK test. Figure 3A shows the sign of CMK’s Tau, with positive values characterizing increasing trends (in blue) and negative values decreasing trends (in red). The statistics of this map are summarized in Table 5: 70.65% of CMK significant areas were found increasing and 29.35% decreasing, which accounts for 28.36% and 11.78% of the total area of oak forests. These results show that more than one-tenth of the area of cork and holm oak stands in mainland Portugal have significant declining trends over the 34-year period. The most affected areas are located in the northern, coastal and central parts of the Alentejo NUTS II region, where decreasing clusters are located (Figure 3A). Conclusions regarding the Algarve region have to be moderated since most of its oak stands have suffered forest fires these last few decades and have been masked in this study (Figure 2B).
Figure 3B gives the annual rate of change of NDVI, obtained by the TS slope estimator, for the 34 years of the time series, for the areas with significant trends (Figure 3A). Rates values followed a normal distribution. For CMK significant trend pixels, 95% of the values were between −0.0065 and 0.0059. The highest negative values were found in Serra de Portel and in the south of Serra de Grândola.
CMK significant increasing and decreasing trends for each land cover class are compared in Figure 4. The proportion of significant trends by land over class varied between 35% for holm oak agroforestry systems and 47% for holm oak forests. Cork oaks agroforestry systems and forests had respectively 41% and 43%, while agroforestry systems with both species presented 36% of the significant trend areas.
According to the maps presented in Figure 4, cork oak agroforestry (C) appears to be in steeper decline in land cover (16.66% of the area of this class presents a significantly decreasing trend), while holm oak agroforestry systems (D) present less than half of this value (7.68% of the class total area). Agroforestry systems with both species (E) are in an average state between those two values (12.49%). Forests of cork oak (A) and holm oak (B) have close values, respectively 13.13% and 12.06% of their areas reveal decreasing significant trends. When several classes coexist, trend values are coherent. The main decreasing NDVI trends are found in the same locations whatever the cultural systems: for example, in the south of Grândola for cork oak forests and agroforestry systems (A and C); around Portel for holm oaks (B and D). Forests of cork and holm oaks (A and B) evolve the same way where they co-occur, both presenting downward trends around Portel and Grândola and upward trends in the east of Castelo Branco.

3.3. Study Sites Trends

Figure 5 presents the trend results at a detailed spatial scale, in windows of 3 km2 around each of the six study sites. Figure 5A shows the significant trend areas (black polygons) according to MK and CMK tests. The pixel-based MK test has a noisy result, with more isolated significant pixels that are ignored with the contextual approach (CMK), while the CMK test produces larger, more homogeneous clusters of pixels with significant trends. Those differences were expected and have already been discussed in CMK test evaluation works [50]. The CMK approach seems more representative and more suitable to apply at 30-meter spatial resolution.
It is possible to notice some large patches of pixels with the same trend (e.g., Bêbeda, Portel), which may correspond to a unique forest stand, submitted to the same pedoclimatic conditions and the same management options (age, density). On the contrary, clusters of opposite significant trends can be observed side by side, for example, around Coruche site. The detailed shape and size of those small units clearly demonstrates the interest of using 30-m spatial resolution satellites images.
The TS slope estimator (Figure 5B) gives the magnitude of the trends, which is useful to unveil the different behaviors of the stands (for example, Lezírias and Grândola show both increasing trends, but with a very different rate of change) or identify peculiar patterns (discernible roads around Bêbeda) for precise monitoring purposes.
Portel stand is not located inside a cluster on the trend maps because it does not have a constant land cover class from 1995 to 2015. A different delimitation of COS patches between 1995 and 2007, and possibly shrub encroachment, made it appear first as cork oak agroforestry systems with pasture, then as cork oak forest.

3.4. Comparison with MODIS Trends

To validate the adjustment of TM and ETM+ NDVI, Landsat trends were compared with the ones obtained with MODIS imagery. The Mann-Kendall test and the Theil-Sen estimate were applied to the July–August average NDVI time series of each study site. Landsat results are presented for the whole period (1984–2017) and for the 17 years (2000–2017), the most common period with the MODIS sensor.
Table 6 presents the Tau and p-value obtained with the Mann-Kendall test. All Landsat long-term trends (1984–2017) were found to be significant with the MK test (p-value < 0.05). Considering the last 17 years (2000-2017), significant monotonic trends appeared in Coruche, Grândola and Lezírias stands with Landsat imagery, and in Grândola and Lezírias using MODIS NDVI. The shorter duration of the time series may explain the lack of significance in some sites. Tau signs were consistent for every site between the two Landsat and MODIS time series, except for Bêbeda.
Rates of change obtained with TS linear slope estimator, expressed in NDVIunit.year−1, are presented in Table 7. Regarding the 34-year period, four stands presented increasing trends and two were significantly decreasing, Contenda and Portel. The most contrasting slope values were found in Portel (-0.0051) and Coruche plots (0.0059). Those results are in accordance with the rates of change observed in each plot map (Figure 5B). Landsat slope values are quite different for the two periods considered. For the 17-year trend, they appear either much lower (Contenda, Portel) or higher (Grândola, Lezírias) than 34-year slopes.
The non-parametric Pettitt change-point test was run for each study site of Landsat trends to detect potential change-points. Figure 6 presents the complete time series of Landsat and MODIS (respectively in blue and black) for each study site, the results of the TS linear estimator (blue lines) and Pettit test (vertical dash lines) for Landsat 34-year trends. The six study sites had a significant change-point of Landsat trend in the same decade, between 1996 and 2005. As seen in Figure 6, NDVI values follow similar annual variations for the two sensors and present coherent trends. However, the magnitude of trends, i.e., the TS slopes (Table 7), reveal some differences that may be a direct consequence of sensors disparities in spatial and spectral resolutions.

3.5. Comparison with Field Data

The mean cork production of each study site (see Section 2.5) was graphically represented as a function of the NDVI TS slope derived from the Landsat imagery (Figure 7). The two different types of data show a large agreement. Portel and Contenda stands, plots of 100 years or more, present low cork production and a clear NDVI decline, which can be the consequence of a reduction either in the number of trees or in their global health. On the contrary, Coruche plot, the youngest, stands out from the others with its high values of productivity and large NDVI increase. Lezírias, Bêbeda and Grândola sites, all on Arenosols, present similar productivities and NDVI trends despite their different proportions of canopy cover. Figure 7 suggests a linear or exponential relation between the two variables, although this study would need to be extended to more sites to strengthen these conclusions.

4. Discussion

Our study aimed at mapping and estimating the direction and magnitude of change in the health and productivity of cork and holm oaks stands in continental Portugal. The methodology was successfully applied to map 34-year NDVI trends of all stands of cork and holm oaks. The 30-m resolution of Landsat images allowed an accurate identification of spatial patterns within the oak forests, efficiently avoiding roads, water surfaces, and adjacent croplands, without its processing being over-time-consuming in GEE.
Monotonic NDVI trends were calculated for 34 consecutive years, which was never done before, especially for the whole country. The comparison of MODIS and Landsat trends for 17 years showed the importance of using long time series, which ensures that the significant trends found are not due to short-term events such as droughts or debarking, but reflect a long-term evolution of the health of the forest stands.
Summer NDVI revealed to be a good proxy of cork oak stands health and productivity for the six study sites, presenting a similar evolution to the one expected regarding field data: the results suggested a good relation between summer NDVI trends and cork growth (Figure 7). The six study sites had a significant change-point of Landsat trend in the same decade, between 1996 and 2005. According to a Portuguese drought report covering the period 1975–2006 [62], when our time series began, 1984–1989, the climate was particularly wet and was followed by a dry episode of 4 years (1990–1994). A more severe episode of drought occurred in 2002–2006. These events could explain the existence of change points in the same decade in study site trends, since cork and holm oak caliper development and photosynthesis are directly correlated to precipitation [18,63].
The trends also matched field knowledge and observations at a larger scale. A model of cork oak crown diameter applied to the Portuguese national forest inventory data from cork oak stands [17] determined a 9% decrease of cork oaks stands with a dense canopy cover (20 to 40%) between 1996 and 2006, in favor of sparse stands (less than 20% of cover). In our study, 12% of oak stands presented a significant decreasing trend.
Cork oak trends were estimated in southern Portugal for a 20 km2 area in Serra de Grândola [27]. Landsat-derived EVI time series between 2000 and 2013 were extracted for five land cover types: dense cork oak stands with and without understory, sparse cork oak stands with and without understory, and grasslands. A Mann-Kendall test was used to determine significant trends of maximum EVI (in spring, from February to May) and minimum EVI (in summer, from July to September). Trends were hardly significant due to the short duration of the series, but the MK’s Tau of each land cover type revealed to be positive (increasing trend) for summer EVI trends, which is consistent with the main NDVI trends found in our study for the Grândola study site (Figure 5A,B). They also support the only-summer-NDVI methodology applied in our work by concluding that, while the upward trend detected in summer seemed to be driven by cork oak canopy, the decreasing trend found for spring was a consequence of shrub signal, even using an index less sensitive to the background like EVI.
Significant increasing NDVI trends (28.36% of total area) were expected in most of the stands, since the tree canopy extends and gets thicker each year. In the case of old stands, natural or artificial regenerations are supposed to efficiently replace dying trees with new ones, equilibrating NDVI values. On the other hand, significantly decreasing NDVI trends (11.78% of total area) can be a consequence of the drier and drier summers observed and predicted in southern Portugal [64]: one of the first tree responses to hydric stress is a decrease in photosynthesis activity [63], causing a progressive decline in NDVI values. Downward NDVI trends can also be due to the death of main branches or of the trees themselves. Besides the more severe climatic conditions such as droughts, some pests and diseases can also cause the decline and death of evergreen oaks, as a consequence of climate change, as the exhaustive review made by De Sampaio e Paiva Camilo-Alves et al. [65] shows. They analyzed the link between oak decline and the oomycete Phytophthora cinnamomi, an oak root pathogen that causes the same symptoms as drought and can lead to the death of the tree. The authors also pointed out the importance of other factors associated with oak decline, first of which was soil depth, compaction, and texture, which are directly impacted by cultural practices. Thus, oaks decline seems a complex multifactor process. By detecting summer stress and potential deaths, NDVI allows following efficiently the health and global yield of the stand.
The five land cover classes presented different proportions of significant decreasing NDVI trends. Costa et al. [66] described through three field examples the land abandonment and fragmentation of Portugal evergreen oak plots, due to the main changes in agricultural economic priorities and labor force availability since 1970, starting in the less productive stands (sparse holm oak stands). Then, the highest percentage of downward significant trends was found for cork oak agroforestry systems and the lower for holm oak agroforestry systems. In those systems, tree density is lower than in a forest structure and the canopy cover does not allow effective water retention in summer [6,67]. A comparison of co-occurring cork and holm oaks’ response to summer drought showed that holm oaks have better adaptability to xerophytic climatic conditions, in accordance with their Portugal geographical distribution (Figure 2), possessing a more effective root system and being less vulnerable to embolism [11]; even though a study comparing tree-rings growth of Q. ilex in hotter and cooler Spanish sites suggests the drought adaption of this species has its limits and trees suffer from heat in the hotter sites [68]. Therefore, in the context of hotter and drier summers [64], physiological differences between the two species could have a long-term damaging effect on cork oaks and may benefit for now to holm oaks, especially in a low-density agroforestry system with low competitivity for water resources compared to a forest.
The Union of the Mediterranean Forest (UNAC, Portugal) classification of oak forests in Portugal based on soil and orography [69] allowed us to compare our trends with soil characteristics (Figure 8). The ‘Tejo and Sado’ river basins class is situated in arenite soils, with a more or less deep argic horizon, with reduced hydric retention and low fertility, with an altitude lower than 400 m. The ‘Alentejo’ type represents the large plains of the region, occupied by metamorphic or eruptive soils, more fertile and deep, with gentle slopes and an altitude level lower than 200 m. The ‘Serra’ (mountain in Portuguese) type has schist soils, an altitude higher than 200 m, and can present steep slopes. Figure 8 reveals that most declining stands are situated on the ‘Serra’. This can be due to the soil quality [65] accentuated by severe climate episodes [62]. The model developed by the UNAC, however, points to another cause: a lower rentability of the ‘Serra’ type, which has more production cost for lower productivity. This decline would thus be a consequence of low levels of management and investment.

5. Conclusions

A methodology was successfully applied to map 34-year trends of NDVI of cork and holm oak stands at a spatial resolution of 30 meters through a whole country. Twelve percent of the total area was found to be declining (thirty percent of the area with significant trends). The use of six study sites permitted the validation of the results, as well as the exploration of the relationship between trends and productivity, and to find change-points in trends around the end of the 20th century. This work can be reproduced, improved, and adapted to other indices, areas, and vegetal species. The resulting trend map is the first step to identify the reasons of oak decline and anticipate the consequences of the inevitable climate change to come.
This work highlights the reliability of Landsat time series and the advantage for the scientific community to have this imagery for free. It also reveals the convenience of GEE to study long time series: a fast and easy acquisition of images and pre-processing (cloud mask), data extraction (calculation of vegetation indices), time series plotting, and pixel value extraction. Nevertheless, the platform remains limited to run specific tests (CMK and TS), to extract large images, or print heavy graphics. The procedure developed in GEE for a region and a specific type of forest, in our case Portugal and the cork and holm oak forests, can be easily adapted and updated for other regions and forest types and be used as a monitoring system of long-term changes in forest ecosystems.

Author Contributions

This investigation was realized by V.A., conceptualized and supervised by J.M.N.S. and J.A.P. Data imagery were collected, processed, and analyzed by V.A. Manuscript preparation was conducted by V.A. with inputs on analysis and content review and editing by J.M.N.S. and J.A.P.

Funding

V. Aubard was funded by the European Union Erasmus+ program (Erasmus+ Stages 2018, Vetagro-Sup, France) and the International student mobility grant of the Auvergne-Rhône-Alpes region, France. J.A. Paulo and J.M.N.Silva were funded by postdoctoral grants (SFRH/BPD/96475/2013 and SFRH/BPD/109535/2015, respectively) from the Fundação para a Ciência e Tecnologia (FCT), Ministério da Ciência, Tecnologia e Ensino Superior, Portugal. The Forest Research Centre is a research unit funded by FCT (UID/AGR/00239/2013).

Acknowledgments

Authors acknowledge Companhia das Lezírias, Herdade da Contenda, and Instituto da Conservação da Natureza e das Florestas, for the access to study sites and for the use of facilities. They also thank Rita Andrews for her final English language editing.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Natividade, J.V. Subericultura; Ministério da Economia–Direcção Geral dos Serviços Florestais e Aquícolas: Lisboa, Portugal, 1950. [Google Scholar]
  2. Pereira, H.; Tomé, M. Cork oak. In Encyclopedia of Forest Sciences; Elsevier: Oxford, UK, 2004; pp. 613–620. [Google Scholar]
  3. Moreno, G.; Aviron, S.; Berg, S.; Crous-Duran, J.; Franca, A.; de Jalón, S.G.; Hartel, T.; Mirck, J.; Pantera, A.; Palma, J.H.N. Agroforestry systems of high nature and cultural value in Europe: Provision of commercial goods and other ecosystem services. In Agroforestry Systems; Springer: Berlin, Germany, 2017; pp. 1–15. [Google Scholar]
  4. Da Silva, P.M.; Aguiar, C.A.; Niemelä, J.; Sousa, J.P.; Serrano, A.R. Cork-oak woodlands as key-habitats for biodiversity conservation in Mediterranean landscapes: A case study using rove and ground beetles (Coleoptera: Staphylinidae, Carabidae). Biodivers. Conserv. 2009, 18, 605–619. [Google Scholar] [CrossRef]
  5. Modica, G.; Solano, F.; Merlino, A.; Fazio, S.D.; Barreca, F.; Laudari, L.; Fichera, C.R. Using Landsat 8 imagery in detecting cork oak (Quercus suber L.) woodlands: A case study in Calabria (Italy). J. Agric. Eng. 2016, 47, 205–215. [Google Scholar] [CrossRef]
  6. Joffre, R.; Rambal, S. Tree-grass interactions in the south-western Iberian Peninsula dehesas and montados. Sci. Chang. Planétaires Sécher. 2006, 17, 340–342. [Google Scholar]
  7. Rigueiro-Rodríguez, A.; McAdam, J.; Mosquera-Losada, M.R. Agroforestry in Europe: Current Status and Future Prospects; Springer Science & Business Media: Berlin, Germany, 2008; Volume 6, ISBN 1-4020-8271-1. [Google Scholar]
  8. Pereira, H. Cork: Biology, Production and Uses; Elsevier: Amsterdam, The Netherlands, 2007; ISBN 978-0-444-52967-1. [Google Scholar]
  9. Associação patronal do setor corticeiro (APCOR). Cork Yearbook; APCOR: Santa Maria de Lamas, Portugal, 2018. [Google Scholar]
  10. Godinho, S.; Guiomar, N.; Machado, R.; Santos, P.; Sá-Sousa, P.; Fernandes, J.P.; Neves, N.; Pinto-Correia, T. Assessment of environment, land management, and spatial variables on recent changes in montado land cover in southern Portugal. Agrofor. Syst. 2016, 90, 177–192. [Google Scholar] [CrossRef]
  11. David, T.S.; Henriques, M.O.; Kurz-Besson, C.; Nunes, J.; Valente, F.; Vaz, M.; Pereira, J.S.; Siegwolf, R.; Chaves, M.M.; Gazarini, L.C.; et al. Water-use strategies in two co-occurring Mediterranean evergreen oaks: Surviving the summer drought. Tree Physiol. 2007, 27, 793–803. [Google Scholar] [CrossRef] [PubMed]
  12. Kim, H.N.; Jin, H.Y.; Kwak, M.J.; Khaine, I.; You, H.N.; Lee, T.Y.; Ahn, T.H.; Woo, S.Y. Why does Quercus suber species decline in Mediterranean areas? J. Asia-Pac. Biodivers. 2017, 10, 337–341. [Google Scholar] [CrossRef]
  13. Brasier, C.M.; Robredo, F.; Ferraz, J.F.P. Evidence for Phytophthora cinnamomi involvement in Iberian oak decline. Plant Pathol. 1993, 42, 140–145. [Google Scholar] [CrossRef]
  14. Romero, M.A.; Sánchez, J.E.; Jiménez, J.J.; Belbahri, L.; Trapero, A.; Lefort, F.; Sánchez, M.E. New Pythium Taxa Causing Root Rot on Mediterranean Quercus Species in South-west Spain and Portugal. J. Phytopathol. 2007, 155, 289–295. [Google Scholar] [CrossRef]
  15. Costa, A.; Pereira, H.; Madeira, M. Landscape dynamics in endangered cork oak woodlands in Southwestern Portugal (1958–2005). Agrofor. Syst. 2009, 77, 83. [Google Scholar] [CrossRef]
  16. Silva, J.S. Os Montados; Fundação Luso-Americana: Público, Lisboa, 2007; ISBN 978-989-619-106-1. [Google Scholar]
  17. Paulo, J.; Faias, S.; Ventura-Giroux, C.; Tomé, M. Estimation of stand crown cover using a generalized crown diameter model: Application for the analysis of Portuguese cork oak stands stocking evolution. IForest Biogeosci. For. 2016, 9, 437–444. [Google Scholar] [CrossRef]
  18. Paulo, J.A.; Pereira, H.; Tomé, M. Analysis of variables influencing tree cork caliper in two consecutive cork extractions using cork growth index modelling. Agrofor. Syst. 2017, 91, 221–237. [Google Scholar] [CrossRef]
  19. Pettorelli, N. Chapter 5: NDVI and environmental monitoring & Chapter 6: NDVI and plant ecology. In The Normalized Difference Vegetation Index; Oxford University Press: Oxford, UK, 2013; pp. 56–80. [Google Scholar]
  20. Wang, J.; Rich, P.M.; Price, K.P.; Kettle, W.D. Relations between NDVI and tree productivity in the central Great Plains. Int. J. Remote Sens. 2004, 25, 3127–3138. [Google Scholar] [CrossRef]
  21. Ogaya, R.; Barbeta, A.; Başnou, C.; Peñuelas, J. Satellite data as indicators of tree biomass growth and forest dieback in a Mediterranean holm oak forest. Ann. For. Sci. 2015, 72, 135–144. [Google Scholar] [CrossRef]
  22. Julien, Y.; Sobrino, J.A.; Mattar, C.; Ruescas, A.B.; Jiménez-Muñoz, J.C.; Sòria, G.; Hidalgo, V.; Atitar, M.; Franch, B.; Cuenca, J. Temporal analysis of normalized difference vegetation index (NDVI) and land surface temperature (LST) parameters to detect changes in the Iberian land cover between 1981 and 2001. Int. J. Remote Sens. 2011, 32, 2057–2068. [Google Scholar] [CrossRef]
  23. Novillo, C.J.; Arrogante-Funes, P.; Romero-Calcerrada, R. Recent NDVI Trends in Mainland Spain: Land-Cover and Phytoclimatic-Type Implications. ISPRS Int. J. Geo-Inf. 2019, 8, 43. [Google Scholar] [CrossRef]
  24. Röder, A.; Hill, J.; Duguy, B.; Alloza, J.A.; Vallejo, R. Using long time series of Landsat data to monitor fire events and post-fire dynamics and identify driving factors. A case study in the Ayora region (eastern Spain). Remote Sens. Environ. 2008, 112, 259–273. [Google Scholar] [CrossRef]
  25. Reichstein, M.; Ciais, P.; Papale, D.; Valentini, R.; Running, S.; Viovy, N.; Cramer, W.; Granier, A.; Ogée, J.; Allard, V.; et al. Reduction of ecosystem productivity and respiration during the European summer 2003 climate anomaly: A joint flux tower, remote sensing and modelling analysis. Glob. Chang. Biol. 2007, 13, 634–651. [Google Scholar] [CrossRef]
  26. Ramos, A.; Pereira, M.J.; Soares, A.; do Rosário, L.; Matos, P.; Nunes, A.; Branquinho, C.; Pinho, P. Seasonal patterns of Mediterranean evergreen woodlands (Montado) are explained by long-term precipitation. Agric. For. Meteorol. 2015, 202, 44–50. [Google Scholar] [CrossRef]
  27. Santos, M.J.; Baumann, M.; Esgalhado, C. Drivers of Productivity Trends in Cork Oak Woodlands over the Last 15 Years. Remote Sens. 2016, 8, 486. [Google Scholar] [CrossRef]
  28. Sulla-Menashe, D.; Friedl, M.A.; Woodcock, C.E. Sources of bias and variability in long-term Landsat time series over Canadian boreal forests. Remote Sens. Environ. 2016, 177, 206–219. [Google Scholar] [CrossRef]
  29. Kovalskyy, V.; Roy, D.P.; Zhang, X.Y.; Ju, J. The suitability of multi-temporal web-enabled Landsat data NDVI for phenological monitoring—A comparison with flux tower and MODIS NDVI. Remote Sens. Lett. 2012, 3, 325–334. [Google Scholar] [CrossRef]
  30. Van Leeuwen, W.J.D.; Orr, B.J.; Marsh, S.E.; Herrmann, S.M. Multi-sensor NDVI data continuity: Uncertainties and implications for vegetation monitoring applications. Remote Sens. Environ. 2006, 100, 67–81. [Google Scholar] [CrossRef]
  31. Teillet, P.; Barker, J.; Markham, B.; Irish, R.; Fedosejevs, G.; Storey, J. Radiometric cross-calibration of the Landsat-7 ETM+ and Landsat-5 TM sensors based on tandem data sets. Remote Sens. Environ. 2001, 78, 39–54. [Google Scholar] [CrossRef]
  32. Liu, Y.; Hill, M.J.; Zhang, X.; Wang, Z.; Richardson, A.D.; Hufkens, K.; Filippa, G.; Baldocchi, D.D.; Ma, S.; Verfaillie, J.; et al. Using data from Landsat, MODIS, VIIRS and PhenoCams to monitor the phenology of California oak/grass savanna and open grassland across spatial scales. Agric. For. Meteorol. 2017, 237–238, 311–325. [Google Scholar] [CrossRef]
  33. Carreiras, J.M.B.; Pereira, J.M.C.; Pereira, J.S. Estimation of tree canopy cover in evergreen oak woodlands using remote sensing. For. Ecol. Manag. 2006, 223, 45–53. [Google Scholar] [CrossRef]
  34. Helman, D.; Lensky, I.M.; Tessler, N.; Osem, Y. A Phenology-Based Method for Monitoring Woody and Herbaceous Vegetation in Mediterranean Forests from NDVI Time Series. Remote Sens. 2015, 7, 12314–12335. [Google Scholar] [CrossRef]
  35. Ju, J.; Masek, J.G. The vegetation greenness trend in Canada and US Alaska from 1984–2012 Landsat data. Remote Sens. Environ. 2016, 176, 1–16. [Google Scholar] [CrossRef]
  36. Cerasoli, S.; e Silva, F.C.; Silva, J.M. Temporal dynamics of spectral bioindicators evidence biological and ecological differences among functional types in a cork oak open woodland. Int. J. Biometeorol. 2016, 60, 813–825. [Google Scholar] [CrossRef] [PubMed]
  37. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 2017, 202, 18–27. [Google Scholar] [CrossRef]
  38. DGTerritório-Carta de Uso e Ocupação do Solo de Portugal Continental (COS). Available online: http://www.dgterritorio.pt/dados_abertos/cos/ (accessed on 11 November 2018).
  39. Caetano, M. Especificações Técnicas da Carta de Uso e Ocupação do Solo de Portugal Continental Para 1995, 2007, 2010 e 2015; Direção-Geral do Território: Lisboa, Portugal, 2018; p. 103. [Google Scholar]
  40. Oliveira, S.L.J.; Pereira, J.M.C.; Carreiras, J.M.B. Fire frequency analysis in Portugal (1975–2005), using Landsat-based burnt area maps. Int. J. Wildland Fire 2012, 21, 48–60. [Google Scholar] [CrossRef]
  41. Informação Geográfica—ICNF. Available online: http://www2.icnf.pt/portal/florestas/dfci/inc/info-geo (accessed on 9 December 2018).
  42. Regulation (EC) No 1059/2003 of the European Parliament and of the Council of 26 May 2003 on the establishment of a common classification of territorial units for statistics (NUTS). Available online: http://data.europa.eu/eli/reg/2003/1059/2018-01-18 (accessed on 11 January 2019).
  43. Food and Agriculture Organization of the United Nations. World Reference Base for Soil Resources 2014: International Soil Classification System for Naming Soils and Creating Legends for Soil Map; FAO: Rome, Italy, 2014; ISBN 978-92-5-108369-7. [Google Scholar]
  44. Sexton, J.O.; Song, X.-P.; Feng, M.; Noojipady, P.; Anand, A.; Huang, C.; Kim, D.-H.; Collins, K.M.; Channan, S.; DiMiceli, C.; et al. Global, 30-m resolution continuous fields of tree cover: Landsat-based rescaling of MODIS vegetation continuous fields with lidar-based estimates of error. Int. J. Digit. Earth 2013, 6, 427–448. [Google Scholar] [CrossRef]
  45. Rouse, J., Jr.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring Vegetation Systems in the Great Plains with ERTS; NTRS: Chicago, IL, USA, 1974. [Google Scholar]
  46. Zhu, Z.; Woodcock, C.E. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sens. Environ. 2012, 118, 83–94. [Google Scholar] [CrossRef]
  47. Steven, M.D.; Malthus, T.J.; Baret, F.; Xu, H.; Chopping, M.J. Intercalibration of vegetation indices from different sensor systems. Remote Sens. Environ. 2003, 88, 412–422. [Google Scholar] [CrossRef]
  48. Mann, H.B. Nonparametric tests against trend. Econ. J. Econ. Soc. 1945, 13, 245–259. [Google Scholar] [CrossRef]
  49. Kendall, M.G. Rank Correlation Methods; American Psychological Association: Washington, DC, USA, 1975; Volume 35. [Google Scholar]
  50. Neeti, N.; Eastman, J.R. A Contextual Mann-Kendall Approach for the Assessment of Trend Significance in Image Time Series. Trans. GIS 2011, 15, 599–611. [Google Scholar] [CrossRef]
  51. Theil, H. A rank-invariant method of linear and polynomial regression analysis. In Henri Theil’s Contributions to Economics and Econometrics; Springer: Berlin, Germany, 1992; pp. 345–381. [Google Scholar]
  52. Sen, P.K. Estimates of the regression coefficient based on Kendall’s tau. J. Am. Stat. Assoc. 1968, 63, 1379–1389. [Google Scholar] [CrossRef]
  53. Pettitt, A.N. A non-parametric approach to the change-point problem. Appl. Stat. 1979, 28, 126–135. [Google Scholar] [CrossRef]
  54. Team, R.C. R: A language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2017; ISBN 3-900051-07-0. Available online: https://www.r-project.org/ (accessed on 1 March 2019).
  55. Bronaugh, D.; Werner, A. Zyp: Zhang+ Yue-Pilon Trends Package (Pacific Climate Impacts Consortium); R Package Version 0.10-1; R Foundation for Statistical Computing: Vienna, Austria, 2013. [Google Scholar]
  56. Pohlert, T. Trend: Non-Parametric Trend Tests and Change-Point Detection; R Package Version 1.1.0; R Foundation for Statistical Computing: Vienna, Austria, 2018. [Google Scholar]
  57. Faias, S.P.; Paulo, J.A.; Palma, J.H.N.; Tomé, M. Understory effect on tree and cork growth in cork oak woodlands. For. Syst. 2018, 27, 2. [Google Scholar] [CrossRef]
  58. Costa, A.; Pereira, H.; Oliveira, A. Variability of radial growth in cork oak adult trees under cork production. For. Ecol. Manag. 2003, 175, 239–246. [Google Scholar] [CrossRef]
  59. Correia, O.A.; Oliveira, G.; Martins-Loução, M.A.; Catarino, F.M. Effects of bark-stripping on the water relations of Quercus suber L. Scientia gerundensis. 1992, 18, 195–204. [Google Scholar]
  60. Paulo, J.A.; Tomé, M. Does debarking intensity during the first cork extraction affect future cork thickness? Ann. For. Sci. 2017, 74, 66. [Google Scholar] [CrossRef]
  61. Häusler, M.; Silva, J.M.; Cerasoli, S.; López-Saldaña, G.; Pereira, J.M. Modelling spectral reflectance of open cork oak woodland: A simulation analysis of the effects of vegetation structure and background. Int. J. Remote Sens. 2016, 37, 492–515. [Google Scholar] [CrossRef]
  62. Vivas, E.; Maia, R. Characterization of Historical Drought situations in mainland Portugal—The importance of indicators use. In 2.as Jornadas de Hidráulica, Recursos Hídricos e Ambiente; Faculdade de Engenharia da Universidade do Porto: Porto, Portugal, 2007; pp. 51–61. ISBN 978-989-95557-1-6. [Google Scholar]
  63. Vaz, M.; Pereira, J.S.; Gazarini, L.C.; David, T.S.; David, J.S.; Rodrigues, A.; Maroco, J.; Chaves, M.M. Drought-induced photosynthetic inhibition and autumn recovery in two Mediterranean oak species (Quercus ilex and Quercus suber). Tree Physiol. 2010, 30, 946–956. [Google Scholar] [CrossRef] [PubMed]
  64. Pal, J.S.; Giorgi, F.; Bi, X. Consistency of recent European summer precipitation trends and extremes with future regional climate projections. Geophys. Res. Lett. 2004, 31. [Google Scholar] [CrossRef]
  65. De Sampaio e Paiva Camilo-Alves, C.; da Clara, M.I.E.; de Almeida Ribeiro, N.M.C. Decline of Mediterranean oak trees and its association with Phytophthora cinnamomi: A review. Eur. J. For. Res. 2013, 132, 411–432. [Google Scholar] [CrossRef]
  66. Costa, A.; Madeira, M.; Lima Santos, J.; Oliveira, Â. Change and dynamics in Mediterranean evergreen oak woodlands landscapes of Southwestern Iberian Peninsula. Landsc. Urban Plan. 2011, 102, 164–176. [Google Scholar] [CrossRef]
  67. Gouveia, A.C.; Freitas, H. Intraspecific competition and water use efficiency in Quercus suber: Evidence of an optimum tree density? Trees 2008, 22, 521. [Google Scholar] [CrossRef]
  68. Gea-Izquierdo, G.; Cherubini, P.; Cañellas, I. Tree-rings reflect the impact of climate change on Quercus ilex L. along a temperature gradient in Spain over the last 100years. For. Ecol. Manag. 2011, 262, 1807–1816. [Google Scholar] [CrossRef]
  69. União da Floresta Mediterrânica (UNAC). Capitulo 4: Caraterização dos modelos de produção de cortiça. In Estudo de Determinação dos Pontos de Valor Acrescentado da Fileira da Cortiça; UNAC: Lisbon, Portugal, 2013; pp. 17–25. [Google Scholar]
Remotesensing 11 00525 g001 550
Figure 1. Flowchart of the methodology, allowing to build long-term NDVI time series (top, in Google Earth Engine), to compare MK and CMK results at Portugal scale (bottom left, in Idrisi) as well as Landsat and MODIS NDVI results for the six study sites (bottom right, in R).
Figure 1. Flowchart of the methodology, allowing to build long-term NDVI time series (top, in Google Earth Engine), to compare MK and CMK results at Portugal scale (bottom left, in Idrisi) as well as Landsat and MODIS NDVI results for the six study sites (bottom right, in R).
Remotesensing 11 00525 g001
Remotesensing 11 00525 g002 550
Figure 2. Study area, consisting of cork and holm oak areas with a constant land cover class in the period 1995–2015, according to the Land Use and Land Cover Maps of Portugal (A); burned areas in mainland Portugal (from 1984 to 2017), which were removed from the analysis, with the locations of the six study sites (B).
Figure 2. Study area, consisting of cork and holm oak areas with a constant land cover class in the period 1995–2015, according to the Land Use and Land Cover Maps of Portugal (A); burned areas in mainland Portugal (from 1984 to 2017), which were removed from the analysis, with the locations of the six study sites (B).
Remotesensing 11 00525 g002
Remotesensing 11 00525 g003 550
Figure 3. NDVI trends for cork and holm oak areas in Portugal for the period 1984–2017, with (A) significant decreasing (in red) and increasing trends (in blue) based on the Contextual Mann-Kendall test; and (B) the NDVI rate of change given by TS estimator.
Figure 3. NDVI trends for cork and holm oak areas in Portugal for the period 1984–2017, with (A) significant decreasing (in red) and increasing trends (in blue) based on the Contextual Mann-Kendall test; and (B) the NDVI rate of change given by TS estimator.
Remotesensing 11 00525 g003
Remotesensing 11 00525 g004 550
Figure 4. Significant NDVI trends for cork and holm oak areas in Portugal (1984–2017), based on the Contextual Mann-Kendall test, per land cover class: (A) cork oak forests; (B) holm oak forests; (C) cork oak agroforestry systems; (D) holm oak agroforestry systems; and (E) agroforestry systems with cork oak and holm oak.
Figure 4. Significant NDVI trends for cork and holm oak areas in Portugal (1984–2017), based on the Contextual Mann-Kendall test, per land cover class: (A) cork oak forests; (B) holm oak forests; (C) cork oak agroforestry systems; (D) holm oak agroforestry systems; and (E) agroforestry systems with cork oak and holm oak.
Remotesensing 11 00525 g004
Remotesensing 11 00525 g005a 550Remotesensing 11 00525 g005b 550
Figure 5. Details around the six study sites: (A) comparison of significant NDVI trend areas (black polygons) determined by MK (left) and CMK tests (right), Tau values allow to visualize increasing (in blue) and decreasing (in red) trends; (B) rate of change of NDVI values (TS slope in NDVIunit.year−1), from red (lowest) to blue color (highest), with the limits of CMK significant trend areas overlaid (in black).
Figure 5. Details around the six study sites: (A) comparison of significant NDVI trend areas (black polygons) determined by MK (left) and CMK tests (right), Tau values allow to visualize increasing (in blue) and decreasing (in red) trends; (B) rate of change of NDVI values (TS slope in NDVIunit.year−1), from red (lowest) to blue color (highest), with the limits of CMK significant trend areas overlaid (in black).
Remotesensing 11 00525 g005aRemotesensing 11 00525 g005b
Remotesensing 11 00525 g006 550
Figure 6. July–August times series per study site for Landsat 34-year trends (blue) and MODIS (black)-derived NDVI values, with the linear estimation of Theil-Sen and the Pettitt’s change point of Landsat significant trends.
Figure 6. July–August times series per study site for Landsat 34-year trends (blue) and MODIS (black)-derived NDVI values, with the linear estimation of Theil-Sen and the Pettitt’s change point of Landsat significant trends.
Remotesensing 11 00525 g006
Remotesensing 11 00525 g007 550
Figure 7. Twenty-year mean cork production in millimeters per tree per year as a function of 34-year Landsat NDVI trends for the six study sites.
Figure 7. Twenty-year mean cork production in millimeters per tree per year as a function of 34-year Landsat NDVI trends for the six study sites.
Remotesensing 11 00525 g007
Remotesensing 11 00525 g008 550
Figure 8. Comparison of the location of three types of oak forest and the TS rate of change (NDVIunit.year−1) found for all CMK significant trends.
Figure 8. Comparison of the location of three types of oak forest and the TS rate of change (NDVIunit.year−1) found for all CMK significant trends.
Remotesensing 11 00525 g008
Table
Table 1. Area covered by each land cover class concerning cork and holm oaks in Portugal before and after fire masking and the proportion of the area removed.
Table 1. Area covered by each land cover class concerning cork and holm oaks in Portugal before and after fire masking and the proportion of the area removed.
Land Cover ClassInitial Area (ha)Fire-Masked Area (ha)Proportion of Area Removed (%)
Cork oak forests844,153663,64421.4
Holm oak forests264,896225,87914.7
Cork oak agroforestry systems292,334273,5986.4
Holm oak agroforestry systems614,577586,1244.6
Agroforestry systems with cork oak and holm oak161,158153,1145.0
Total2,177,1171,902,36012.6
Table
Table 2. Image collections used, available in Google Earth Engine (GEE).
Table 2. Image collections used, available in Google Earth Engine (GEE).
SensorRange of DatesTime ScaleSpatial ResolutionGEE Image Collection ID
Landsat-5 TM1984-01-01–2012-05-0516-day cycle30 × 30 mLANDSAT/LT05/C01/T1_SR
Landsat-7 ETM+1999-01-01–2017-12-3116-day cycle30 × 30 mLANDSAT/LE07/C01/T1_SR
MODIS Terra V6
Vegetation Indices		2000-02-18–2017-12-3116-day cycle250 × 250 mMODIS/006/MOD13Q1
Table
Table 3. Site-specific linear relations between ETM+ and TM NDVI values.
Table 3. Site-specific linear relations between ETM+ and TM NDVI values.
LocationEquationR2
BêbedaETM+ = 1.0205 * TM0.9866
ContendaETM+ = 1.0101 * TM0.9754
CorucheETM+ = 1.0124 * TM0.9922
GrândolaETM+ = 1.0109 * TM0.9821
LezíriasETM+ = 1.0039 * TM0.9878
PortelETM+ = 1.0345 * TM0.9772
Table
Table 4. Comparison of TM and ETM+ NDVI values for each study site: proportion of mean and median errors using raw data and three different adjustments.
Table 4. Comparison of TM and ETM+ NDVI values for each study site: proportion of mean and median errors using raw data and three different adjustments.
SiteTypeRaw DataEquation (2)Equation (3)Site Specific
BêbedaMean error (%)−5.275−3.492−1.770−3.329
ContendaMean error (%)−8.785−7.085−5.410−7.868
CorucheMean error (%)−3.499−1.6440.072−2.303
GrândolaMean error (%)−9.750−8.074−6.411−8.764
LezíriasMean error (%)−1.7070.1771.930−1.324
PortelMean error (%)−7.582−5.880−4.162−4.391
Average mean error (%)−6.099−4.392−3.292−4.663
BêbedaMedian error (%)−5.010−3.224−1.496−3.060
ContendaMedian error (%)−13.134−11.540−9.920−12.260
CorucheMedian error (%)−3.293−1.4350.285−2.094
GrândolaMedian error (%)−10.412−8.772−7.097−9.433
LezíriasMedian error (%)−1.8150.0661.818−1.432
PortelMedian error (%)−6.667−4.963−3.213−3.444
Average median error (%)−6.722−5.000−3.971−5.287
Table
Table 5. Area and proportion of the cork and holm oak total area found for NDVI trends in Portugal, obtained with the Contextual Mann-Kendall test.
Table 5. Area and proportion of the cork and holm oak total area found for NDVI trends in Portugal, obtained with the Contextual Mann-Kendall test.
Area (ha)Proportion of Total Area (%)
Significant trends763,57340.14
of which:
-
increasing trends
539,45128.36
-
decreasing trends
224,10911.78
Not significant trends1,138,78759.86
Total area1,902,360100
Table
Table 6. Comparison of Landsat and MODIS Mann-Kendall test results for each study site.
Table 6. Comparison of Landsat and MODIS Mann-Kendall test results for each study site.
Study SiteLandsat (1984–2017)Landsat (2000–2017)MODIS (2000–2017)
Taup-ValueTaup-ValueTaup-Value
Bêbeda0.45818.64 × 10−5−0.00651.00000.03270.8814
Contenda−0.48313.14 × 10−5−0.22880.2008−0.05250.7617
Coruche0.65422.96 × 10−90.34640.04770.12420.5009
Grândola0.47594.23 × 10−50.47710.00510.59480.0003
Lezírias0.60544.99 × 10−70.55080.00150.44590.0100
Portel−0.63649.83 × 10−9−0.20260.2599−0.22880.2008
Table
Table 7. Comparison of Landsat and MODIS Theil-Sen linear slope estimate for each study site.
Table 7. Comparison of Landsat and MODIS Theil-Sen linear slope estimate for each study site.
Study SiteLandsat (1984–2017)Landsat (2000–2017)MODIS (2000–2017)
Bêbeda0.0034 *−0.000250.0002
Contenda−0.0042 *−0.00199−0.0005
Coruche0.0059 *0.00339 *0.001
Grândola0.002 *0.00284 *0.0043 *
Lezírias0.0043 *0.00478 *0.0034 *
Portel−0.0051 *−0.00148−0.0014
Significant trends according to Mann-Kendall test (p-value < 0.05) are marked with an asterisk (*).

Share and Cite

MDPI and ACS Style

Aubard, V.; Paulo, J.A.; Silva, J.M.N. Long-Term Monitoring of Cork and Holm Oak Stands Productivity in Portugal with Landsat Imagery. Remote Sens. 2019, 11, 525. https://doi.org/10.3390/rs11050525

AMA Style

Aubard V, Paulo JA, Silva JMN. Long-Term Monitoring of Cork and Holm Oak Stands Productivity in Portugal with Landsat Imagery. Remote Sensing. 2019; 11(5):525. https://doi.org/10.3390/rs11050525

Chicago/Turabian Style

Aubard, Valentine, Joana Amaral Paulo, and João M. N. Silva. 2019. "Long-Term Monitoring of Cork and Holm Oak Stands Productivity in Portugal with Landsat Imagery" Remote Sensing 11, no. 5: 525. https://doi.org/10.3390/rs11050525

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop