Open Access is an initiative that aims to make scientific research freely available to all. To date our community has made over 100 million downloads. It’s based on principles of collaboration, unobstructed discovery, and, most importantly, scientific progression. As PhD students, we found it difficult to access the research we needed, so we decided to create a new Open Access publisher that levels the playing field for scientists across the world. How? By making research easy to access, and puts the academic needs of the researchers before the business interests of publishers.
We are a community of more than 103,000 authors and editors from 3,291 institutions spanning 160 countries, including Nobel Prize winners and some of the world’s most-cited researchers. Publishing on IntechOpen allows authors to earn citations and find new collaborators, meaning more people see your work not only from your own field of study, but from other related fields too.
To purchase hard copies of this book, please contact the representative in India:
CBS Publishers & Distributors Pvt. Ltd.
www.cbspd.com
|
customercare@cbspd.com
The study employs a Fourier transform analysis approach to assess the land-cover changes in a mountainous Mediterranean protected area using multi-temporal satellite images. Harmonic analysis was applied to a time series of Landsat satellite images acquired from 1984 to 2008 to extract information about land cover status with a vegetation spectral index, the Normalized Difference Vegetation Index (NDVI). Ancillary cartographic information depicting land cover classes and the enlargement of the protected area over time (i.e., maps showing the original delineation in 1995 and subsequent enlargement in 2007) were employed as additional factors to understand vegetation-cover changes. Significant differences in the NDVI and harmonic components values were observed with respect to both factors. The application of the Fourier transform was particularly successful to extract subtle information. The harmonic analysis of the NDVI time series revealed valuable information about the evolution of the landscape. The initially protected area (northern sector) seems more affected by human activities than the southern sector (enlarged area in 2007) as revealed by the analysis of the first harmonic component that was closely related with vegetation coverage. Rural abandonment is a major driver of land-cover changes in the study area.
Department of Agrochemistry and Environment, University Miguel Hernández of Elche. Av/ Universidad s/n, Edificio Alcudia. Elche (Alicante), Spain
Encarni I. Hernández
Department of Agrochemistry and Environment, University Miguel Hernández of Elche. Av/ Universidad s/n, Edificio Alcudia. Elche (Alicante), Spain
Jose Navarro-Pedreño
Department of Agrochemistry and Environment, University Miguel Hernández of Elche. Av/ Universidad s/n, Edificio Alcudia. Elche (Alicante), Spain
Ignacio Gómez
Department of Agrochemistry and Environment, University Miguel Hernández of Elche. Av/ Universidad s/n, Edificio Alcudia. Elche (Alicante), Spain
Magaly Koch
Center for Remote Sensing, Boston University, Boston, Massachusetts, USA
*Address all correspondence to: imelendez@umh.es
1. Introduction
Land cover (i.e., biophysical attributes of the Earth’s surface) changes play an important role in the global environmental changes [1,2]. Land-cover changes may result from human-induced land-use changes or natural processes such as climatic variability and natural disturbances [3]. Studying land-cover change is a challenge because it represents a global dynamic phenomenon affecting the Earth’s large surface areas and continuously evolves throughout time. Fortunately, the status and rate of land-cover changes may be estimated and monitored through the analysis of remotely sensed data [4]. Remote sensing allows the study of the role of terrestrial vegetation in large-scale global processes (e.g., the carbon cycle) [5]. The scientific community often relies on data gathered by satellite sensors because of their synoptic view, worldwide coverage, and repeated temporal sampling capability of the Earth’s surface enabling monitoring of vegetation dynamics from regional to global scales [6,7].
Optical multispectral remote sensing sensors record the reflected portion of electromagnetic radiation-illuminating land covers. They record reflected light as digital counts that can be later transformed into comparable physical measurements such as surface reflectance, which facilitates the comparison with spectrometer measurements of surface covers. Many data transformation techniques (i.e., mathematical operations to reexpress the information content of multispectral and hyperspectral data) exist to perform spectral comparison across sensor platforms, which greatly facilitate the interpretation of land cover types [8]. In this context, vegetation indices have been used to derive measures that correlate with surface biophysical properties [9] and as such facilitate the analysis of large amounts of satellite data for spatial- and temporal-scale analyses [10]. Vegetation indices are directly related to plant vigor, density, and growth conditions and can be used to detect environmental conditions and human activities [4,11]. The most frequently used vegetation index is the Normalized Difference Vegetation Index (NDVI) as described by Rousse et al. [12], which is employed in this study.
Vegetation indices are frequently employed for the characterization of vegetation dynamics and landscape changes using satellite time-series data. They are fundamental data for monitoring vegetation phenology or land-cover changes associated with events such as fire, drought, land-use conversion, and climate fluctuation [13,14]. However, using vegetation-index time series for the analyses of vegetation changes involves a highly complex task and requires the application of specific techniques and methodologies that have been developed over time [11]. The application of Fourier Transform (FT) technique (harmonic analysis) facilitates the extraction of valuable and interpretable characteristics from the time series, which are usually distorted by atmospheric noise, sensor instability, and/or orbit deviations [15]. Harmonic analysis of vegetation-index time series can be used to analyze changes in land covers and vegetation status by examining the altered values of amplitude, phase, or the additive term over a period of years [11,16].
The objective of this study was the evaluation of landscape dynamics by the Fourier Transform technique for the analysis of land-cover changes using a satellite image time series. Harmonic analysis was applied to a time series of Landsat satellite images acquired from 1984 to 2008 to extract information about landscape dynamics using a vegetation index. Landscape dynamics were assessed by combining harmonic analysis results with ancillary cartographic information about the types of land cover classes of the protected area.
The study area is located in the Rodeno pine forest protected reserve (Paisaje Protegido de los Pinares de Rodeno or PPPR in Spanish) of the southern Iberian System in the Sierra de Albarracín shire (Teruel province, Spain). The study area is positioned around 40.33 °N and 1.36°W. The PPPR was established in 1995 covering an area of 3355.34 ha (Figure 1). In 2007, the protected landscape was enlarged to about 6829.05 ha. The PPPR occupies part of the municipalities of Albarracín, Gea de Albarracín and Bezas, with a current population (2015) of 1049, 387, and 67 inhabitants, respectively.
Figure 1.
Location map of the study area, the Rodeno pine forest protected reserve (PPPR). Original and enlarged perimeters of the protected area are shown.
It is a mountainous area with maximum altitudes ranging from 1063 to 1601 m and a mean altitude of 1315 m. The climate is characterized by significant seasonal temperature variations between summer (up to 35°C) and winter (down to −22°C). It is classified as wet sub-Mediterranean with an average annual temperature of 9–10°C and annual precipitation of 600–700 mm [17].
The name Rodeno is attributed to the characteristic Permian–Triassic reddish sandstones and conglomerates found in this area (Figure 2a). Other important lithological materials are Ordovician quartzite and Jurassic limestone and dolomites. Ordovician materials were affected by Hercynian orogeny, and the totality of the study area was shaped by the Alpine orogeny [18,19]. Lithology is an important factor affecting plant species distribution. Extensive Pinus pinater Arr. forest (Figure 2b) dominates in sandstones and quartzite areas (i.e., more acid soils), while Quercus ilex L. and Juniperus thurifera L. are more prone to calcaric soils (i.e., soils developed from limestone and dolomite parent materials). In addition to the environmental value of PPPR, the protected area contains an interesting cultural heritage. The presence of numerous rock paintings is particularly noteworthy, which are fortunately protected by law and are recognized as a world heritage by UNESCO (Rock Art of the Mediterranean Basin on the Iberian Peninsula URL: http://whc.unesco.org/en/list/874)
Figure 2.
Pictures of the study area: (a) north view from the Peña de la Cruz peak (see location in Figure 1) and (b) characteristic view of a Pinus pinaster forest.
Geographical information systems (GIS) and remote-sensing techniques are excellent tools for the identification and analysis of landscape spatial patterns [20]. Thus, satellite imagery and cartographic data were analyzed with digital image processing methods and spatial analysis techniques to detect spatial–temporal land-cover changes. A multitemporal Landsat satellite dataset formed the basis for the change detection procedure. Images were acquired by the multispectral Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) sensors onboard Landsat satellites. Images were selected at regular intervals of 8 years from 1984 to 2008. The dates of available scenes were 1984-09-21 (Landsat 4-TM), 1992-04-20 (Landsat 5-TM) the only available scene without cloud cover, 2000-08-08 (Landsat 7-ETM+), and 2008-09-15 (Landsat 7-ETM+ SLC-off). An additional Landsat 7-ETM+ SLC-off acquired on 2008-10-01 was employed for filling gaps in the 2008-09-15 Landsat 7 image. Digital image processing procedure included preprocessing of satellite multispectral images to ensure the temporal comparability (spatial and radiometric coherence and filling the gaps of SLC-off images) between scenes, followed by the calculation of a vegetation index and the application of Fourier Transform analysis to assess temporal changes of land covers.
In addition, several cartographic data sources were employed. For instance, the official protected area boundaries were obtained from the Aragon Territorial Information System (SITAR http://sitar.aragon.es) of the Aragon autonomous community government. The data comprised vector perimeters of the protected area in 1995 and after its enlargement in 2007. Furthermore, a land-cover map was derived from the National Forest Map of Spain [21]. Seven land-cover classes were defined according to dominant vegetation classes included in the vector cartography (Table 1). Pinus pinaster L stands is the most characteristic vegetation type within the protected area.
Land cover class
Description
Pinus sp.
Pine forest dominated by Pinus pinaster L at non-calcareous soils. Pine forest at calcareous soils is denoted by the presence of Pinus nigra J.F.Arnold stands.
Juniperus thurifera
Pure stands of small trees of Juniperus thurifera L developed at calcareous soils.
Quercus ilex
Quercus ilex L. stands preferentially developed at calcareous soils. Usually in the presence of Juniperus thurifera.
Juniperus mixture
Formations of medium-sized shrubs dominated by Juniperus thurifera and Juniperus communis L. They have been intensely grazed and are under natural vegetation processes.
Shrubs
Xerophyte shrubs and pastureland at recursively grazed areas.
Rocky
Non-calcareous rocky areas with the presence of small herbs and shrubs such as Cistus laurifolius L. or Erica scoparia L.
Agricultural
Non-irrigated cereal crops.
Table 1.
Descriptions of land cover classes.
2.1. Preprocessing
Satellite image preprocessing included geometric and atmospheric corrections with the aim to ensure the spatial comparability with other data sources and to obtain at-ground reflectance pixel spectra. Various geo-referenced data types were used for the geometric correction: aerial ortho-photos (1 m pixel resolution) and digital cartography (scale = 1:25000). The 2000 Landsat 7 ETM+ scene was selected because of its very good visual quality. It was geometrically corrected using Ground Control Points (GCP) identified on the ortho-photos and cartographic maps. A quadratic mapping function of polynomial fit and the nearest-neighbor resampling method were used for the correction. The nearest-neighbor resampling method was selected because it ensures that the original (raw) pixel values are retained in the resulting output image, which is an important requirement in any change detection analysis [8]. The maximum allowable root mean square error (RMSE) of the geometric correction was less than half a pixel, a reference value frequently cited as acceptable [5,22,23].
Since May 2003, Landsat 7 images suffer a radiometric problem caused by a malfunctioning scan line corrector (SLC), which compensates for the forward motion of the satellite, and subsequent efforts to recover the SLC were unsuccessful. Fortunately, Landsat 7 is still capable of acquiring useful image data with the SLC turned off, particularly within the central portion of any given scene, and various interpolation and compositing schemes have been developed to expand the coverage of useful data [24]. A local linear histogram-matching method using two SLC-off images (i.e., 15th September 2008 and 1st October 2008 Landsat 7-ETM+ scenes) was applied to fill the gaps [25]. A rescaling function is derived after applying a local linear histogram matching in a moving window over each missing pixel. A rescaling function is then used to convert the radiometric values of a single input scene into equivalent radiometric values of the scene being gap filled, and the transformed data are then used to fill the gaps of that scene [26]. The adaptive local linear histogram adjustment algorithm can yield good results if the scenes have comparable seasonal conditions (the scenes were taken within the shortest revisit time of the satellite over the study area, i.e., 16 days) and do not exhibit radical differences in target radiance (i.e., presence of clouds, snow, sun glint, or changes in solar illumination geometry) [25].
Atmospheric correction involves the estimation of the atmospheric optical characteristics at the time of image acquisition [27]. Radiometric calibration was applied prior to the atmospheric correction. The conversion of raw digital numbers (DNraw) of Landsat level 1 (L1) image products to at-satellite radiance values (Lsat) required the application of current rescaling values [28] by applying the following expression [24,28,29]:
Lsat=(LMAXλ−LMINλ255)⋅DN+LMINλE1
Where Lsat is the at-satellite radiance [W/(m2 sr μm)]; Lminλ is the spectral radiance that is scaled to Qcalmin [W/(m2 sr μm)] (Qcalmin is the minimum quantized calibrated pixel value, i.e., DN = 0, corresponding to Lminλ); Lmaxλ is the spectral radiance that is scaled to Qcalmax [W/(m2 sr μm)] (Qcalmax is the maximum quantized calibrated pixel value, i.e., DN=255, corresponding to Lmaxλ); and DN are digital numbers of the L1 image product. Surface reflectance values (ρ) were computed using the image-based COST method [30] and applied according to Melendez-Pastor et al. [31]. They computed the path radiance (Lp) values using the equation reported in Song et al. [32] and assuming a 1% surface reflectance for dark objects [30,33,34]. The optical thickness for Rayleigh scattering (τr) was estimated according to the equation given in Kaufman [27].
2.2. Vegetation index
The spectral reflectance of the vegetated surface covers recorded by remote sensors can be compressed into vegetation indexes [9], which are directly related to plant vigor, density, and growth conditions and may also be used to detect unfavorable environmental conditions [4]. Vegetation indexes provide an effective way to perform valuable large spatial-and temporal-scale analyses with large amounts of satellite data [10]. The normalized difference vegetation index (NDVI) [12] has been frequently employed for vegetation status assessment and is formulated as follows:
NDVI=ρNIR−ρredρNIR+ρredE2
Where ρNIR is the reflectance of the near infrared spectral band and ρRED is the reflectance of the red spectral band. NDVI has been related with several vegetation parameters such as changes in the amount of green biomass and chlorophyll content, leaf water content, CO2 net flux, absorbed photosynthetically active radiation (APAR), leaf area index (LAI), and many others [35]. It is also employed as vegetation status data source in many environmental modeling approaches [36,37].
2.3. Fourier transform
Satellite NDVI time series can be analyzed by means of the Fourier Transform in order to obtain its frequency domain components [11]. Harmonic or Fourier analysis allows the decomposition of temporal data to the frequency domain, by which frequency information is represented in terms of the sum of a set of sine and cosine functions [38]. Harmonic analysis allows a complex curve to be expressed as the sum of a cosine waves (terms) and an additive term [39]. The Fourier Transform converts a function f(x) to the frequency domain by a linear combination of trigonometric functions [40].
Input variables were checked for normal distribution with the Kolmogorov–Smirnov test. Required variable transformations were applied to ensure the normality of input variables. A one-way analysis of variance (ANOVA) test was employed to detect the NDVI or harmonic component differences using the protected area declaration (PAD) as factor. The same procedure was applied to detect NDVI or harmonic component differences using land-cover classes (LC) as a factor. In addition, a two-way ANOVA test was employed to assess the harmonic component differences using protected area declaration (PAD) and land-cover classes (LC) as factors.
The PPPR study area is a mountainous area with a landscape dominated by the presence of large coniferous forest (Pinus sp.). Figure 3a shows a false color composite RGB = 742 of the 1984 Landsat TM image. Green areas are associated with coniferous forest, whitish areas correspond to almost bare soils of agricultural fields, and reddish areas correspond to sandstone areas with different vegetation classes. A digital elevation model (DEM) shown in Figure 1 allowed the topographic characterization of the area. The mean altitude and average slope of the 1995 protected area and the 2007 enlargement area are very similar. The mean altitude and average slope of the PPPR are 1315 m and 13.7°, respectively.
Figure 3.
Satellite imagery and cartographic data: (a) false color composite RGB = 742 of the 1984 Landsat satellite image and (b) land cover map from the National Forest Map.
The major land cover class is Pinus sp. forest that accounts for 80% (5439.69 ha) of the PPPR (Figure 3b and Table 2). The proportion was slightly higher for the 2007 enlargement area (88% or 3147.39 ha) compared to that for the 1995 perimeter (70% or 2292.30). Shrub (432.9 ha), rock (405.99 ha), and agriculture (371.52 ha) classes represent about 6% each. Juniperus thurifera land cover class is not present in the 1995 perimeter. Meanwhile, Quercus ilex, Juniperus mixture, and shrub land-cover classes are no present in the 2007 enlargement.
Land cover clases
Current PPPR (ha)
1995 perimeter (ha)
2007 enlargement (ha)
Pinus sp.
5439.69
2292.30
3147.39
Juniperus thurifera
18.72
No presence
18.72
Quercus ilex
97.02
97.02
No presence
Juniperus mixture
65.43
65.43
No presence
Shrubs
432.90
432.90
No presence
Rocky
405.99
108.72
297.27
Agriculture
371.52
268.65
102.87
TOTAL
6831.27
3265.02
3566.25
Table 2.
Land cover surface (ha) distribution in the current PPPR, the original protected area (1995) and 2007 park enlargement.
The average NDVI values for each land-cover class (Table 3) were computed for each satellite image (from 1984 to 2008). Pinus sp. and rocky areas had the largest average NDVI values with 0.53 and 0.51, respectively. On the contrary, shrubs and agricultural fields had average NDVI values less than 0.36. Juniperus thurifera, Quercus ilex, and Juniperus mixture stands had sequentially decreasing NDVI values from 0.43 to 0.37.
Variables
x¯±σ
PAD
Land cover
NDVI 1984
0.50 ± 0.31
0.000 ***
0.000 ***
NDVI 1992
0.59 ± 0.35
0.000 ***
0.000 ***
NDVI 2000
0.46 ± 0.28
0.000 ***
0.000 ***
NDVI 2008
0.51 ± 0.32
0.000 ****
0.000 ***
A0
0.51 ± 0.30
0.054 N.S.
0.008 **
A1
0.04 ± 0.59
0.000 ***
0.000 ***
A2
0.05 ± 0.62
0.000 ***
0.000 ***
P1
1.95 ± 1.92
0.060 N.S.
0.220 N.S.
P2
5.30 ± 0.92
0.000 ***
0.000 ***
Table 3.
One-way analysis of variance results (p-values) for NDVI and harmonic components images.
** = p ≤ 0.010 ; *** = p ≤ 0.001 ; N.S. = not significant.
Independent ANOVA test for protected area declaration (PAD) and land cover class factors. Descriptive statistics (mean ± standard deviation) of variables are included. Phase values are in radians.
One-way ANOVA of the NDVI values of the sampled points allowed the independent comparison of the effect of protected area declaration and land-cover classes on vegetation index values (Table 3). Significant differences of NDVI values were observed for both factors. Different land cover types tend to report different values of NDVI, thus allowing the identification of temporal patterns of the vegetation index for each land cover type. In addition, the average NDVI values for both protected area declaration (PAD) subzones could be different by their dissimilar landscape structure. These changes in the vegetation index values are related to the phenological changes in the land covers being characteristic of a particular type of vegetation and may be affected by human actions [11]. Several authors have noted the usefulness of vegetation time series for classifying land cover classes in mountainous areas such as the study area [42] and for detecting changes in forest stands [43,44].
The first two harmonic terms and the additive term were derived from the Fourier Transform of the NDVI time series (Figure 4). Harmonic analysis allowed the conversion of the NDVI time series from the time domain to the frequency domain, which is very useful to mitigate noise effects and gain information about the time at which NDVI changes occurred. Such information could be obtained after a combined interpretation of the phase and amplitude components. The first three or four harmonic terms concentrate vegetation phenology patterns, while the high-frequency noise is partitioned into the higher harmonic terms [41]. The additive term (A0) corresponds with the average NDVI of the time series.
Figure 4.
Cartography depiction of NDVI-based Fourier Transforms: (a) additive or amplitude 0 image, (b) amplitude of the first harmonic component, and (c) phase of the first harmonic component.
The average values of the amplitude of the first harmonic component (A1) and the second harmonic component (A2) are shown for each land cover class (Figure 5a). Highest amplitude values were observed for Pinus sp., rocky and agricultural land cover classes suggesting larger temporal variations of NDVI values compared to those for the other land cover classes. Conversely, Quercus ilex stands seem to be highly stable throughout the time that is particularly important for a low-growth rate species. Juniperus thurifera and Juniperus mixture land cover classes had moderate A1 values suggesting a progressive but slow regrowth of vegetation due to cessation of rural activities.
Figure 5.
NDVI-based Fourier Transform results: (a) mean value of the amplitude (with standard deviation error bars) for the first two harmonic components for each land cover class; and (b) mean phase value of the first harmonic component phase (in years) for each land cover class.
A two-way ANOVA test with the sampled values of the first two harmonic components and the additive term was performed. Protected area declaration and land cover class were employed as factors. Only land cover classes present in both PPPR protected areas (1995 original area and 2007 enlargement) were considered. A total of 2502 pixel of Pinus sp., rocky and agricultural land cover classes were analyzed (Table 4). In relation with the protected area declaration factor, all amplitude and phase variables exhibited significant variations except for the second harmonic component. In relation with the land cover factor, significant differences were observed for zero term and the phase values but not for the other amplitude terms (A1 and A2). Finally, the effects of the interaction of both factors were assessed. Only significant differences were observed for the phase of the first component. A two-way ANOVA test with NDVI values was discarded because such variables did not show any relevant information. Significant differences of NDVI values for PAD and LC factors were always detected.
Factor
Variable
df
F
p-value
PAD
A0
1
14.374
0.000 ***
A1
1
8.221
0.004 *
A2
1
2.145
0.143 NS
P1
1
4.746
0.029 *
P2
1
1.985
0.159 NS
Land cover
A0
2
260.635
0.000 ***
A1
2
2.867
0.057 NS
A2
2
1.661
0.190 NS
P1
2
14.583
0.000 ***
P2
2
4.833
0.008 **
PAS | Land cover
A0
2
2.958
0.052 NS
A1
2
1.297
0.274 NS
A2
2
1.735
0.177 NS
P1
2
24.030
0.000 ***
P2
2
3.697
0.025*
Table 4.
Two-way analysis of variance results for the harmonic component images.
* = p ≤ 0.050; ** = ≤ 0.010; *** = p ≤ 0.001 ; N.S. = not significant.
Protected area declaration (PAD) and land cover class are the factors. Degrees of freedom, F statistic, and p-values are included.
A final interpretation of the Fourier Transform results was obtained after considering the results of the visual inspection of the images complemented with ancillary geographical information and the two-way ANOVA results. The amplitude of the first harmonic component (A1) was primarily associated with changes in the forest stands by logging and the significant cessation of grazing. We analyzed the A1 image in a GIS and located the highest numeric values areas (dark red areas in Figure 4b). Such areas were identified on aerial photography from 1997 and 2008 available at the aforementioned SITAR. We observed a progressive vegetation re-growth of such areas after the cessation of past clear-cutting (currently logging activities are less aggressive, i.e., horses are even used to transport logs to minimize environmental impact) and the progressive decrease of grazing pressure by the rural abandonment. The phase of the first harmonic component (P1; Figure 4c) was less than π (in radians or 180 degrees) for such high A1 values suggesting that the peak of the first harmonic component for such areas was before 1996 (Figure 5). The average values of the phase 1 component were computed as a function of the land cover classes within the entire PPPR. The peak of the first harmonic component was for the late 1980s and early 1990s with later values for Pinus sp., rocky and agricultural land cover classes. The evolution of the population is a factor commonly used to explain the dynamics of the landscape [45]. The current landscape of the Spanish Mediterranean mountain areas is the result of multiple biophysical agents, such as demographic trends, and the political and economic dynamics [46]. The population of the Sierra de Albarracín is continuously decreasing since the early XX century due to rural abandonment. This process of decline in rural population results in a reduction of an extensive livestock activity that is largely responsible for the maintenance of shrublands and pastures [47] and promoting forest regeneration [7]. In the last 25 years, Spain has undergone deep social and administrative changes (rural abandonment, new forestry policies, and administrative changes) promoting vegetation recovery by the reduction of grazing pressure as a result of depopulation of rural areas [48]. Nevertheless, agricultural and logging activities still persist within the park but are well regulated.
Finally, the application of the Fourier Transform was particularly interesting to extract subtle information. The solely NDVI analysis suggested large temporal and land cover-dependent differences on the vegetation cover, but subtle changes seemed to be undetectable. The harmonic analysis of the NDVI time series revealed highly valuable and easily interpretable information. Higher amplitude values (A0, and A1) were observed for the three land cover classes (Pinus sp., rocky and agricultural) within the 2007 enlargement area. The two-way ANOVA revealed significant differences of the additive term and the first harmonic component (A1 and P1) as a function of the PAD factor (Table 4). The additive term also have significant differences as a function of the land cover factor (higher values for the forest cover). The application of the two-way ANOVA was particularly useful to understand the combined effect of land protection and land-cover classes on vegetation dynamics. We suggest that our first harmonic component was a good indicator of vegetation recovery processes. The fact that no significant differences of the A1 component were obtained for the land cover factor but were significant for the PAD factor suggests that the expansion of the park had a positive impact on vegetation coverage in the reclaimed areas. Such differences could not be directly attributable to a poor management of the protected area because the northern sector of the park is closest to the major population centers of the shire (Albarracín and Gea de Albarracín), while the southern sector has greater inaccessibility (worst road network). The presence of a lumber mill in Albarracín (northern study area) suggests larger logging in the northern sector (better accessibility) compared to that in the southern sector of the protected area. The distance to urban areas or roads is a critical factor affecting land-cover changes induced by processes such as deforestation [49,50] land abandonment [51] or urban growth [52].
Future research should be focused on updating the Fourier Transform time series analysis with the inclusion of additional images for September-October 2016 (in order to maintain the 8-year interval). In order to continue the Landsat time-series analysis, Landsat-8 imagery (available from April 2013 on) would be the most appropriate choice as this sensor has similar characteristics to Landsat 7. However, their comparability is somewhat complex by the reflectance differences expected by their notably dissimilar specifications (i.e., narrower bands, 16 bits instead of 8 bits for the OLI (operational land imager) sensor) [53]. It is difficult to comprehensively model their spectral reflectance differences, which also depend on the surface reflectance and atmospheric state [54]. Recent studies have modeled the relationship between the vegetation indices images obtained from ETM+ and OLI [54–56]. This information is very important for future applications of remote sensing to assess landscape temporal dynamics. Vegetated areas had better NDVI agreement than non-vegetated surfaces (especially water areas) [55], and the seasonal agreement of both sensors is better for forest and shrub areas during growth periods than for others land covers [56]. In addition, OLI images are very infrequently saturated [54], and hence more precise quantitative applications of remote sensing are expected for complex landscapes such as wetlands, mountainous areas and forest.
This study applied the Fourier Transform technique to analyze land-cover changes based on a NDVI time series of satellite images in a mountain forestland, considering the expansion stages of a protected reserve. Although the complexity of the study area (geomorphology) could affect the expected results, the harmonic analysis provided a highly efficient method to minimize the time-series noise effect and to extract valuable information for interpreting spatial and temporal changes of land covers. The first two harmonic components plus the additive term were selected for the evaluation of the landscape dynamics. Higher amplitude values were obtained for very dynamic land cover classes or drastic landscape changes (e.g., intensive logging areas), while lower amplitude values were associated with slow regeneration vegetation classes. Further statistical analysis revealed significant differences of NDVI harmonic components for the expansion stages of the protected reserve and land-cover factors. Those results suggested a complex landscape dynamics greatly influenced by the land management and vegetation classes. The interpretation of such landscape changes and underlying driving factors was supported by previous research of key socioeconomic and environmental factors. The study area is suffering serious depopulation problems manifested by the progressive expansion of forest areas at the expense of former grazing land. Awareness of environmental and cultural values of the study area has led to the progressive protection of the area to improve conservation and management. In this sense, the use of remote sensing images to explore and understand the temporal evolution of ecosystems was an excellent source of knowledge to improve land management. The Landsat scenes selected to analyze vegetation cover changes between the four time periods in this Spanish Mediterranean mountainous protected area provided valuable information usable to understand landscape dynamics. The combined use of land cover and protected area development maps with the first harmonic components enabled the analysis of landscape dynamics as influenced by rural abandonment. The multitemporal NDVI analysis with the Fourier transform allows assessing the places and the time that logging activities have taken place. The use of these techniques can help improve the study and ecological conservation in remote rural areas throughout the world.
The authors would like to thank the U.S. Geological Survey for providing the Landsat satellite imagery used in the study.
References
1.IGBP. Land-Use and Land-Cover Change (LUCC) Implementation Strategy. IGBP Report No. 48/IHDP Report No 10. Stockholm, Sweden: International Geosphere-Biosphere Programme; 1999. 126 p.
2.Sala OE, Chapin III FS, Armesto JJ, Berlow E, Bloomfield J, Dirzo R. Global biodiversity scenarios for the year 2100. Science. 2000;287:1770–1774.
3.Lupo F, Linderman M, Vanacker V, Bartholomé E, Lambin EF. Categorization of land-cover change processes based on phenological indicators extracted from time series of vegetation index data. International Journal of Remote Sensing. 2007;28:2469–2483.
4.Li J, Lewis J, Rowland J, Tappan G, Tieszen LL. Evaluation of land performance in Senegal using multi-temporal NDVI and rainfall series. Journal of Arid Environments. 2004;59:463–480.
5.Jensen JR. Introductory Digital Image Processing (3rd Edition). Upper Saddle River, NJ, USA: Prentice Hall; 2005. 544 p.
6.Zhang X, Friedl MA, Schaaf CB, Strahler AH, Hodges JCF, Gao F. Monitoring vegetation phenology using MODIS. Remote Sensing of Environment. 2003;84:471–475.
7.Melendez-Pastor I, Hernández EI, Navarro-Pedreño J, Gómez I. Socioeconomic factors influencing land cover changes in rural areas: the case of the Sierra de Albarracín (Spain). Applied Geography. 2014;52:34–45.
8.Mather PM, Koch M. Computer Processing of Remotely-Sensed Images: An Introduction (4th Edition). Chichester, UK: Wiley-Blackwell; 2011. 460 p.
9.Myneni RB, Hall FG, Sellers PJ, Marshak AL. The interpretation of spectral vegetation indexes. IEEE Transactions on Geoscience and Remote Sensing. 1995;33:481–486.
11.Melendez-Pastor I, Navarro-Pedreño J, Koch M, Gómez I, Hernández EI. Land-cover phenologies and their relation to climatic variables in an anthropogenically impacted Mediterranean coastal area. Remote Sensing. 2010;2:697–716.
12.Rouse JW, Hass RH, Schell JA, Deering DW. Monitoring vegetation systems in the great plains with ERTS. In: Proceedings, Third Earth Resources Technology Satellite-1 Symposium. Washington, DC, USA: NASA SP-351; 1974. p. 3010–3017.
13.Justice CO, Holben BN, Gwynne MD. Monitoring East African vegetation using AVHRR data. International Journal of Remote Sensing. 1986;7:1453–1474.
14.Moody A, Johnson DM. Land-surface phenologies from AVHRR using the discrete fourier transform. Remote Sensing of Environment. 2001;75:305–323.
15.Immerzeel WW, Quiroz RA, De Jong SM. Understanding precipitation patterns and land use interaction in Tibet using harmonic analysis of SPOT VGT-S10 NDVI time series. International Journal of Remote Sensing. 2005;26:2281–2296.
16.Jakubauskas ME, Legates DR, Kastens JH. Crop identification using harmonic analysis of time-series AVHRR NDVI data. Computers and Electronics in Agriculture. 2002;37:127–139.
17.Government of Aragón. Atlas Climático de Aragón. Zaragoza, Spain: Department of Enviornment Government of Aragón; 2007. 291 p.
18.IGME. Mapa Geológico de España. E 1:50000. Hoja 566 Cella. Madrid, Spain: Geological and Mining Institute of Spain (IGME); 1983. p. 70.
19.IGME. Mapa Geológico de España. E 1:50000. Hoja 589 Terriente. Madrid, Spain: Geological and Mining Institute of Spain (IGME); 1983. p. 82.
20.Gumbricht T, Mccarthy J, Mahlander C. Digital interpretation and management of land cover - a case study of Cyprus. Ecological Engineering. 1996;6:273–279.
21.MMA. Mapa Forestal de España 1:200.000. Teruel 7-6. Madrid, Spain: Ministry of Environment (MMA); 1991.
22.Townsend PA, Walsh SJ. Remote sensing of forested wetlands: application of multitemporal and multispectral satellite imagery to determine plant community composition and structure in southeastern USA. Plant Ecology. 2001;157:129–151.
23.Melendez-Pastor I, Navarro-Pedreño J, Koch M, Gómez I. Applying imaging spectroscopy techniques to map saline soils with ASTER images. Geoderma. 2010;158:55–65
24.Irish RR. Landsat 7 Science Data Users Handbook. Landsat Project Science Office, Goddard Space Flight Center-NASA [Internet]. 2011. Available from: http://landsathandbook.gsfc.nasa.gov [Accessed: 2016-03-21]
26.Chen J, Zhu X, Vogelmann JE, Gao F, Jin S. A simple and effective method for filling gaps in Landsat ETM+ SLC-off images. Remote Sensing of Environment. 2011;115:1053–1064.
27.Kaufman YJ. The atmospheric effect on remote sensing and its corrections. In: Asrar G, editor. Theory and Applications of Optical Remote Sensing. New York, USA: Wiley-Interscience; 1989. p. 336–428.
28.Chander G, Haque MO, Micijevic E, Barsi JA. A Procedure for radiometric recalibration of Landsat 5 TM reflective-band data. IEEE Transactions on Geoscience and Remote Sensing. 2010;48:556–574.
29.Chander G, Markham B. Revised Landsat-5 TM radiometric calibration procedures and postcalibration dynamic ranges. IEEE Transactions on Geoscience and Remote Sensing. 2003;41:2674–2677.
30.Chavez Jr PS. Image-based atmospheric corrections - Revisited and improved. Photogrammetric Engineering & Remote Sensing. 1996;62:1025–1036.
31.Melendez-Pastor I, Hernández EI, Navarro-Pedreño J, Gómez I. Mapping soil salinization of agricultural coastal areas in southeast spain. In: Escalante B, editor. Remote Sensing—Applications. Rijeka, Croatia: InTech; 2012. p. 117–140.
32.Song C, Woodcock CE, Seto KC, Lenney MP, Macomber SA. Classification and change detection using Landsat TM data: when and how to correct atmospheric effects? Remote Sensing of Environment. 2001;75:230–244.
33.Moran MS, Jackson RD, Slater PN, Teillet PM. Evaluation of simplified procedures for retrieval of land surface reflectance factors from satellite sensor output. Remote Sensing of Environment. 1992;41:169–184.
34.Chavez Jr PS. Radiometric calibration of Landsat Thematic Mapper multispectral images. Photogrammetric Engineering & Remote Sensing. 1989;55:1285–1294.
35.Jensen JR. Remote Sensing of the Environment: An Earth Resource Perspective. Upper Saddle River, NJ, USA: Prentice Hall; 2007. p. 608
36.Chen SH, Su HB, Tian J, Zhang RH, Xia J. Estimating soil erosion using MODIS and TM images based on support vector machine and à trous wavelet. International Journal of Applied Earth Observation and Geoinformation. 2011;13:626–635.
37.Melendez-Pastor I, Navarro-Pedreño J, Koch M, Gómez I, Hernández EI. Evaluation of land degradation after forest fire with a fuzzy logic model. Environmental Engineering and Management Journal. 2013;12:2087–2096.
38.Briggs WL, Henson VE. The DFT: An Owners’ Manual for the Discrete Fourier Transform. Philadelphia, PA, USA: Society for Industrial Mathematics; 1995. 450 p.
39.Rayner JN. An Introduction to Spectral Analysis. London, UK: Pion Ltd.; 1971. p. 174
40.Sakamoto T, Yokozawa M, Toritani H, Shibayama M, Ishitsuka N, Ohno H. A crop phenology detection method using time-series MODIS data. Remote Sensing of Environment. 2005;96:366–374.
42.Evans JP, Geerken R. Classifying rangeland vegetation type and coverage using a Fourier component based similarity measure. Remote Sensing of Environment. 2006;105:1–8.
43.Parent MB, Verbyla D. The browning of Alaska’s Boreal Forest. Remote Sensing. 2010;2:2729–2747.
44.Yen P, Ziegler S, Huettmann F, Onyeahialam AI. Change detection of forest and habitat resources from 1973 to 2001 in Bach Ma National Park, Vietnam, Using Remote Sensing Imagery. International Forestry Review. 2005;7:1–8.
45.Izquierdo AE, De Angelo CD, Aide TM. Thirty years of human demography and land-use change in the Atlantic Forest of Misiones, Argentina: an evaluation of the forest transition model. Ecology and Society. 2008;13:3.
46.Gómez Vargas FJ, Boada Juncà M, Sànchez Mateo S. Local models of land use & land cover change analysis. The case of Ridaura sessile oak forestland (Natural Park of Montseny. Barcelona). Boletín de la Asociación de Geógrafos Españoles. 2008;47:379–386.
47.García-Ruiz JM, Lasanta T, Ruiz-Flano P, Ortigosa L, White S, González C, et al. Land-use changes and sustainable development in mountain areas: a case study in the Spanish Pyrenees. Landscape Ecology. 1996;11:267–277.
48.Valbuena-Carabaña M, López de Heredia U, Fuentes-Utrilla P, González-Doncel I, Gil L. Historical and recent changes in the Spanish forests: a socio-economic process. Review of Palaeobotany and Palynology. 2010;162:492–506.
49.Chomitz KM, Gray DA. Roads, land use, and deforestation: a spatial model applied to Belize. The World Bank Economic Review. 1996;10:487–512.
50.Southworth J, Marsik M, Qiu Y, Perz S, Cumming G, Stevens F, et al. Roads as drivers of change: trajectories across the tri-national frontier in MAP, the Southwestern Amazon. Remote Sensing. 2011;3:1047–1066.
51.Díaz GI, Nahuelhual L, Echeverría C, Marín S. Drivers of land abandonment in Southern Chile and implications for landscape planning. Landscape and Urban Planning. 2010;99:207–217.
52.Su S, Xiao R, Zhang Y. Multi-scale analysis of spatially varying relationships between agricultural landscape patterns and urbanization using geographically weighted regression. Applied Geography. 2012;32:360–375.
53.USGS. Landsat 8 (L8) Data Users Handbook. Version 1.0. Sioux Falls, SD, USA: Earth Resources Observation and Science (EROS) Center. US Geological Survey; 2015. p. 106.
54.Roy DP, Kovalskyy V, Zhang HK, Vermote EF, Yan L, Kumar SS, et al. Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity. Remote Sensing of Environment. 2016;in press. doi:10.1016/j.rse.2015.12.024
55.Ke Y, Im J, Lee J, Gong H, Ryu Y. Characteristics of Landsat 8 OLI-derived NDVI by comparison with multiple satellite sensors and in-situ observations. Remote Sensing of Environment. 2015;164:298–313.
56.She X, Zhang L, Cen Y, Wu T, Huang C, Baig HM. Comparison of the continuity of vegetation indices derived from Landsat 8 OLI and Landsat 7 ETM+ data among different vegetation types. Remote Sensing. 2015:7:13485–13506
Written By
Ignacio Melendez-Pastor, Encarni I. Hernández, Jose Navarro-
Pedreño, Ignacio Gómez and Magaly Koch
Open access peer-reviewed chapter
