MODELING PINUS ELLIOTTII GROWTH WITH MULTITEMPORAL LANDSAT DATA : A STUDY CASE IN SOUTHERN BRAZIL

Remote sensing data are a key proxy to forest monitoring and management at local, regional and global scales. Considering the hypothesis that NDVI and EVI can be used at least during one decade to monitor Pinus elliottii in Southern Brazil, the objective of this study was to identify saturation time after planting of these vegetation indices in a Pinus elliottii plantation and the most suitable index by adjusting theoretical functions to each one of them. Based on Landsat Surface Reflectance Higher-Level Data Products, 32 scenes were selected between 1984 to 2015. A set of theoretical polynomial, gaussian and logistic mathematical functions were applied to fit the experimental data on vegetation indices. The determination coefficient (R2) and RMSE at 95% probability were also used. Finally, EVI efficiency was tested by changing the L parameter. The logistic model was the one that best explained the data resulting from NDVI and EVI over time. NDVI was more effective than EVI for this forest monitoring, identifying the forest growth pattern until its 18 years of age. EVI may have been saturated after 14 years and the L factor may be set to near to zero to achieve a higher coefficient of determination.


Introduction
Planted forest containing Pinus spp.achieved near 2 million hectares in Brazil.In the South of Brazil (Rio Grande do Sul State), because of soil and climate characteristics similar to those of the southeastern region of the United States, Pinus elliottii found adequate conditions of growth (Gholz and Fisher 1982, Izumi et al. 2008, Yuan et al. 2013).Such plantations are mainly intended to be a raw material to pulp and paper industries, mechanical processing and resin extraction (Shimizu 2008).As forest production is related to dynamic variables such as ecological conditions of the planting site, silvicultural treatments, and environmental/climate condition, monitoring systems based on remote sensing techniques are recommended to cover large areas (Günter et al. 2009, Restrepo andOrrego 2015).
Even on a regional scale, the conditions of trees and stands change over time.Therefore, forest monitoring using traditional in-situ methods is a challenge when time and cost effectiveness are considered.Remote sensing data were successfully used to monitoring changes in forests growth and productivity (Park et al. 2016).Moreover, there are terrestrial ecosystem models which can be performed almost totally with satellite-derived data input (Milne and Cohen 1999, Prince and Goward 1995, Bradley et al. 2007).By monitoring forest growth, the commercial wood and resin production can be estimated, allowing a more accurate management of plantations.
To identify and understand possible future forecasts changes in forests, vegetation indices are widely used, which consist of spectral transformations of two or more bands designed to enhance the response of vegetation properties.These indices allow reliable spatial and temporal intercomparisons of terrestrial photosynthetic activity, and canopy structural variations (Wiegand et al. 1991, Huete 2002, Solano et al. 2010).The normalized difference vegetation index (NDVI) (Rouse et al. 1973) and enhanced vegetation index (EVI) (Huete et al. 2002) have been designed for vegetation studies and improve upon the detection of vegetation changes through remote sensing data and extraction of canopy biophysical parameters (Huete 2002).However, both NDVI and EVI show scaling problems, asymptotic (saturated) signals over high biomass conditions (Huete 1988, Solano 2010).EVI was proposed to enhance the vegetation signal and improve sensitivity in high biomass environments.Although EVI decouples the soil and atmospheric influences from the vegetation signal (Huete 2002), the high dependence on near infrared band is an important constraint to be considered (Galvão et al. 2011).
Despite several the adjustments and advances proposed in the literature, most optical reflectance vegetation indices still show saturation problems when in high biomass conditions.Considering this limitation, our hypothesis is that both NDVI and EVI can be used at least for one decade to monitor Pinus elliottii in Southern Brazil and that a theoretical function can be adjusted to simulate the greening process over the time.Such equation can be locally applied to obtain forest quantities.

Study area
A commercial plantation of Pinus elliottii of six thousand hectares located in the municipality of Santa Vitória do Palmar, in Rio Grande do Sul State, Brazil was selected to the study because of the availability of old plantations with planting date archive at a per-plot basis (Figure 1).This region belongs to the coastal region, characterized by flat landforms and formed by unconsolidated sediments under the influence of marine transgression/regression process.The soil is mainly formed by Quartzarenic Neosol (Ustic Quartzipsamments in the Soil Taxonomy).According to the Köppen-Geiger climate classification, the region has a Cfa climate type (Wrege et al. 2011) with an annual average temperature between 16.42 °C and 18.11 °C, and average annual rainfall of 1258.8 mm well distributed along the year.
The slash pine (Pinus elliottii Engelm) plantation was established in 1988 on initial spacing of 2 x 2 m (2,500 trees ha-1).Thinning and pruning silvicultural treatments were not made due their expensive costs and, on the other hand, because the Pinus resin harvesting became most lucrative (FLOPAL 2009).Currently, according to the company forest inventory for exploration, still remains 1,338 trees ha-1 and the average diameter at breast height of 22.4 cm (±5.79) and the total height of 21.4 m (±3.30) (FLOPAL 2017).

Data acquisition
Data on atmospherically corrected surface reflectance from "Surface Reflectance Higher-Level Data Products" were acquired (for visible, near infrared and short wave infrared bands).The surface reflectance images (Landsat 5 Thematic Mapper -TM and Landsat 7 Enhanced Thematic Mapper Plus -ETM+) were produced using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm, a specialized software originally developed through a National Aeronautics and Space Administration (NASA) Making Earth System Data Records for Use in Research Environments (MEaSUREs) grant by NASA Goddard Space Flight Center (GSFC) and the University of Maryland (Masek et al. 2006).
The software applies Moderate Resolution Imaging Spectroradiometer (MODIS) atmospheric correction routines to Level-1 data products.Water vapor, ozone, geopotential height, aerosol optical thickness, and digital elevation are input with Landsat data to the Second Simulation of a Satellite Signal in the Solar Spectrum (6S) radiative transfer models to generate top of atmosphere (TOA) reflectance, surface reflectance, brightness temperature, and masks for clouds, cloud shadows, adjacent clouds, land, and water.
On the other hand, Landsat 8 Surface Reflectance data were generated from the Landsat Surface Reflectance Code (LaSRC) that makes use of the coastal aerosol band to perform aerosol inversion tests, and uses auxiliary climate data from MODIS and uses a unique radiative transfer model.Additionally, LaSRC hardcodes the view zenith angle to "0" (zero), and the solar zenith and view zenith angles are used for calculations as part of the atmospheric correction.As all products were obtained using the EarthExplorer tool (http://earthexplorer.usgs.gov/)which provides products with all the mentioned applied processes, therefore, it was not necessary to perform other atmospheric corrections in the images.
The available surface reflectance scenes were inspected and 32 cloud-free scenes were selected to cover the period from 1984 to 2015 (USGS 2016).To reduce the seasonal (intra-annual) variations of vegetation and view-illumination geometry, only scenes from May to October of each year were selected, in this case, the solar elevation amplitude between the data was 27º, which may affect the EVI and less the NDVI (Galvão et al. 2011).
Monthly accumulated precipitation and average temperature data were acquired from the field meteorological station of the Brazilian National Institute of Meteorology (INMEThttp://www.inmet.gov.br/portal/)located 89 km away from the study area (station coordinates: 33º31'S; 53º21'W and 6 m a.s.l.).The meteorological data were used to identify abnormal periods of precipitation and temperature as a complement to the vegetation analysis (Morellato et al. 2000).
Complementarily, field data were used in association with data on planting date, forest inventory and ancillary information (photos, interview, and applied forest management techniques).The field data were contrasted to imagery-derived information.Unfortunately, inventory data were only available for the year of 2017 (FLOPAL 2017).

Data analysis
The analysis was conducted using 500 pixels randomly sampled within the forest plantation.For that purpose, first all forest plots were mapped and used as a mask to hide other land cover types.As a reference, a SPOT 2015 high resolution spatial image was used, which was assisted with field work.The image was acquired on May 4, 2015, in the multispectral form (IMG_SPOT6_MS_201505041303399_ORT_1357488101_R1C1, with 6 meters of spatial resolution) and panchromatic (IMG_SPOT6_P_201505041303399_ORT_1357488101_R1C1, with 1.5 meters).In order to improve the spatial quality of multispectral image, the Gram-Schmidt fusion was applied (Laben et al. 1998).Ultimately, the average and standard deviation of spectral derived parameters were calculated (Exelis Vis 2013).
After performing the meteorological and spectral characterization of the study area, NDVI and EVI were calculated according to Equations 1 and 2, respectively (Rouse et al. 1973, Huete et al. 2002). (2) where NIR refers to the near infrared band reflectance, R represents the red band reflectance and B, the blue band reflectance of the Landsat series data.L is the canopy background adjustment for correcting nonlinear, differential NIR and red radiant transfer through a canopy; C1 and C2 are the coefficients of the aerosol resistance term (which uses the blue band to correct aerosol influences in the red band); and G is a gain or scaling factor.The coefficients adopted in the EVI algorithm were, L=1, C1=6, C2=7.5, and G=2.5 (Solano 2010).
Average vegetation indices were analyzed over time by graphing time series  and a spatial per-pixels analysis.To discuss the inter-annual variations, meteorological and field data were used.Such approach allowed understanding the dynamic of tree growth within the study area.
A set of theoretical polynomial, Gaussian and logistic mathematical functions were applied to fit the data on experimental vegetation indices.To evaluate the performance of the functions, the determination coefficient (R²) at 95% probability was used.An evaluation was also made by determination of RMSE (Root Mean Square Error), which accounts for the total error of measurement for a given model, defined by the square root of the sum of the variances.This measure assumes that the biggest mistake in estimating the income has a greater proportional weight that reduces errors (Equation 3) (Rojas 2007).
Where, Y is the average output; Ŷ is the average yield estimated by the mathematical model expression; and N is the sample number used in the setting.
Finally, EVI efficiency was tested by changing the canopy background adjustment factor (L) of its equation (Equation 2).Thus, the EVI was re-evaluated for all images considering different values for the parameter (1.00; 0.75; 0.50; 0.25; 0.10) and the resulting coefficient of determination was analyzed and discussed.

Data analysis
Analysis of monthly-averaged accumulated rainfall over 32 years  (Figure 2) presented a mean of 104.1 mm (± 14.47 mm) and a smooth tendency of higher precipitation in the first semester rather than in the second semester.The temperature followed the annual pattern expected for middle latitude with well-defined seasons.Although the average rainfall is well distributed along the year, with at least 80 mm per year per month, Kuinchtner and Buriol (2001) found that the South Coast is one of the regions with less rain than the other regions of the State of Rio Grande do Sul (Figure 2).
Bulletin of Geodetic Sciences, 24(3): 286-, Jul-Sept, 2018 Forest plantation growth was measured through the spectral analysis (Figure 3).In spite of the multispectral resolution of the sensors, clearly the forest growth could be analyzed.The spectra of 1986 represent typical spectral mixture between vegetation and background (sandy soil) and agrees on field data regarding the starting planting year (1985).The initial spacing of 2 x 2 m (2,500 trees ha-1) presented a gently absorption feature in the red band highlighting the presence of vegetation.The five-year step time sampling process caught the variations of the vegetation and displayed a decrease of reflectance in the visible waveband region and an increase in the near infrared in the first 15 years.
A positive relationship exists between vegetation increase and near-infrared reflectance due to strong scattering and weak absorptance by leaves in the canopy (Spanner et al. 1990, Jensen 2009).In addition, as the forest canopy development occurs, the vegetation absorbs more in the red due to the increase of the photosynthesis (Jensen 2009).Nevertheless, after 2001 the spectral behavior changes, since there is no a single factor dominating the forest reflectance.As the vegetation cover reaches a high density, the shadows presence may cause unexpected fluctuations in the reflectance values, presenting smaller reflectance values in denser forests (Carreiras et al. 2006).Therefore, there is threshold above which there is no longer any change in the reflectance values of the canopy.In the shortwave infrared (1300 to 2500nm) the vegetation has a lower reflectance due to strong water content absorption (Jensen 2009).
Figure 3: Spectral response variability of Pinus elliottii after planting (1985) based on Landsat surface reflectance Higher-Level product.A smooth trace fitting algorithm was applied over the 500 pixels average reflectance.

Time series and theoretical modeling
The inter-annual analysis of the vegetation indices from 1984 to 2015 showed a range of 0.22 to 0.88 and 0.16 to 0.53 for NDVI (Figure 4a) and EVI (Figure 4b), respectively.These values show the transition from bare soil predominance to close green canopy.Clearly, NDVI does saturate after EVI, and indicating that this index can better represent the variation of Pinus elliottii over the time.
Conversely, NDVI presented higher fluctuation between adjacent years than EVI.NDVI tends to saturate 17 years after planting.Thereafter it seems insensitive to variations in increased plant biomass (Sellers 1989).In contrast to its formulation, EVI saturated 13 years after planting.In fact, this index is strongly correlated with near-infrared band (Galvão et al. 2011) reducing its applicability to monitoring of this forest.Another important point refers to the red band influence on the indices.In general, the red band reflectance has higher influence over NDVI than EVI, especially in low reflectance conditions (higher photosynthetically active biomass).In such cases, NDVI tends to saturate (Huete et al. 2002).
Bulletin of Geodetic Sciences, 24(3): 286-, Jul-Sept, 2018 Because of the fluctuations in the experimental models (around the average of 500 pixels per year), the adjustment of a theoretical mathematical function allows to better describe the time dynamics of the forest until the saturation of the vegetation indices.Therefore, the performance of some models was evaluated and the best one was selected.The linear and second order polynomial functions returned acceptable R² values; however, for the linear model, higher RMSE has been achieved.The theoretical mathematical function which best described NDVI over the time was a logistic function, R² of 0.95 and the lowest RMSE (Table 1).
For modeling of EVI over time, the same functions were evaluated, and the logistic model yielded the highest coefficient of determination (Table 2).By contrasting the NDVI and EVI functions (Figure 4a and 4b, respectively), it can be clearly seen that EVI saturated before NDVI.Our results suggest that NDVI can be used for almost 18 years while EVI reaches the saturation after 14 years.By obtaining these functions, now one can estimate the growth of Pinus elliottii even few year after planting (considering same or similar environmental conditions).
Another important result refers to the near linear behavior of NDVI before saturation in contrast to EVI.This pattern reduces the influence of external factors (view/illumination geometry, atmospheric correction accuracy, local topographic parameters) (Breunig et al. 2015).Some factors highlighted by Huete et al. (2002), as main influences in the near-infrared band, such as the type of canopy and plant shape of the leaves, may have not lead to many changes in this forest.
In addition, the blue band is used in order to mitigate atmospheric interference; however, plants of the understory to photosynthesis also use it (Kaufman and Tanré 1992).
The saturation of the indices could also be analyzed spatially (Figure 5).In this case, EVI reached the saturation point long before NDVI.A similar pattern was observed throughout the study area where Pinus elliottii was planted.When the EVI efficiency was tested by changing the L parameter, the results presented and increase of the determination coefficient with the decrease of the L parameter used in the EVI equation.The biggest factor L (1.00) had the lowest coefficient of determination.This coefficient increased as the L factors approached zero (Table 3) (Galvão et al. 2011).The L factor near zero enables R² values similar to those presented by NDVI data, which showed a better fit.Thus, this result may suggest that for forests with commercial purposes, factor L may not be necessary close to one as the default set in many products.

Conclusions
The potential of vegetation indices for forest monitoring must be constantly investigated.In this paper we explored the use of remotely sensed time series to study the growth of a Pinus elliottii plantation in southern Brazil.The curve fit methodology presented here allowed temporal intercomparisons of terrestrial photosynthetic activity, and canopy structural variations.It was verified that the logistic model was the one that best explained the data resulting from NDVI and EVI over time.Such equation may be locally applied to obtain forest quantities when in similar environmental conditions.
For monitoring of this forest, NDVI was more effective than EVI, as it was more sensitive to increased photosynthetic biomass over the years, taking longer to saturate.NDVI has identified the forest growth pattern until its 18 years of age; after that, the index saturated.EVI may have been saturated after 14 years.The L factor used in EVI calculation may be set to near zero to achieve a higher coefficient of determination.

Figure 1 :
Figure 1: Location of study area in Rio Grande do Sul State, south of Brazil.The area belongs to the company Flopal Ltda.A false-color composition (RGB654) acquired by Landsat 8 OLI sensor in 2015 is shown.Patch/row 221-083.

Figure 2 :
Figure 2: Average monthly accumulated precipitation and temperature (1984 and 2015), measured in Santa Vitória do Palmar, RS meteorological station.For rainfall, the error bars were displayed and for temperature, a b-spline interpolation was applied to the standard deviation bars.

Figure 5 :
Figure 5: The spatial assessment of vegetation index for the first 20 years.

Table 1 :
Theoretical mathematical models adjustment to NDVI over the time.NDVI refers to an average of 500 pixels for each year.N refers to the number of years with valid images.

Table 2 :
Theoretical mathematical models adjustment to EVI over the time.EVI refers to an average of 500 pixels for each year.N refers to the number of years with valid images

Table 3 :
Coefficient of determination (R²) from different canopy background adjustment factors (L) of the EVI equation.