Spatiotemporal Analysis of LANDSAT Satellite Imagery for Change Detection in Česma Forest Ecosystem

: Development of remote sensing and increased availability of satellite imagery of different spatial, spectral, temporal and radiometric characteristics makes information obtained from such sources of vital importance for studying and mapping vegetation. Vegetation indices have a significant role in vegetation change detection and tracking, whether in quantity or quality terms. Each index has specific significance and performance characteristics. A multiple regression statistical analysis of average vegetatio n index values (NDVI, NDWI, GNDVI, EVI and SAVI) was performed for 2012, 2013 and 2014 periods of a wider Česma forest area n ear Vrbovec, Croatia. Further, rasters using three-year average values, sums, variances and standard deviations for all five indices were created. Differencing of average NDVI index values for years 2005 and 2014 was also performed. Imagery chosen was from the active vegetation period and used as a basis for cluster analysis de tection of significant change areas. Česma forest area was selected due to previous field monitoring and point analysis conducted (2012, 2013 and 2014) that serve as validation for this research. Finally, a raster analysis of select areas, surrounding accumulation dam and encompassing Česma forest area exclusivel y, was conducted. The intention was determining vegetation index change dynamics. Spaciotemporal analysis around accumulation dams determined vegetation changes in dam areas. The advantage of the applied method is that, by using the Principal Components Analysis - PCA, it allows change detection, tracking and monitoring on wide areas more promptly than with other methods.The analysis itself was made using 92 LANDSAT images acquired over a 10-year period. of satellite imagery in assessing and monitoring natural resources, determining land use, tree types and development stages


INTRODUCTION
Remote sensing has a significant role in numerous diverse research fields. Its use is next to mandatory in environmental studies, hydrology, minerology, ecosystem studies and others [1]. Forestry application research was conducted to prove applicability of satellite imagery in assessing and monitoring natural resources, determining land use, tree types and development stages of forests, forest classifications, forest inventories, biomass assessments, monitoring forest conditions and determining changes, forest damage evaluations, early detection and monitoring forest fires, studies and condition assessments following natural events and disasters, wildlife and hunting applications, studying of protected areas and the like [2]. Data obtained from raster analysis were, along with the results of field research, integrated into a comprehensive spatiotemporal database as a basis for geostatistical analysis. This provides an increase in accuracy for change detection tracking in a forest ecosystem when combined with analysis from other relevant scientific areas, and finally results in determining impact on an entire forest ecosystem. Previous field research of the area includes ground water levels survey for 2012-2014 period, and a point raster analysis on 22 piezometer locations in a tenyear time span (more in [3]). A semiautomatic selection and download method of imagery for the 2005-2016 period was implemented using QGIS followed by manual removal of images containing cloud cover or those whose quality did not meet analysis requirements. A total of 92 images form Landsat 7 and Landsat 8 sensors was collected for a 10 year period and used in analysis.
This paper describes how the trends coincide with vegetation indices not atmospherically corrected, what, unlike previous researches, enables tracking changes on wider areas in real time. Prior investigations [3] were used to verify results obtained in this research.

Research Area
Economic unit "Česma" is located between longitudes 16°36' and 16°49' east and latitudes 45°49' and 45°53' north, covering an area of 3423.03 ha. It is positioned in a valley surrounded by mountains Papuk and Psunj east, Kalnik and Medvednica west, Moslavačka gora south and Bilogora north of the valley (Fig. 1). The area was protected in 1982 as a special reserve of forest vegetation. Česma forest encompasses two forest sections. The natural differentiators of those sections are the forest communities of oak trees and hornbeam trees estimated to be 130 years old [4]. Such forest communities are extremely susceptible to climate fluctuations and, thus, to change.
Lowland Česma oak forest was chosen as a test filed due to its specific location of test accumulations dams (barriers placed in 2007 and 2009). The drainage canal network built alongside forest transit lines was connected to old waterways as the main recipients, which caused increased draining of surface waters from the forest, resulting in change of ecological conditions. This resulted in drying out the area, change in vegetation structure and decreased forest regeneration.
As a possible relief measure, foremost for the draught, building and placing dams in dried out watercourses was implemented to contain the water within the forest. To eliminate negative effects of the implemented measure, monitoring of seasonal tree growth, crop yield and leaf mass production was established in the accumulation areas [5]. Thus, a network of piezometric stations was established and ground water levels and regimes were monitored. Current analysis of piezometric measurements shows changes in ground water levels, which led to change manifested in the forest ecosystem [6].

Data Collection and Methodology 1.2.1 LANDSAT
Landsat 7 and Landsat 8 images were used in this research. Landsat ETM+ offers images in eight spectral ranges, wavelengths from 0.45 µm up to 0.9 µm. Surface area covered by one image is around 179 km north-south and 183 east-west direction [7]. Resolution of all but the panchromatic channels is 30 m. Panchromatic channel has a 15 m resolution and serves in spatial resolution improvement and detections.
Landsat 8 mission has scientific characteristics and, unlike its predecessor, has no operational mandate but an open and free access to data. The main goal of the mission is to maintain continuity of observation so that the data is consistent and comparable to previous Landsat implementations. Landsat OLI and TIRS images contain nine spectral ranges of 30 m resolution, with one channel (the eighth) having a 15 m resolution. The new channel 1, of expressly blue spectrum, is used in coastal research, while the new channel 9 is used for cirrus cloud detection. Channels 10 and 11 provide higher accuracy of surface temperatures in 100 m resolution. The approximate scene size is 170 km north-south and 183 east-west direction acquired daily (550 per day on average) [8].

Empirical Vegetation and Water Indices
Both the empirical vegetation and water indices were used in raster analysis of Česma forest. Empirical vegetation indices encompass the SR -simple ratio, NDVI -Normalized Difference Vegetation Index, GNDVI -Green Normalized Difference Vegetation Index, SAVI -Soil-Adjusted Vegetation Index and EVI -Enhanced Vegetation Index. Empirical water indices encompass NDWI -Normalized Difference Water Index and MNDWI -Modified Normalized Difference Water Index.
Precision Leaf Area Index (LAI) surveys are vital to understanding biological and physical processes involving vegetation and are data used in ecosystem modelling [9]. LAI is a nondimensional variable quantifying the leaf to ground surface ratio.
According to research, SR and LAI show the greatest correlation ratio. Its significance was proven in vegetation detection, where it reduces the atmospheric and topographic effects, but has a low ground, ice and water detection response. SR is calculated using [10]: NDVI is the most commonly used indicator in vegetation cover change detection. It is often used as an indicator in tracking vegetation health and dynamics thus allowing spaciotemporal comparisons in vegetation [11]: GNDVI uses the green channel instead of the red one. It shows greater correlation to LAI than NDVI [12]: SAVI is a modified NDVI. It is used in cases when it is necessary to remove the atmospheric and topographic effects in scarce vegetation cover areas [13]: EVI is also a modified NDVI that has increased biomass sensitivity with minimal atmospheric reflection and ground surface effect. Thus, it is a valuable index for monitoring forest vegetation cover [14]: NDWI index was designed to accomplish the following effects: maximum water refraction using the green channel, minimization of low reflection in near-infrared (NIR) spectral range off water bodies along with utilization of high vegetation and earthen objects reflection in the NIR spectral range. The formula for this index is [15]: MNDWI reduces noise effects and improves clear water bodies delineation in imagery. MNDWI is obtained using the following formula [13]: where MIR is the middle infrared spectral range.

Raster Analysis of Vegetation Indices
This research used a multiple regression average value statistical analysis of vegetation indices for years 2012, 2013 and 2014. It is a type of average index values linear regression analysis using each 2012 raster as a nonindependent variable and corresponding 2013 and 2014 rasters as independent (predictor or regression) variables. A new raster based on the regression was derived. The initial assumption for using this analysis was linear variable correlation. The analysis was made for the stated period as it corresponds to the period of piezometer observations. Thus, analysis trends could be compared to in-field conditions and interaction between groundwater and treetop appearance could be determined.
A statistical analysis based on average vegetation indices was made for the 2012-2014 period.
Three-year average value, sums, variances and standard deviations rasters for all five indices were derived. NDVI average value minimum and maximum values are 0.3-0.61, sums 0.34-1.85, variance 0-0.004 and standard deviations −0.013-0.059 (Fig. 2). Vegetation index standard deviation is a value that statistically quantifies the variance and shows the deviation of data from the average. Unlike variance it uses the same units as the data itself. This analysis considers the wider Česma forest area, which affects the average values of the rasters and the results. Hence, subsequent analysis involved the forest area exclusively. An average NDVI index value difference was calculated between 2005-and 2014-year periods. Images used were from the active vegetation period, and the differences were used for a cluster analysis representing a statistical technique of determining relatively homogenous pixel groups. Ten pixel group categories were defined that clearly showcase the greatest changes in eastern and central forest areas (Fig. 3) [16].

Raster Analysis in Accumulation Dam Areas
Finally, a limited area raster analysis including just the accumulation dams' areas within the Česma forest was made. The intention was to show the vegetation index change dynamics on vegetated surfaces exclusively. The accumulation dam area analysis revealed vegetation changes in the dam areas. A mask of exclusively forest areas was created. All the images were then masked to enable a spatiotemporal analysis.
Masked images were stacked to calculate the five vegetation indices for each pixel. A recalculation of NDVI, NDWI, GNDVI, EVI and SAVI was done for each masked image in the dam areas. First raw images where the pixel value was defined as a digital number (also known as DN i.e. numerical value describing grey tone value), and then DOS1 (a well-known and described method) atmospherically corrected images were used. Afterward, mean index values in all 92 images along with yearly averages were calculated (Tab. 1).    Fig. 4 shows monthly EVI values. The monthly values start from top left corner image in December, followed by EVI results for May, June, July, August and September. Those are also the only Landsat 7 images available for that year, as other have too much cloud cover. The highest index values are in spring for all indices as the highest concentration of chlorophyll in plants is in that period, and the same goes for EVI. Fig. 5 shows the NDWI index which indicates water saturation in plants [16]. The vegetation index correlation matrix shows a high correlation ratio of indices in thick vegetation areas (Fig.  6).

Figure 6 Correlation matrix of average vegetation indices values [16]
A mean vegetation index value was derived from all 92 images, limited to the masked area, and depicted in an annual index trend graph in Fig. 7 (with DOS1 correction) and Fig. 8 (without DOS1 correction) for all indices combined. The graphs show the same index trends as the ones derived from a point analysis (more in [3]).

Figure 7
Annual vegetation index trend including DOS1 correction [16] A Pearson analysis and main component analysis on all indices were made exclusively in the forest area. Pearson linear correlation coefficient shows a positive correlation between indices. Fig. 9 shows correlation of all index pairs in 2015, and diagonally positioned are histograms of each index. A visual inspection confirms that the EVI and SAVI indices are completely correlated when DOS1 correction is employed (r = 1). The greatest correlation in each year was obtained between EVI and SAVI, followed by NDVI and GNDVI indices [16].  Main components analysis represents one of the simplest multivariant techniques. It is used when there are multiple same dimension variables that lack any additional information compared to those derived from other variables. It makes visual inspection of derived vegetation indices possible, as well as simultaneously obtaining information on all indices and discerning differences in PCA analysis application for each year and index. As the correlation matrix in Fig. 9 shows, indices are correlated which is beneficial to the PCA analysis. Results of PCA analysis on NDVI can be seen in Fig. 10 for the 2005-2015 period. The nine images represent the initial nine components of all NDVI images and depict representative changes in forest vegetation cover. The first PCA component represents the main variance trying to describe the original data trend Fig. 10. That variance shows the vegetation amount but also the quality of vegetation changes in the environment [16].
The PCA analysis was conducted for 92 images per index to reduce the number of images to a lesser component number. Fig. 11 shows results of a PCA analysis on one of the indices (NDWI). The characteristics of indices in the entire data set were ultimately downsized to the first nine components [16].  [16] This was done as it is one approach applied in vegetation index time series analysis that proved appropriate in reducing large data sets and detecting variations in cycle vegetation trends is PCA. Its most valuable characteristic is efficiency in ground-vegetation delineation allowing automatic extraction of forest areas in satellite imagery. It is a fast and less demanding method compared to image classification that requires numerous parameters input. Spatiotemporal vegetation index variances need to be emphasized as long-term changes, where large datasets are grouped into intercorrelated variables.
Results confirm the great potential of PCA remote sensing data analysis in forest cover research.

CONCLUSION
Vector analysis, index differentiation and pixel grouping were the methods applied to the Česma forest area for change detection aimed at detection and visualization of changes occurring in a select area.
The analysis showed a high correlation between vegetation indices and confirmed the continuity of vegetation trends in all five indices considered in the research. Correlation differences were determined between the indices calculated for a wider and those of the exclusive forest areas. Vegetation indices of forest vegetation areas had double the vegetation index pixel values than those in agricultural areas. The highest correlations were determined in NDVI to GNDVI and SAVI on a wider area, as well as in EVI and SAVI ratios in exclusively forest areas (0,99 -1,00), followed by GNDVI and SAVI and similarly NDVI and GNDVI where the correlation was above 0,8, while NDWI has the lowest correlation ratio in this scenario compared to other indices. Spatial vegetation index trend is primarily manifested in seasonal vegetation activity. For, when vegetation indices are concerned, the result is biased i.e. dependent on spectral characteristics of the observed object.
Vegetation indices were determined using both the DOS1 atmospherically corrected and uncorrected satellite imagery.
The research confirmed that EVI is comprised of values not applicable when raw images without atmospheric correction are used (Fig. 8). The same is not the case in the other four indices that can use uncorrected imagery, thus enabling easier processing of large data sets on wider areas.
Unlike younger forests, tree tops of older stands cover greater surface areas and absorb more water which is clearly identifiable in annual NDWI mean values. EVI shows variations in vegetation structure and is a good type of forest stand indicator as well as forest topography visualization asset because it allows easy identification of young forest areas. Vector change analysis is clearest in EVI. It clearly identified vegetation changes in accumulation dam areas.
We can conclude that prediction of forest ecosystem development based on detailed spatiotemporal vegetation index analysis is possible.
Data processing was done using open source programs like Saga, Semi-Automatic Classification Plugin for QGIS and statistical programming language R. R proved to be a valuable resource in processing 92 raster images needed for a quality geostatistical analysis.
Compared to the previous point analysis research [3], we can conclude that the point and raster analysis show matching trends in NDVI, NDWI, GNDVI and SAVI. Value dissipation of all indices in both point and raster analysis shows a slight downward trend of vegetation activity. This proves that spatiotemporal analysis on an entire raster can produce same trend results like those of the point analysis.
It allows change detection, tracking and monitoring on wide areas more promptly than with other methods Also, comparative advantage of the applied method to point analysis for wide area is that it allows prompt change detection and tracking of vegetation index trends using atmospherically non-corrected imagery.