Modeling leaf area index using time-series remote sensing and topographic data in pure Anatolian black pine stands.

We aimed to map and analyze LAI by using Landsat 8 and Sentinel-2 time series and the corresponding ground measurements collected in pure Anatolian black pine [Pinus nigra J.F. Arnold ssp. pallasiana (Lamb.) Holmboe] stands within seven-month (from June to December) period. A total of 30 sample plots were selected and seven-month changes of LAI values were determined through hemispherical photography for each sample plot. Remote sensing (reflectance values and vegetation indices obtained from Landsat-8 and Sentinel-2) and topographic (elevation, aspect, and slope) data were used to model the LAI for each month using multiple linear regression (MLR) method. Additionally, the data for all months were combined and modeled. In this case, autoregressive modeling techniques were used to solve the temporal autocorrelation problem. Our study indicated that the models developed from Sentinel-2 give more successful results than Landsat 8 on monthly LAI models. The most successful models were obtained in June by using the reflectance values ( R adj 2  = 0.39, RMSE = 0.3138 m2 m-2), reflectance values-topographic data ( R adj 2  = 0.59, RMSE = 0.3174 m2 m-2), vegetation indices-topographic data ( R adj 2  = 0.82, RMSE = 0.2126 m2 m-2), and reflectance values-vegetation indices-topographic data ( R adj 2  = 0.93, RMSE = 0.1060 m2 m-2). Among the autoregressive modeling techniques, the highest success was obtained with the Landsat 8 OLI using the moving average (2) procedure (R 2 = 0.56). This study is significant that it is the first to analyze the monthly effect on LAI modeling and mapping in pure Anatolian black pine stands using both reflectance values, vegetation indices, and topographic data.

Introduction

Forests, one of the world’s most important terrestrial ecosystems, provide a variety of ecosystem services and contain more than three-quarters of the world’s biodiversity (FAO 2010). In the last quarter, climate change has a serious impact on forest ecosystems (Birdsey and Pan 2011; Xie et al. 2021). Today, global climate change is one of the most important concerns affecting humanity. Climate change is influenced by vegetation on both global and regional scales. Because the leaf area index (LAI) is a crucial component in the global environmental change response process model for vegetation canopy, it is critical to estimate LAI efficiently and precisely in forest ecosystems (Zhang et al. 2002).

One-half of the total green leaf area per unit ground surface area is described as the LAI (Chen and Black 2010). LAI is a significant vegetation structural variable that regulates energy exchanges, water and carbon flows, and net primary productivity in forest ecosystems (Fang et al. 2019; Kang et al. 2021). Also, LAI is a biophysical variable that affects photosynthesis, interception, evapotranspiration, nutrient, and energy balance. Therefore, LAI is an important ecosystem characteristic, and it is closely related to productivity in many terrestrial ecosystems (Atwell et al. 1999; Tuzet et al. 2003). In addition, the total amount of LAI is a significant vegetation parameter that may be used to quantify and model the role of vegetation in a variety of Earth surface processes, including rainfall interception, aboveground biomass and carbon outflow, as well as primary productivity (Meyer et al. 2019).

There exist two ways to calculate LAI: direct and indirect approaches (Chen et al. 1997). LAI produces more accurate findings in the direct approach than in the indirect approach; however, the direct method is exceedingly difficult, time-consuming, and costly (Kamal et al. 2016; Sumnall et al. 2016). Also, LAI is difficult to estimate from ground observations over broad areas and long time periods (Brown et al. 2019).

However, remote sensing data have been utilized to calculate LAI, which is more efficient than field measurements, especially for large forest areas. Because of the capability of remote sensing data to cover a broad area, it reduces field measuring activity (Gray and Song 2012; Pu and Cheng 2015; Günlü et al. 2017; Zhou et al. 2017). In the last two decades, many studies have been conducted to estimate LAI using various satellite images with different spatial and temporal resolutions (Tian et al. 2017; Lu and He 2019; Brede et al. 2020; Brown et al. 2021; Wang et al. 2022). Landsat satellite images, particularly Landsat 7 and 8, have demonstrated the ability to predict forest LAI (Soudani et al. 2006; Brede et al. 2020; Kang et al. 2021). In addition, due to higher spatial and temporal resolution than Landsat satellite images, Sentinel-2 imagery, provided by the European Space Agency (ESA), has recently been utilized to detect various vegetation characteristics such as LAI at local and regional scales (Addabbo et al. 2016; Wang et al. 2022). Furthermore, there are studies in which Landsat and Sentinel-2 images are evaluated together in the estimation of LAI (Ganguly et al. 2012; Korhonen et al. 2017; Meyer et al. 2019). In studies on the estimation of LAI in forest ecosystems, band reflectance (Stenberg et al. 2008; Günlü et al. 2017; Vafaei et al. 2021), vegetation indices (Kodar et al. 2008; Vafaei et al. 2021; Kinane et al. 2021), texture (Zhou et al. 2017), or combined (Wang et al. 2022) data are mainly used as remote sensing data.

When looking through the literature, it was revealed that there are some recent studies to assess LAI and topographic data (elevation, aspect, slope etc.) with remote sensing data such as reflectance values, vegetation indices, and texture features. However, there are a few studies in the literature that combine topographic and remote sensing data to estimate some stand parameters. Wang et al. (2010) produced good results in predicting the stand volume by combining the spectral properties of the ZiYuan 3 (ZY-3) image with topographic information. Another study conducted by Hilbert and Schmullius (2012) using Pol-InSAR data indicated that slope had a considerable impact on tree height estimation. Similar study by Xie et al. (2017) estimated the stand basal area and stand volume by using the reflectance and vegetation indices obtained from the Spot-5 satellite image. In the same study, it was observed that the success of the model increased when topographic features were included in the model.

Many statistical models are used to estimate LAI based on remote sensing data, and they can be classified into two groups: (1) parametric linear models (Günlü et al. 2017; Meyer et al. 2019) and nonlinear models (Zhang and Song 2021), and (2) nonparametric models such as artificial neural networks (Xie et al. 2021), random forest (Vafaei et al. 2021), and support vector machines (Verrelst et al. 2012; Cohrs et al. 2020). Multiple linear regression (MLR) is the most commonly used method in LAI estimation studies (Soudani et al. 2006; Günlü et al. 2017; Guo et al. 2021; Wang et al. 2022). The MLR analysis is often utilized to estimate future values using features of a particular time series or other relevant time-series data (Chatfield 2004). When there is temporal and spatial dependence between the data, autocorrelation problem occurs, and it is not appropriate to use the MLR method. Autoregressive (AR), moving average (MA), and autoregressive moving average (ARMA) models (Thomas and Fiering 1962; Box and Jenkins 1976) were used to solve this problem. The ARMA modeling has become one of the most common methods of time-series modeling in recent years. ARMA is the most widely applied to time-series predicting, and it offers complementary methods to the problem (Zhang 2003; Box et al. 2008).

This study aims to investigate the relationships between LAI (from May to December) and the reflectance values and vegetation indices obtained from Landsat 8 OLI and Sentinel-2, and topographic data (elevation, aspect, and slope). However, it is unknown whether combining topographic data may improve the accuracy or how topographic data can change the LAI prediction accuracy by remote sensing dataset. The relationships between the reflectance values and monthly vegetation indices, as well as the topographical data of the study area were investigated, and the LAI values modeled for these seven months using MLR method. Three models developed for each month to estimate LAI: (1) the model only using reflectance values, (2) the model using both reflectance values and vegetation indices, and finally (3) the model using reflectance values, vegetation indices, and topographic data. Another aim of this study is to evaluate the performance of the Landsat 8 and Sentinel-2 images on LAI modeling.

Materials and methods

Study area

This study was carried out in pure Anatolian black pine [Pinus nigra J.F. Arnold ssp. pallasiana (Lamb.) Holmboe] stands, which are distributed within the borders of Yenice and Üçpınar Forest Planning Units, in Ilgaz Forest Enterprise, Ankara Regional Directorate of Forestry (41° 0′ 36ʺ–41° 2′ 22ʺ N, 33° 42′ 21ʺ–33° 44′ 37ʺ E) (Fig. 1). The average slope of the study area is 42%, which varies from 21 to 91%. Also, the average elevation is 1385 m above sea level, with minimum and maximum values of 1285 and 1585 m, respectively. The climate of study area is mainly temperate climate type. The highest temperature is 36 °C (July), and the lowest temperature is − 27 °C (January). The annual mean temperature is 12.0 °C, whereas the annual mean precipitation is 474 mm.

Fig. 1

Location of the study area

Ground measurements

Using the forest cover type map, a total of 30 sample plots were selected by considering different growth stages, crown closure, site index and topographic conditions (Table 1). LAI was determined by the indirect method in this study. A total of five hemispherical photographs from the north, south, east, west, and center of each sample plot were taken to determine the LAI in each sample plot. At the beginning of the study, it was aimed to determine the LAI data in the study area for each month of 2020. However, due to Covid-19 and unfavorable weather conditions, LAI data could not be obtained for the first five (from January to May 2020) months. Therefore, the hemispherical photographs were taken from each sample area for seven months (from June to December). The photographs were then analyzed by using Hemiview software, and the LAI value for each sample plot was obtained by taking the average. The flowchart for the calculation of the LAI is shown in Fig. 2.

Table 1

Distribution of sample plots with respect to some stand characteristics and topographic data

Number of sample plot	Development stage	Crown closure class	Age class	Site index class	Aspect (°)	Elevation (m)	Slope (%)
1	c	3	4	II	234	1322	26
2	d	2	6	II	322	1316	57
3	d	3	5	II	326	1312	33
4	cd	3	5	II	291	1318	34
5	c	3	4	II	303	1303	52
6	bc	3	3	II	270	1285	39
7	bc	3	3	II	290	1291	29
8	d	2	6	II	326	1402	92
9	d	2	6	II	337	1394	61
10	d	2	6	II	10	1414	45
11	bc	3	2	III	176	1380	37
12	bc	3	3	III	163	1392	22
13	bc	3	2	III	175	1379	33
14	b	3	2	III	120	1428	32
15	b	3	2	III	115	1428	42
16	b	3	2	III	113	1462	60
17	cd	2	5	III	133	1458	37
18	cd	3	5	III	105	1476	51
19	d	2	6	II	171	1585	32
20	c	3	4	II	185	1571	39
21	cd	2	6	II	163	1560	32
22	d	1	6	III	343	1333	38
23	bc	3	2	III	235	1310	42
24	d	2	6	III	238	1366	38
25	d	1	6	III	230	1336	69
26	cd	1	5	III	343	1313	48
27	cd	2	5	III	273	1424	34
28	bc	3	3	III	300	1407	40
29	c	3	3	II	262	1312	47
30	c	3	3	II	277	1293	32

Development stage (cm) a: 0–7.9, b: 8–19.9, c: 20–35.9, d: 36–51.9, e: ≥ 52

Crown closure class (%) 1: 10–40, 2: 41–70, 3: 71–100

Site index class (m) I: 30–35, II: 25–29.99, III: 20–24.99, IV: 15–19.99, V: 10–14.99

Age class (year) 1: 0–20, 2: 21–40, 3: 41–60, 4: 61–80, 5: 81–100, 6: 101–120

Fig. 2

Flowchart for the calculation of the LAI values: a photo-acquisition plan in the sample plots b taking a photograph c analyzing the photographs and calculating the LAI

Remote sensing data

Landsat 8 data

Landsat 8 orthorectified images with cloud cover less than 5% were downloaded from the dedicated web page of United States Geological Survey (USGS—https://earthexplorer.usgs.gov/). A total of seven images (on May 9, July 12, August 29, September 19, October 16, November 17, and December 3, 2020) were used (path/row: 177/31 for June, July, and August; 177/32 for September, October, November, and December). Since all the images in June were cloudy, the image acquired on May 3, 2020, was used for June. Six spectral bands (Band 2-blue, Band 3-green, Band 4-red, Band 5-NIR, Band 6-SWIR-I and Band 7-SWIR-II) with 30-m spatial resolution were used in this study. All bands were extracted according to the study area boundary. Using the radiance rescaling parameters in the sensor metadata file, the radiance values were converted to reflectance using QGIS. Then, buffer zones (with radius of 11.28 for 400 m2, 13.82 for 600 m2, and 15.96 m for 800 m2) were added to the sample plots according to the size of each sample plot. The reflectance value of each sample plot was calculated by taking the average of the reflectance values remaining in the buffer for each sample plot using GIS.

Sentinel-2 data

We utilized freely accessible Sentinel-2 images (Level-2A) with cloud cover less than 17% downloaded from the USGS web page. No cloud removal procedure was applied since there was no significant cloud contamination over the study area. A total of seven images (on June 27, July 12, August 6, September 10, October 3, November 2, 2020, and January 3, 2021) were used. However, since all the images available for December were cloudy, the image (on January 3, 2021) was used as the closest one to December. We used all bands with 10- and 20-m spatial resolutions (Band 2-blue, Band 3-green, Band 4-red, Band 5-vegetation red edge, Band 6-vegetation red edge, Band 7-vegetation red edge, Band 8-NIR, Band 8a-vegetation red edge, Band 11-SWIR, and Band 12-SWIR). All bands were extracted according to the study area boundary. Since Sentinel 2A images directly provide reflectance values, no conversion was necessary. The reflectance value of each sample plot was calculated by taking the average of the reflectance values remaining in the buffer for each sample plot using Geographic Information Systems (GIS). Using the reflectance for the bands generated from Landsat 8 and Sentinel-2 imagery, the vegetation indices in Table 2 were calculated for each sample plot.

Table 2

Vegetation indices used in modeling the monthly LAI

Vegetation indices	Formula	References	Vegetation indices	Formula	References
Albedo	(B + G + R + NIR + SWIR1 + SWIR2)	Lu et al. (2004)	MNDVI (Modified normalized difference vegetation index)	(NIR − SWIR2)/(NIR + SWIR2)	Jurgens (1997)
ARVI (Atmospherically resistant vegetation index)	(NIR) − 2(R) + (G)/(NIR) + 2(R) − (G)	Kaufman and Tanre (1992)	MSAVI (Modified soil adjusted vegetation Index)	(2 × NIR + 1 − sqrt ((2 × NIR + 1)2 − 8 × (NIR − R)))/2	Qi et al. (1994)
BNDVI (Blue-normalized difference vegetation index)	(NIR − B)/(NIR + B)	Wang et al. (2007)	MSI (Moisture stress index)	(SWIR1/NIR)	Hunt and Rock (1989)
CTVI (corrected transformed vegetation index)	(NDVI + 0.5)abs(NDVI + 0.5)xsqrt(abs(NDVI + 0.5))	Perry and Lautenschlager (1984)	ND53	(SWIR1 − R)/(SWIR1 + R)	Lu et al. (2004)
CVI (Chlorophyll vegetation index)	(NIR/ G) × (R/G)	Vincini et al. (2007)	ND73	(SWIR2 − R)/(SWIR2 + R)	Lu et al. (2004)
EVI (Enhanced vegetation index)	2.5 × ((NIR − R)/(NIR + 6 × R − 7.5 × B + 1))	Liu and Huete (1995)	NDVI (Normalized difference vegetation index)	(NIR − R)/(NIR + R)	Rouse et al. (1974)
GEMI (Global environment monitoring index)	(n × (1 − 0.25 × n) − (NIR − 0.125)/(1 − NIR)); n = 0.5	Pinty and Verstraete (1992)	NDWI (Normalized difference water index)	(NIR − SWIR1)/(NIR + SWIR1)	McFeeters (1996)
DVI (Difference vegetation index)	(NIR − R)	Tucker (1980)	NLI (Nonlinear index)	((NIR2) − R)/((NIR2) + R)	Goel and Qin (1994)
GNDVI (Green normalized difference vegetation index)	(NIR − G)/(NIR + G)	Gitelson and Merzlyak (1996)	PSSR (Pigment-specific simple ratio)	(NIR/R)	Blackburn (1998)
GOSAVI (Green optimized soil adjusted vegetation index)	(NIR − G)/(NIR + G + Y); Y = 0.723	Rondeaux et al. (1996)	PVR (Photosynthetic vigour ratio)	(G − R)/(G + R)	Metternicht (2003)
GRNDVI (Green red NDVI)	(NIR − (G + R))/(NIR + (G + R))	Gitelson and Merzlyak (1996)	RVI (Ratio vegetation index)	(NIR/R)	Jordan (1969)
GSAVI (Green soil adjusted vegetation index)	(NIR − G)/(NIR + G + L) × (1.0 + L); L = 0.752	Sripada (2005)	SARVI (Soil and atmospherically resistant vegetation index)	(1.0 + L) × (NIR − (Rr − y × (RB − Rr)))/(NIR—(Rr − y × (RB − Rr)) + L) y = 0.735; Rr = 0.740; L = 0.487; RB = 0.560	Kaufman and Tanre (1992)
GVMI (Global vegetation moisture index)	((NIR + 0.1) − (SWIR2 + 0.02))/((NIR + 0.1) + (SWIR2 + 0.02))	Ceccato et al. (2002)	VIS123	(B + G + R)	Lu et al. (2004)
IPVI (Infrared percentage vegetation index)	(NIR/(NIR + R))	Crippen (1990)	WDRVI (Wide dynamic range vegetation Index)	(0.1 × NIR − R)/(0.1 × NIR + R)	Wang et al. (2007)
MID (Middle infrared wavelengths)	(SWIR1 + SWIR2)	Kaufman and Remer (1994)	WDVI (Weighted difference vegetation index)	(NIR − a × R); a = 0.460	Clevers (1989)

Topographic data

The digital elevation model (DEM) produced from the Alos-Palsar image of the study area was used to extract the topographic data such as aspect, slope, and elevation of the sample plots. The elevation of each sample plot was determined from the DEM data. Additionally, the slope (percentage) and aspect (degree) maps were produced from DEM. Then, the slope and aspect for each sample plot’s were calculated. Descriptive statistical values for the aspect, slope and elevation for sample plots are given in Table 3.

Table 3

Descriptive statistical values for the aspect, slope, and elevation for sample plots

Topographic data	Number of sample plot	Minimum	Maximum	Mean	Standard deviation
Aspect (0)	30	10.30	343.03	227.59	87.1948
Slope (%)	30	21.54	91.58	42.37	14.3251
Elevation (m)	30	1285.00	1585.00	1385.67	84.2711

Modeling process

The response variable was the LAI predicted from the hemispherical photographs, whereas the predictor variables were reflectance values, vegetation indices, and topographic data. Firstly, reflectance values and vegetation indices were processed individually for each month and the associated image to build the models. In the second stage, LAI was modeled using reflectance values–topographic data, vegetation indices–topographic data, and reflectance values–vegetation indices–topographic data independent variable combinations. As a result, five models were produced using different independent variables and their combinations for each month. A total of 70 models were developed for two different sensors (i.e., Landsat 8 and Sentinel-2) and seven months to predict LAI. Therefore, the impacts of considering reflectance values, vegetation indices separately and in combination on the estimation success of LAI values were evaluated. To assess the accuracy of the obtained models, the coefficient of determination (), root-mean-square error (RMSE), and correlation coefficient (r) were calculated. Additionally, LAI images were produced for Sentinel-2 using ESA’s SNAP biophysical products (Weiss and Baret, 2019). Using these LAI images, the LAI values were obtained for 30 sample plots and seven months. Then, paired t test was used to check whether the LAI values obtained from the Sentinel-2 and the ground measurements were comparable.

Multiple linear regression analysis

The multiple linear regression (MLR) algorithm was employed to determine the relationship among LAI and reflectance values, vegetation indices, and topographic data. MLR was also used to develop the monthly LAI estimations with stepwise procedure. IBM SPSS 23 software was used to fit the MLR models. The structure of the MLR model applied in the study is as follows:where LAI is the leaf area index, X1 … X are the variable vectors corresponding to Landsat 8 and Sentinel-2 (reflectance values and vegetation indices) and topographic data (aspect, slope, and elevation), β1 … βn characterize the model coefficients, and finally, ε is the additive bias.

Autoregressive modeling

Multiple measurements in time were carried out to model the LAI in the study. Therefore, the temporally dependent data structure was inherently occurred in dataset, resulting in autocorrelation problem. This data structure not only leads to increase in the variance of the estimated parameters, but also affects the sign and the magnitude. Ultimately, the erroneous predictions can be obtained, and hence model results can be misleading (Whittingham et al. 2006; Zuur et al. 2010). To alleviate this issue, autoregressive (AR), moving average (MA), and autoregressive moving average (ARMA) (Thomas and Fiering 1962; Box and Jenkins 1976), which are time-series models, were used. In the AR(p) model, the Z value is the linear function of the weighted sum of the p-period past values of the series and the random error term. In the MA(q) model, the Z value is the linear function of the backward error terms of series over q-period and its mean. ARMA(p, q) model is a general stationary stochastic process model, and a linear function of past observations and past error terms. These models can be expressed as:where Z is standardized of the observed time series, p and are the order and coefficient of the AR model, q and are the order and coefficient of the MA model, , , and represent the stochastic series of predictions from AR, MA, and ARMA models, respectively. In the study, the LAI values from each month were combined into a dataset. Then, the month variable was used as a dummy variable to incorporate the temporal dependency in the data. For determining the best fit, the autoregressive models having different degree of temporal dependence (i.e., AR(1), AR(2), AR(3), MA(1), MA(2), MA(3), ARMA(1), ARMA(2) and ARMA(3)) were run on the final dataset. The increase in degree of the temporal dependency means that temporal dependency among the monthly datasets was stronger. All analysis were performed in SAS software with PROC MODEL procedure.

Model evaluation criteria

In order to assess the performance of the obtained models, r, R2, , and RMSE metrics were used. For high model success, r, R2, and values were expected to be high, on the other hand, RMSE to be low.where h is residual, X is primary variable, Y is secondary variable, y is observation, y is estimation, y is mean of measured values, n is the number of samples, and p is the number of independent variables in the model. The general flowchart of the method used in the study is shown in Fig. 3.

Fig. 3

Flowchart indicating the overall methodology used in this study

Results and discussion

The LAI values were obtained from ground measurements for each sample plots in the months of June, July, August, September, October, November, and December 2020. The distribution of the observed LAI values was normal according to Kolmogorov–Smirnov normality test (p > 0.05). Also, the LAI values for each month were derived from the SNAP toolbox using Sentinel-2. The paired t test was used to determine whether there was a significant difference between the ground-sampled and SNAP toolbox-derived LAI values, and the result revealed that there was no significant difference in the LAI values for the months (p > 0.05). Descriptive statistical values for monthly calculated LAI values from ground measurements and SNAP toolbox are represented in Table 4.

Table 4

Descriptive statistical values of LAI obtained from ground measurements and SNAP toolbox (Sentinel-2) from June to December 2020

Month	Number of sample plot	Minimum		Maximum		Mean		Standard deviation
		Observed	SNAP	Observed	SNAP	Observed	SNAP	Observed	SNAP
June	30	0.570	0.580	2.160	2.560	1.312	1.539	0.401	0.360
July	30	0.570	0.660	2.320	1.820	1.357	1.308	0.439	0.264
August	30	0.590	0.640	2.060	2.410	1.259	1.422	0.371	0.360
September	30	0.560	0.490	1.910	2.090	1.161	1.210	0.331	0.403
October	30	0.470	0.410	1.410	1.890	0.991	1.118	0.251	0.393
November	30	0.450	0.200	1.300	1.700	0.909	0.932	0.242	0.417
December	30	0.410	0.350	2.410	1.180	0.906	0.826	0.370	0.192

Multiple linear regression (MLR) analysis was used to predict the relationships between the seven-month leaf area index (LAI) value and the reflectance values and vegetation indices obtained from Landsat 8 and Sentinel-2 time series in this study. Also, it was investigated how the model’s success changed when topographic data were included in the MLR models. A total of ten MLR equations were generated for each month in the study. As a result, 70 MLR models were developed for seven months (Table 5). When Table 6 is examined, Sentinel-2, the reflectance values model (except for October and November month) performed better than Landsat 8 models in predicting the LAI. The best model successes were produced using Landsat 8 for August, October, and November; and Sentinel-2 for June, July, September, and December by combining topography data with reflectance values. When the success of the models obtained using vegetation indices was evaluated, the best model success were obtained for five months (from June to October) in Sentinel-2, for two months (November and December) in Landsat 8. When the RMSE values of the models were found to be the same for both satellite images, the model obtained from the Sentinel-2 with a smaller RMSE value can be considered as more successful. The best model successes were obtained from Sentinel-2 for six months (June, August, September, October, November, and December) and Landsat 8 for one month (July) by combining topography data with vegetation indices. For all months, except September, Sentinel-2 outperformed the Landsat 8 in the models combining reflectance values, vegetation indices, and topographic data. The model information for the best model among the LAI models produced for seven months were available in Table 6.

Table 5

Model evaluation criteria of seven-month MLR analysis for estimating LAI

Month	Auxiliary data	Landsat 8 OLI			Sentinel-2
		r	Radj2	RMSE	r	Radj2	RMSE
June	Reflectance	0.36	0.13	0.362	0.62	0.39	0.290
	Reflectance—Topography	0.61	0.37	0.297	0.77	0.59	0.281
	Vegetation indices	0.47	0.22	0.400	0.72	0.52	0.222
	Vegetation indices—Topography	0.70	0.49	0.256	0.91	0.82	0.169
	Reflectance—Vegetation indices—Topography	0.67	0.45	0.245	0.96	0.93	0.064
July	Reflectance	0.42	0.18	0.379	0.55	0.30	0.382
	Reflectance—Topography	0.63	0.40	0.334	0.68	0.46	0.324
	Vegetation indices	0.39	0.15	0.365	0.63	0.40	0.317
	Vegetation indices—Topography	0.75	0.56	0.252	0.72	0.52	0.238
	Reflectance—Vegetation indices—Topography	0.76	0.58	0.225	0.87	0.75	0.192
August	Reflectance	0.45	0.20	0.317	0.55	0.30	0.324
	Reflectance—Topography	0.65	0.42	0.275	0.66	0.43	0.287
	Vegetation indices	0.51	0.26	0.287	0.81	0.65	0.170
	Vegetation indices—Topography	0.70	0.49	0.253	0.88	0.78	0.146
	Reflectance—Vegetation indices—Topography	0.75	0.56	0.230	0.94	0.88	0.087
September	Reflectance	0.37	0.14	0.311	0.48	0.23	0.298
	Reflectance—Topography	0.53	0.28	0.301	0.71	0.50	0.251
	Vegetation indices	0.55	0.30	0.265	0.69	0.47	0.219
	Vegetation indices—Topography	0.78	0.61	0.188	0.79	0.62	0.184
	Reflectance—Vegetation indices—Topography	0.92	0.85	0.111	0.81	0.65	0.174
October	Reflectance	0.57	0.33	0.208	0.49	0.24	0.218
	Reflectance—Topography	0.66	0.43	0.189	0.62	0.38	0.193
	Vegetation indices	0.49	0.24	0.182	0.77	0.60	0.117
	Vegetation indices—Topography	0.62	0.39	0.155	0.82	0.68	0.100
	Reflectance—Vegetation indices—Topography	0.70	0.49	0.151	0.84	0.70	0.090
November	Reflectance	0.52	0.27	0.220	0.43	0.18	0.216
	Reflectance—Topography	0.71	0.51	0.178	0.49	0.24	0.199
	Vegetation indices	0.52	0.27	0.199	0.39	0.15	0.184
	Vegetation indices—Topography	0.75	0.57	0.144	0.85	0.72	0.121
	Reflectance—Vegetation indices—Topography	0.81	0.65	0.127	0.91	0.82	0.090
December	Reflectance	0.47	0.22	0.191	0.54	0.29	0.229
	Reflectance—Topography	0.62	0.38	0.141	0.71	0.51	0.123
	Vegetation indices	0.77	0.60	0.162	0.74	0.55	0.152
	Vegetation indices—Topography	0.76	0.58	0.101	0.85	0.73	0.111
	Reflectance—Vegetation indices—Topography	0.81	0.65	0.085	0.92	0.84	0.079

Table 6

Independent variables and coefficients for the most successful LAI model among months

Independent variable	Coefficient	t value	Sig
Constant	31.107	2.636	0.023
Band 3	414.905	4.027	0.002
Band 5	− 53.693	− 4.728	0.001
Band 6	37.757	5.691	0.000
Band 11	230.205	6.395	0.000
Band 12	21.985	2.640	0.023
ARVI	− 90.012	− 6.322	0.000
VIS123	− 197.085	− 4.283	0.001
ND53	106.133	5.986	0.000
PSSR	− 1.509	− 5.277	0.000
NDWI	81.115	5.591	0.000
MSI	− 41.603	− 4.749	0.001
DVI	155.886	2.402	0.035
CVI	− 1.851	− 2.978	0.013
SARVI	5.786	3.888	0.003
WDVI	− 436.636	− 6.135	0.000
Elevation	− 0.002	− 3.936	0.002
Aspect	− 0.004	− 8.531	0.000
Slope	0.005	2.341	0.039

In the study, the best model successes were obtained by combining reflectance values, vegetation indices, and topographic data together. The most successful MLR model in the reflectance values were obtained using Sentinel-2 ( = 0.39), and with topographic data included, where the model specification coefficient increased ( = 0.59) for June. The highest success with the vegetation indices was obtained for August with the Sentinel-2 ( = 0.65). When vegetation indices and topographic data were combined, the highest success was found for June with the Sentinel-2 ( = 0.82). When reflectance values, vegetation indices and topographic data sets were combined, the most successful MLR model was obtained using Sentinel-2 for June ( = 0.93, RMSE = 0.064). The scatterplots of LAI model with remote sensing and topographic data are shown in Figs. 4 and 5. Also, the variations of the coefficient of determination with respect to months according to Landsat 8 and Sentinel-2 are shown in Fig. 6.

Fig. 4

Predicted LAI values against observed LAI values for the MLR models from Landsat 8

Fig. 5

Predicted LAI values against observed LAI values for the MLR models from Sentinel-2

Fig. 6

Monthly R2 change; a reflectance values b reflectance values–topographic data c vegetation indices d vegetation indices–topographic data e reflectance values–vegetation indices–topographic data

In this study, whole LAI datasets without separating them by months were modeled with independent variables obtained from Landsat 8 (cf. Table 7), Sentinel-2 (cf. Table 8), and topographic data using MLR method. It was determined that there was an autocorrelation problem in the MLR models obtained for Landsat 8 (DW = 1.58) and Sentinel-2 (DW = 1.42) (cf. Tables 7, 8). In addition, a positive autocorrelation was observed in the MLR model obtained with both satellite images (Landsat 8 Pr < DW = 0.0003, Sentinel-2 Pr < DW =  < 0.0001).

Table 7

Parameters of the “best fit” MLR model of LAI based on the reflectance values, vegetation indices derived from Landsat 8 and topographic data for seven months

Parameter	Coefficient	Approx SE	t value	Pr >|t|
Constant	2.235937	0.688500	3.25	0.0014
Band 2	− 5.567590	2.729300	− 2.04	0.0427
MSI	− 4.153290	1.285000	− 3.23	0.0014
EVI	5.484900	1.978600	2.77	0.0061
DVI	− 10.856400	3.143200	− 3.45	0.0007
ARVI	7.350215	1.810700	4.06	< .0001
ND53	11.555490	2.243100	5.15	< .0001
ND73	− 5.632200	1.455800	− 3.87	0.0001
CVI	− 0.785410	0.219700	− 3.57	0.0004
Elevation	− 0.000800	0.000297	− 2.68	0.0080
Aspect	− 0.001740	0.000299	− 5.82	< .0001
Slope	− 0.009390	0.001450	− 6.48	< .0001
Month	− 0.114940	0.022300	− 5.15	< .0001

DW, Durbin–Watson statistics; Pr < DW, positive autocorrelation; Pr > DW, negative autocorrelation

Table 8

Parameters of the “best fit” MLR model of LAI based on the reflectance values, vegetation indices derived from Sentinel-2 and topographic data for seven months

Parameter	Coefficient	Approx SE	t value	Pr >|t|
Constant	1.890271	0.576400	3.28	0.0012
Band 5	3.919959	1.096500	3.57	0.0004
ND73	2.575112	0.759100	3.39	0.0008
MNDVI	3.393199	0.875400	3.88	0.0001
GRNDVI	− 2.657380	0.759400	− 3.50	0.0006
Elevation	− 0.000890	0.000310	− 2.89	0.0043
Aspect	− 0.001780	0.000303	− 5.87	< .0001
Slope	− 0.008890	0.001490	− 5.99	< .0001
Month	− 0.093050	0.011600	− 8.03	< .0001

DW, Durbin–Watson statistics; Pr < DW, positive autocorrelation; Pr > DW, negative autocorrelation

The positive autocorrelation problem was handled by applying AR (p), MA (q) and ARMA (p, q) autoregressive modeling procedures. Three procedures were applied for AR (p), MA (q), and ARMA (p, q) techniques, and all Durbin–Watson statistical values were approximately 2 with the application of autoregressive modeling (cf. Table 9).

Table 9

Durbin–Watson statistics values of applied autoregressive modeling results

Satellite	Procedure	DW	Pr < DW	Pr > DW
Landsat 8 OLI	AR(1)	2.06	0.6009	0.3991
	AR(2)	1.95	0.2840	0.7160
	AR(3)	2.01	0.4436	0.5564
	MA(1)	1.92	0.2085	0.7915
	MA(2)	2.01	0.4688	0.5312
	MA(3)	2.00	0.4393	0.5607
	ARMA(1, 1)	2.03	0.5085	0.4915
	ARMA(2, 2)	2.02	0.4703	0.5297
	ARMA(3, 3)	2.00	0.4566	0.5434
Sentinel-2	AR(1)	2.07	0.6387	0.3613
	AR(2)	1.98	0.3859	0.6141
	AR(3)	2.00	0.4455	0.5545
	MA(1)	1.90	0.1898	0.8102
	MA(2)	1.98	0.4016	0.5984
	MA(3)	2.00	0.4556	0.5444
	ARMA(1, 1)	2.03	0.5004	0.4996
	ARMA(2, 2)	2.00	0.4229	0.5771
	ARMA(3, 3)	2.01	0.4286	0.5714

DW, Durbin–Watson statistics; Pr < DW, positive autocorrelation; Pr > DW, negative autocorrelation

The highest success was obtained with the MA(2) technique, and all independent variables included the MA(2) model were significant (cf. Table 10). The fact that the month parameter was significant in the model indicates that there was a significant variability according to the months (p < 0.05). Negative coefficient of the month variable indicated that LAI decreased in the following months, and there was an inverse relationship with respect to time (coefficient of the month variable = − 0.08876). LAI surface maps for pure Anatolian black pine stands and predicted–observed graphs were represented by applying the model obtained for each month using the MA(2) technique in Figs. 7 and 8.

Table 10

Parameters of the “best fit” MA(2) model of LAI based on Landsat 8 and topographic data for seven months

Parameter	Coefficient	Approx SE	t value	Pr >|t|
Constant	3.175782	0.588900	5.39	< .0001
MSI	− 3.467960	1.232000	− 2.81	0.0054
DVI	− 4.003740	1.509300	− 2.65	0.0086
ARVI	3.817366	1.379600	2.77	0.0062
ND53	7.894709	2.253200	3.50	0.0006
ND73	− 3.042050	1.554600	− 1.96	0.0518
CVI	− 0.381480	0.211500	− 1.80	0.0728
Elevation	− 0.000660	0.000337	− 1.94	0.0534
Aspect	− 0.001460	0.000302	− 4.83	< .0001
Slope	− 0.008330	0.001370	− 6.08	< .0001
Month	− 0.088760	0.024500	− 3.63	0.0004
ρ1	− 0.293580	0.070700	− 4.15	< .0001
ρ2	− 0.285030	0.073500	− 3.88	0.0001
Method	DW	Pr < DW	Pr > DW	Radj2	r	RMSE
MA(2)	2.01	0.4688	0.5312	0.56	0.77	0.263

DW, Durbin–Watson statistics; Pr < DW, positive autocorrelation; Pr > DW, negative autocorrelation; 1,2 the temporal autocorrelation parameters

Fig. 7

Surface maps of the monthly LAI values obtained using the autoregressive model based on MA(2) procedure with Landsat 8 and topographic data

Fig. 8

Scatterplots for monthly LAI values obtained using the autoregressive model based on MA(2) procedure with Landsat 8 and topographic data

Discussion

This study aimed to establish models to predict monthly LAI in pure Anatolian black pine stands. We evaluated the MLR models between ground-based monthly LAI and reflectance values and vegetation indices, and topographic data. Five models for each month were established to estimate LAI, involving model only using reflectance values, model using both reflectance values and topographic data, model only using vegetation indices, model using both vegetation indices and topographic data, and finally, model combined using reflectance values, vegetation indices and topographic data. A total of 70 MLR models were developed, including 35 Landsat 8 and 35 Sentinel-2. The study indicated that the models (i.e., 27 MLR model) developed from Sentinel-2 give more successful results than the ones from Landsat 8. Similar to the findings of this study, Meyer et al. (2019) found that the Sentinel-2 estimates LAI with comparable accuracy to Landsat 8. Astola et al. (2019) also reported the superior performance of Sentinel-2 for stand parameter as compared to Landsat 8. However, in a study conducted by Korhonen et al. (2017) for the modeling of LAI, it was seen that the model successes obtained from Landsat 8 and Sentinel-2 were almost comparable. The success of the models developed by using vegetation indices in the Sentinel-2 (except for November) was more successful than the models in which the reflectance values were used. Models with reflectance values as independent variables in July and October and the models with vegetation indices as independent variables in the other months produced better results for Landsat 8. In November, the model successes of both models were found to be the same. This finding agrees with the results of Günlü et al. (2017) which indicated that LAI was successfully estimated in August month using MLR approach when the vegetation indices generated from Worlview-2 and Aster satellite image were utilized. However, better model results (R2 = 0.65, RMSE = 0.170 m2 m−2) were found using vegetation indices from Sentinel-2 in August compared to the study of Günlü et al. (2017) (R2 = 0.45, RMSE = 0.2932 m2 m−2 for WorldView-2; and R2 = 0.41, RMSE = 0.3245 m2 m−2 for Aster). Our findings on the performance of the Landsat 8 for LAI prediction differ than those of Dube et al (2019), while they reported the superiority of Landsat 8 vegetation indices and reflectance values for the model prediction (R2 = 0.72, RMSE = 0.60 m2 m−2) LAI in May month; in our study, the best model success (R2 = 0.60, RMSE = 0.162 m2 m−2 was achieved in December month (cf. Table 6 and Fig. 6). Chrysafis et al. (2020) presented a study that showed similar results to our study. They founded that vegetation indices generated from Sentinel-2 produced the most successful LAI estimations with R2 = 0.854. In the same study, the NDVI was provided as the most important variable for the LAI estimation. On the other hand, our study demonstrated that NDVI variable was not included as an independent variable in the models. Cohrs et al. (2020) found that simple ratio indices were the most important variable for LAI prediction. Their study indicated good performance (R2 = 0.81) at predicting the LAI for loblolly pine stands. This is similar to the findings of Yu et al. (2019), which respectively indicated that LAI was successfully estimated using an MLR when RVI (R2 = 0.63), RVI (R2 = 0.63), NMDI (R2 = 0.75), DVI (R2 = 0.79), and WDRVI (R2 = 0.7439) vegetation indices were used. Wang et al. (2022) reported that vegetation indices (i.e., simple ratio) outperformed the best predictive performance for LAI prediction. The success levels of the models have generally improved in which the vegetation indices are included as an independent variable according to the studies conducted to model the LAI.

In addition, we assessed the combined use of reflectance values, vegetation indices and topographic data for monthly LAI modeling. Our results indicate that incorporating topographic data improves models depending on Landsat 8 and Sentinel-2 using both reflectance values and vegetation indices. In other words, when the reflectance values and vegetation indices were used to estimate monthly models, the addition of topographic data considerably increases the accuracy of the prediction (cf. Table 5). Topographic data played a significant role in predicting LAI in our study. There have been few studies looking into the role of topographic data in LAI and other stand parameters prediction; however, similar studies have improved results by combining topographic data with remote sensing variables such as reflectance values, vegetation indices, and texture features (Xie et al. 2017; Jiang et al. 2020; Chen et al. 2021; Bhattarai et al. 2022). In a study conducted by Bhattarai et al. (2022), the relationships between the variables generated from Sentinel-1, Sentinel-2 and topographic data, and LAI values were investigated in mixed spruce-fir forest area. According to the results, it was seen that the model successes increased by adding topographic data (especially elevation data) to the variables obtained from the Sentinel-2 satellite image, as in our study. They reported that the R2 value of 0.78 is comparatively lower than the R2 achieved in this study (i.e., 0.93). In addition, the RMSE found in this study (RMSE = 0.064 m2 m−2) is better than they reported (RMSE = 0.95 m2 m−2). Xie et al. (2017) utilized a combination of Spot-5 satellite image variables and topographic data to achieve an overall accuracy R2 of 0.658–0.820 and RMSE of 24.457–33.694 m3 ha−1 for modeling stand volume. The factors were ranked according to their importance, and the results revealed that topographic data such as elevation played a significant influence in predicting some stand parameters. In general, the inclusion of topographic data in the final models (Landsat 8 and Sentinel-2) improved the R2 and RMSE for LAI models in all months (cf. Table 5 and Fig. 6). These are in line with the studies of Luther et al. (2019) and Bhattarai et al. (2022), which emphasized the importance of topographic data for improved LAI prediction performance and accuracy. Another study carried out by Moradi et al. (2021) demonstrated that the highest accuracy (R2 = 0.781, RMSE = 20%) in estimating LAI was obtained using the combination with reflectance values, vegetation indices, and topographic data.

Also, we evaluated the relationships between the LAI values and, the reflectance values and vegetation indices variables obtained from both satellite images for seven months, and topographic data by MLR. Autocorrelation problem occurred due to the use of temporal data in the modeling processing. To overcome this problem, AR, MA, and ARMA models, which are time-series models, were used. As a result, the autocorrelation problem was solved, and model success was obtained for both satellite images. When the model successes are compared, the model produced from the Landsat 8 (R2 = 0.56) outperforms the Sentinel-2 (R2 = 0.51), according to the model successes evaluated monthly. In the models for both images, topographic data are included as an independent variable. In addition, it was observed that band 2 and band 5 variables were included in the model developed with the Landsat 8 and Sentinel-2, respectively. However, the ND73 vegetation indices were the common variable in both models excluding topographic data (cf. Tables 7, 8). While the independent variable band 5 was included in the model developed for Sentinel-2 in our study, it was seen that it was also included as an independent variable in the model developed in Sibanda et al. (2021). As contrary to our results, Lee et al. (2004) investigated that the relationships between LAI and reflectance values are generated from multispectral and hyperspectral data. Their study demonstrated that the band spectral values of Landsat 7 were found to be important in the estimation of LAI. Tillack et al. (2014) observed seasonal relationships between LAI predictions using vegetation indices from a RapidEye image. The vegetation indices that showed the best relationships changed according to the seasons. However, when all seasons were added, the NDVI overperformed.

We also compare the ground measurements LAI values with SNAP toolbox-derived LAI values (Sentinel-2). The validity of the in situ LAI was tested with paired t test using LAI values derived from SNAP for 30 sample plots. According to the test results, there were non-significant differences between ground measurement LAI and SNAP LAI values (p > 0.05). Therefore, it was decided that the use of SNAP toolbox LAI values could be used for LAI prediction for pure Anatolian black pine stands.

Many researchers have predicted LAI using a single satellite image (Pu and Cheng 2015; Günlü et al. 2017). However, for accurate estimation and mapping of seasonal LAI, it is important to monitor and quantitatively analyze the structures and functions of forest ecosystems with time-series data (Chen et al. 2002). When the research in this area is evaluated, it can be concluded that the models developed with the temporal texture features (different orientation and window size) obtained from satellite images provide better results in estimating the LAI (Zhou et al. 2014; Pu and Landry 2019). Zhou et al. (2014) indicated that combining vegetation indices with textures improved the LAI estimate accuracy (R2 value from 0.68 to 0.84). Pu and Landry (2019) stated that textural features were more important than band spectral features and vegetation indices generated from pléiades image in predicting and mapping LAI. As a result, it should be highlighted that using texture features extracted from different satellite images in combination with vegetation indices and topographic data in future studies on LAI estimation in forest ecosystems will improve the model success.

Conclusion

In this study, we aimed to predict the potential of monthly LAI by MLR analysis using time series of Landsat 8 and Sentinel-2. Also, we tried to determine how the model successes changed by adding topographic data (elevation, slope, and aspect) to the models produced with the reflectance values and vegetation indices. Additionally, another aim of this study was to evaluate the performance of the Landsat 8 and Sentinel-2 on LAI estimation. The results showed that the Sentinel-2 gave better results than the Landsat 8 in estimating the monthly LAI. However, it has been demonstrated that the success levels of models produced by combining topographic data with variables derived from satellite images produce better outcomes than models developed only from satellite images. In our study, the best result in predicting the LAI was obtained for the month of June using a model combined use of reflectance and vegetation indices values from the Sentinel-2 as well as topographic data. Also, it was seen that the use of topographic data improved the success of models developed for both monthly and seven-month LAI estimation. The models developed for each month in this study will provide a basis for future studies. It may be useful to use these models in similar forest ecosystems. In addition to, the LAI values obtained by ground measurements and the LAI values obtained from SNAP toolbox (Sentinel-2) were compared with paired t test and, there were non-significant differences between ground measurement LAI and SNAP LAI values. Therefore, it has been concluded that the LAI values obtained from SNAP can be used in pure Anatolian black pine stands. For this reason, it would be useful to carry out these studies in different forest ecosystems consisting of different tree species and to compare the results obtained with the results obtained from SNAP (Sentinel-2). When the literature is examined, there are limited studies that model the temporal change of LAI by evaluating the remote sensing and topographic data together. Therefore, we suggested that future studies assess the different modeling techniques (random forest, deep learning, support vector machines, etc.) and different satellite images (Sentinel-1, Sentinel-3, SAR L band, and P band) with topographic data for LAI modeling.