Descriptions of land cover classes.
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.
- landscape dynamics
- forest management
- rural abandonment
- remote sensing
- Fourier transform
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 . 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 . Remote sensing allows the study of the role of terrestrial vegetation in large-scale global processes (e.g., the carbon cycle) . 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 . In this context, vegetation indices have been used to derive measures that correlate with surface biophysical properties  and as such facilitate the analysis of large amounts of satellite data for spatial- and temporal-scale analyses . 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. , 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 . 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 . 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.
2. Material and methods
The study area is located in the
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 .
Geographical information systems (GIS) and remote-sensing techniques are excellent tools for the identification and analysis of landscape spatial patterns . 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 . Seven land-cover classes were defined according to dominant vegetation classes included in the vector cartography (Table 1).
|Land cover class||Description|
|Pine forest dominated by |
|Pure stands of small trees of |
|Formations of medium-sized shrubs dominated by |
|Shrubs||Xerophyte shrubs and pastureland at recursively grazed areas.|
|Rocky||Non-calcareous rocky areas with the presence of small herbs and shrubs such as |
|Agricultural||Non-irrigated cereal crops.|
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 . 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 . 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 . 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 . 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) .
Atmospheric correction involves the estimation of the atmospheric optical characteristics at the time of image acquisition . 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  by applying the following expression [24,28,29]:
2.2. Vegetation index
The spectral reflectance of the vegetated surface covers recorded by remote sensors can be compressed into vegetation indexes , which are directly related to plant vigor, density, and growth conditions and may also be used to detect unfavorable environmental conditions . Vegetation indexes provide an effective way to perform valuable large spatial-and temporal-scale analyses with large amounts of satellite data . The normalized difference vegetation index (NDVI)  has been frequently employed for vegetation status assessment and is formulated as follows:
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 . 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 . Harmonic analysis allows a complex curve to be expressed as the sum of a cosine waves (terms) and an additive term . The Fourier Transform converts a function
The correct interpretation of harmonic analysis results requires the consideration of several key concepts [11,41]: (1) each wave is defined by a unique amplitude and phase angle, where amplitude is the difference between the maximum value of a wave and the overall average of the time series (usually 0 because the mean is equal to the harmonic term with
2.4. Statistical analyses
Statistical analyses were used to assess the influence of the time of landscape protection (1995 original perimeter of the protected area vs. 2007 park enlargement) and the land-cover classes on the vegetation cover status as denoted by the NDVI images and derived harmonic components. Data for statistical analysis were obtained by randomly selecting pixels on the NDVI and harmonic components images. A total of 2742 pixels within the current perimeter of the PPPR were obtained to extract the NDVI and harmonic component values, and further statistical analyses were performed with the IBM-SPSS statistics 23 © software (IBM Corporation).
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.
3. Results and discussion
The PPPR study area is a mountainous area with a landscape dominated by the presence of large coniferous forest (
The major land cover class is
|Land cover clases||Current PPPR (ha)||1995 perimeter (ha)||2007 enlargement (ha)|
The average NDVI values for each land-cover class (Table 3) were computed for each satellite image (from 1984 to 2008).
|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 ***|
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 . Several authors have noted the usefulness of vegetation time series for classifying land cover classes in mountainous areas such as the study area  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 . The additive term (A0) corresponds with the average NDVI of the time series.
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
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
|Land cover||A0||2||260.635||0.000 ***|
|PAS | Land cover||A0||2||2.958||0.052 NS|
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
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 (
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) . It is difficult to comprehensively model their spectral reflectance differences, which also depend on the surface reflectance and atmospheric state . 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) , and the seasonal agreement of both sensors is better for forest and shrub areas during growth periods than for others land covers . In addition, OLI images are very infrequently saturated , 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.
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.
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.
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.
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.
Jensen JR. Introductory Digital Image Processing (3rd Edition). Upper Saddle River, NJ, USA: Prentice Hall; 2005. 544 p.
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.
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.
Mather PM, Koch M. Computer Processing of Remotely-Sensed Images: An Introduction (4th Edition). Chichester, UK: Wiley-Blackwell; 2011. 460 p.
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.
Tucker CJ, Townshend JRG, Goff TE. African land-cover classification using satellite data. Science. 1985; 227:369–375.
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.
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.
Justice CO, Holben BN, Gwynne MD. Monitoring East African vegetation using AVHRR data. International Journal of Remote Sensing. 1986; 7:1453–1474.
Moody A, Johnson DM. Land-surface phenologies from AVHRR using the discrete fourier transform. Remote Sensing of Environment. 2001; 75:305–323.
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.
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.
Government of Aragón. Atlas Climático de Aragón. Zaragoza, Spain: Department of Enviornment Government of Aragón; 2007. 291 p.
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.
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.
Gumbricht T, Mccarthy J, Mahlander C. Digital interpretation and management of land cover - a case study of Cyprus. Ecological Engineering. 1996; 6:273–279.
MMA. Mapa Forestal de España 1:200.000. Teruel 7-6. Madrid, Spain: Ministry of Environment (MMA); 1991.
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.
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
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]
USGS. Phase 2 gap-fill algorithm: SLC-off gap-filled products gap-fill algorithm methodology [Internet]. 2004 Available from: http://landsat.usgs.gov/documents/L7SLCGapFilledMethod.pdf [Accessed: 2016-03-21]
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.
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.
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.
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.
Chavez Jr PS. Image-based atmospheric corrections - Revisited and improved. Photogrammetric Engineering & Remote Sensing. 1996; 62:1025–1036.
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.
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.
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.
Chavez Jr PS. Radiometric calibration of Landsat Thematic Mapper multispectral images. Photogrammetric Engineering & Remote Sensing. 1989; 55:1285–1294.
Jensen JR. Remote Sensing of the Environment: An Earth Resource Perspective. Upper Saddle River, NJ, USA: Prentice Hall; 2007. p. 608
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.
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.
Briggs WL, Henson VE. The DFT: An Owners’ Manual for the Discrete Fourier Transform. Philadelphia, PA, USA: Society for Industrial Mathematics; 1995. 450 p.
Rayner JN. An Introduction to Spectral Analysis. London, UK: Pion Ltd.; 1971. p. 174
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.
Jakubauskas ME, Legates DR, Kastens JH. Harmonic analysis of time-series AVHRR NDVI data. Photogrammetric Engineering & Remote Sensing. 2001; 67:461–470.
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.
Parent MB, Verbyla D. The browning of Alaska’s Boreal Forest. Remote Sensing. 2010; 2:2729–2747.
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.
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.
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.
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.
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.
Chomitz KM, Gray DA. Roads, land use, and deforestation: a spatial model applied to Belize. The World Bank Economic Review. 1996; 10:487–512.
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.
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.
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.
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.
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
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.
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