2.1.

Experimental Locations and Workflow

The study area of this research is the rubber plantation of Indonesian Rubber Research Institutes (IRRI) in Sembawa, South Sumatra, Indonesia. The fields have been infected by Pestalotiopsis diseases since 2016 with varying degrees of canopy infection severity. Three fields with different rubber tree clones, namely BPM 24, RRIC 100, and GT 1, and aged 10 to 11 years were selected as experimental areas in this research. Figure 1 and Table 1 present the location and digital elevation model image of the study area.

Fig. 1

Location of the study area.

Table 1

Extent of each field for three rubber tree clones.

The workflow of the proposed work is shown in Fig. 2. The method starts with data acquisition, pre-processing, NDVI calculation, imputation and smoothing, analysis of phenological changes, prediction of vegetation indices, and evaluation.

Fig. 2

General workflow of the proposed method: (1) data acquisition, (2) pre-processing, (3) NDVI calculation, (4) imputation and smoothing, (5) analysis of phenological changes, (6) prediction of vegetation indices, and (7) evaluation.

2.2.

Data Acquisition

In this study, we obtained a 4-year time series data from the Sentinel-2 level 2A product of the study areas from January 1, 2019, to December 31, 2022. The 291 imagery data were acquired from the Google Earth Engine (GEE) platform. Sentinel-2 provides reflectance data for 13 bands and has a 5-day temporal resolution. For this study, we used four spectral bands from a top-of-atmosphere product with a spatial resolution of 10 m: B02 (blue), B03 (green), B04 (red), and B08 (infrared). We also collected UAV multispectral data and digital hemispherical data during disease outbreak as the ground truth data for disease severity at two time periods, between January 18–19, and June 15–16, 2022. The use of UAV data enables a broader understanding of the affected area, providing increased resolution and enhanced accuracy through hemispherical data. These data were used by experts to assess the health of vegetation and determine the severity of Pestalotiopsis disease. Daily climate data were obtained from local weather stations in the rubber plantation and the Indonesian Central Bureau of Statistics and the Indonesian Meteorology, Climatology, and Geophysical Agency from 2019 to 2022.²⁶

2.3.

Pre-Processing Sentinel-2 Imagery

Cloud cover is a common problem in evaluating time-series data using optical imaging. Cloud removal is vital for the subsequent data analysis and to increase the integrity and availability of the remote sensing data. There are three approaches of cloud removal, i.e., spatial information-based method, temporal information-based methods, and multi-source auxiliary information-based method.²⁷ Many studies on cloud removal have been conducted.^28–31

To eliminate cloud cover in this experiment, we used s2cloudless, a machine learning-based cloud identification technology developed by Sinergise’s EO Research team. This technique leverages the light GBM algorithm, which was trained on a massive dataset of cloudy and non-clouded samples from around the world.³² Before applying the gradient boost-based technique, all bands are upsampled to 10 m resolution using bilinear interpolation. The algorithm employs Sentinel-2 Level-1C TOA reflectance values from the following 10 bands as input: B01, B02, B04, B05, B08, B8A, B09, B10, B11, B12; and the output is a cloud probability map. In this study, we used a threshold value of 0.4 to reduce cloud omission errors.

2.4.

Calculation of NDVI

To explore the rubber trees phenology stages and define vegetation behavior during Pestalotiopsis infestation, we studied the NDVI-time series data acquired from Sentinel-2 datasets. The time series graph was created using the NDVI from random sites in the three rubber fields. The NDVI is determined using the equation as follows:

2.5.

Imputation and Smoothing NDVI Time Series

The missing information generated by cloud causes spatial and temporal gaps in satellite earth observation data, as well as biases in subsequent image processing and application. Various methods for dealing with missing data have been proposed. Methods are classified into four types: spatial-based methods, spectral-based methods, temporal-based methods, and hybrid methods.³⁴ In this study, we used a temporal-based method or “gap filling” technique. Recent temporal techniques, such as data fusion,^35,36 pixel unmixing,³⁷ data interpolation,³⁸ or best-pixel selection,^39,40 have been examined for large-area land cover mapping. To obtain the gap-filled satellite image, we used the linear interpolation method in this work. The linear interpolation function equation is as follows:⁴¹

The noise of the signals received by the satellite sensors can occur due to cloud cover, atmospheric disturbance, sensor failure, and so on.⁴² The time series smoothing method can be used to reduce noise and assess plant phenological trends.⁴² Smoothing methods were investigated by Atkinson et al. in reconstructing consecutive time series of vegetation indices to estimate vegetation phenology.⁴³ There were three types of smoothing methods:⁴⁴ temporal methods, frequency methods, and hybrid methods. In this study, we used temporal-based method, a moving average (MA) filter for smoothing NDVI time series data. The MA filter reconstructs NDVI time series by replacing the central point with the average of all points in an equation-based moving window of fixed size. The advantage of this strategy is that the seasonal peak of the curve is preserved.

2.6.

Analysis of Phenological Changes

The analysis of phenological changes of NDVI was investigated using harmonic analysis. The harmonic analysis of NDVI time series aims to characterize the spectral behavior of NDVI by decomposing the time series into a sequence of basic sinusoids and a set of sine or cosine terms with unique coefficients and intercept (harmonic terms).⁴⁵ Time-lag cross-correlation analysis was used to understand the relationship and potential lag (lead effects) between NDVI and climate variables. Cross-correlation measures the similarity between two time series as a function of the time lag applied to one of them. Time lag indicates the time delay between precipitation and temperature changes and their impact on NDVI.

2.7.

Prediction of Vegetation Indices

In this study, we analyzed the relationship between NDVI variations and climatic variables to understand the impact of the Pestalotiopsis diseases on vegetation dynamics in a rubber plantation using RFR. RFR was first proposed by Breiman (1996) as an ensemble learning method that can be used in both classification and regression task.⁴⁶ RFR is a non-linear ensemble technique that focuses solely on decision trees for classification or regression. The concept behind RF for time series is to replace the standard bootstrap with a dependent bootstrap in order to subsample time series during the tree construction phase in order to account for time dependence.⁴⁶

2.8.

Evaluation

To validate the performance of the model, multi-step-ahead monthly NDVI forecasts were simulated to walk-forward validation. The walk-forward validation was utilized to improve the model’s precision. The dataset time series was divided into training, validation, and testing sets, with 90% of the dataset used for training and validation, and 10% retrained for testing. Means absolute percentage error (MAPE) performance evaluation was conducted to validate the accuracy of the proposed model. The following expression is the MAPE equation:

3.1.

Environmental Factors Affecting Pestalotiopsis Disease

The incidence and severity of Pestalotiopsis sp. disease are intricately correlated with various environmental factors that influence the fungus’s growth, spread, and impact on plant populations.⁴⁷ Pestalotiopsis sp. prefers specific temperatures for growth and reproduction. Extreme temperatures cause stress in plants, making them more susceptible to infections.⁴⁸ Temperature fluctuations can also affect the fungus’s growth rate and development, as well as the host plants’. High humidity can create an ideal microclimate for fungal growth. Humidity, along with temperature, affects spore production and dispersal. Pestalotiopsis sp. thrives in humid environments, as do many other fungi.

Environmental stressors, such as intense rainfall or flooding, have the potential to accelerate disease spread by creating favorable conditions for spore dispersal and germination.^48,49 In contrast, prolonged periods of drought have the potential to diminish the vitality of plants and render them more vulnerable to infections.^50,51

The vegetation index is an effective tool for assessing plant health in relation to Pestalotiopsis disease. A declining vegetation index may indicate stress on plant populations, which could be caused by Pestalotiopsis spp. infection. Monitoring this index enables the early detection of potential outbreaks. A thorough understanding of the interactions between temperature, precipitation, vegetation index, and other environmental factors is required for predicting, preventing, and managing Pestalotiopsis disease in rubber plantations. The analysis of the relationships between these variables allows for the identification of patterns and correlations, which aids in the prediction of Pestalotiopsis sp. outbreaks, the implementation of preventive measures, and the development of disease management strategies in affected ecosystems.⁴⁷

3.2.

Phenological Changes Based on Harmonic Analysis of NDVI Time

In this study, we observed three rubber clones (GT 1, BPM24, and RRIC 100) aged 10 to 11 years to study the NDVI trend. Figures 3(a)–3(c) show a harmonic model of NDVI time series from 2019 to 2022 for three rubber clones. We found that the overall trend of the NDVI decreased significantly from 2019 to 2022 for all clones.

Fig. 3

Harmonic model of NDVI time series over four years (from 2019 to 2022) for three clones, (a) GT1, (b) BPM24, and (c) RRIC 100. The red curve illustrates the fitted model, while the blue points denote the actual NDVI data.

Pestalotiopsis was reported in rubber plantations in North Sumatra and South Sumatra since 2016.¹ Figure 4 depicts defoliation and refoliation trends of each clone. According to previous research phenology stage of defoliation and refoliation period for rubber wintering occurs in the dry season, which in Sumatra is between July and September¹⁶ (as shown in Fig. 4 in gray area). Defoliation takes place over a 4-week span from July to August, signifying the conclusion of the growing season. Simultaneously, it marks the onset of the subsequent growing season, transitioning into the refoliation period spanning August to September. During 2019 to 2022, we found that defoliation and refoliation events shifted due to Pestalotiopsis disease. Early wintering was reported in all clones. The first defoliation occurred in March to July 2019, followed by secondary leaf fall in November and January 2020. In 2020 and 2021, the defoliation occurred repeatedly or extensively in all clones, indicating physiological stress, nutrient depletion, and a compromised defense system. These factors create favorable conditions for pests and pathogens to thrive, leading to increased vulnerability and the potential for further damage to the trees.

Fig. 4

Defoliation and refoliation trends of several rubber clones (BPM 24, GT 1, and RRIC 100) with Pestalotiopsis infection in 2019 to 2022. The gray area denotes the annual wintering period in Sumatra (between July and September).

In 2022, the NDVI gradually decreased and defoliation significantly and occurred earlier in March. All clones suffered from extensive defoliation between March and July in 2022. According to weather statistics data in Fig. 5, in 2022, precipitation increased and caused a decrease in minimum temperature compared to 2019 to 2021. The increase in precipitation can create a cascade of unfavorable conditions for rubber trees, including the spread of fungal diseases, waterlogged soil, nutrient deficiencies, disruption of photosynthesis, pest infestations, and physiological stress. These factors collectively contribute to the extensive defoliation observed in rubber trees during periods of heavy rainfall and this can be reflected in temporal NDVI profiles with shifts in the timing of NDVI peaks and valleys.

Fig. 5

Temperature (Tmin and Tmax) and precipitation trends during Pestalotiopsis outbreak from 2019 to 2022.

The phenology of rubber trees is affected by the period of defoliation. A longer defoliation period would consume more conserved carbohydrates, reducing available reserves for disease defense and latex production.²⁴ This condition is associated with significant defoliation and decrease in rubber yield. The shifts in the starting and ending periods of the annual dry season and wintering are likely the result of a complex interplay between Pestalotiopsis disease and broader climate change impacts.

The proliferation of Pestalotiopsis leaf disease has caused extensive defoliation. This symptom indicates the presence of Pestalotiopsis disease in all rubber clones. To prove the occurrence of disease in rubber plantations, in 2022 we collected UAV multispectral data and digital hemispherical data as the ground truth data for disease severity at two time periods, between January 18–19, and June 15–16, 2022. Pathologist experts from IRRI also conducted research to determine the severity of Pestalotiopsis disease in rubber plantations. Figure 6 shows the vegetation index maps of disease severity at different fields (BPM 24, RRIC 100, and GT 1) on January 18, 2022. The colors green, yellow, and red on the map represent the values of NDVI which are high, moderate, and low, respectively. Higher values indicate healthy plants and lower values indicate a plant infested with pests, damaged plants, or unhealthy plants.

Fig. 6

Time-lag cross correlation between NDVI and climate variables (precipitation and temperature).

3.3.

Analysis and Prediction of NDVI Based on Climate Variable Time Series Using Random Forest Algorithm

The values of NDVI, temperatures, humidity, precipitation, radiation, and wind speed were used to assess phenological changes in response to disease outbreaks and climate changes for three rubber clones. We applied the augmented Dickey–Fuller test (ADF test) to test the stationary time series. If the $p$-value is >0.05, then the data series is non-stationary. Table 2 shows that the $p$-values of NDVI for all clones (GT1, BPM24, RRIC), precipitation, maximum temperature, humidity, and radiation are greater than the significance level of 0.05, and the ADF statistic is higher than any of the critical values. Therefore, the sequence of NDVI and climate data were non-stationary.

Table 2

ADF test to test the stationary of time series.

We applied time-lag cross-correlation analysis to understand the relationship and potential lag (lead effects) between NDVI and climate variables. We found a time lag which indicates the time delay between precipitation and temperature changes and their impact on NDVI. Figure 7 shows the lag value associated with the peak of the correlation function of precipitation, temperature, and NDVI. Results show that vegetation has a 1- to 3-month lag response to precipitation and temperature anomalies. It means the precipitation and temperatures highly correlate with a time shift of 1 to 3 months. A negative peak cross-correlation between NDVI and precipitation and a positive peak cross-correlation between NDVI and temperature anomaly correlation are crucial for assessing the impact of changing climate patterns. Anomaly correlation between NDVI and climate variables indicates how variations in climate patterns affect vegetation growth and disease prevalence. The time lag between climatic conditions and changes in NDVI reflects the health conditions of vegetation caused by Pestalotiopsis disease. The significant reduction in NDVI values resulting from disease outbreaks and climatic conditions provides evidence for the development of Pestalotiopsis.

Fig. 7

Time-lag cross correlation between NDVI and climate variables (precipitation and temperature) for (a) GT, (b) BPM24, and (c) RRIC 100.

We applied the lagging time series to the features for predicting the NDVI changes using RFR based on time-lag cross-correlation of temperature and precipitation. Table 3 shows the features of the prediction model.

Table 3

List of features for prediction model.

In order to evaluate the performance of the predictive model, we utilized time series cross-validation, a process that involves conducting multiple validation iterations. The prediction model was built using data from 2019 to 2023. The entire time period was divided into training and testing sets. The expanding window walk-forward validation method, specifically a fivefold time series cross-validation, was used in our evaluation, as shown in Fig. 8. As the cross-validation process progressed, the training data expanded to encompass all data that had preceded it, while the size of the testing data remained constant. The objective of walk-forward validation was to enhance the accuracy of the model in predicting the initiation of defoliation, with forecasting intervals varying between 1 and 5 months.

Fig. 8

Illustration of time series cross-validation.

Before starting the model’s training and testing phases with the RFR model, we used randomized search and grid search hyperparameter tuning. Randomized search is used as the first step in hyperparameter optimization. It aids in quickly gaining an understanding of the hyperparameter space and identifying promising areas for further investigation through grid search. The training set was subjected to grid-search time-series cross-validation in order to identify optimal hyperparameters from the search space and to ensure unbiased tuning. We ran a randomized search for each fold in the time series cross-validation. This experimental setup was used for three different rubber clones. Table 4 shows the description and optimal hyperparameters used in the RFR model for each clone.

Table 4

The description and optimal hyperparameter used in the RF regression model.

Figures 9(a)–9(c) show the predicted NDVI results for each clone, with time horizons ranging from 1 to 5 months. The performance of the prediction model is measured by the mean absolute percentage error (MAPE) which yields the following results for each rubber clone: 8.18% for RRIC 100, 6.07% for BPM24, and 5.96% for GT1. The model had accuracies of 90.80%, 93.65%, and 93.35% for RRIC 100, BPM24, and GT1, respectively.

Fig. 9

The predicted result of NDVI with horizons of 1 to 5 months for (a) RRIC 100, (b) BPM24, and (c) GT1 clones.

Figure 10 depicts the important variables in the correlation between NDVI and climate variables for individual rubber clones using RFR model. Incorporating time lag effects of vegetation and climate variables is especially useful in predictive analysis, particularly when vegetation responses to changing conditions follow temporal patterns. Our findings show that time lag effects recognize that the current NDVI value is influenced not only by the immediate past value but also by values from specific time intervals prior to the current observation. In this study, the time-lag effects of NDVI (lag_ndvi1, lag_ndvi2, rolling_mean_ndvi) and climate variables, such as precipitation (lag_pre2, rolling_mean_pre), temperature (lag_diff temp1), radiation (lag_rad1, lag_rad2), humidity (lag_hum2), and wind speed (wind_speed, lag_win1, lag_win2) influenced the prediction of NDVI values in time series data. Significantly, the 1-month interval preceding NDVI features (lag_ndvi1) emerges as the most important NDVI predictor for all rubber clones.

Fig. 10

Importance of variables in relationship between NDVI and climate variables.

Precipitation and temperature have a much greater impact on the phenological pattern of all rubber clones than humidity, wind speed, and radiation. Important characteristics for NDVI value prediction include the significance of time-lag effects, particularly precipitation (lag_pre2), the variance in maximum and minimum temperatures (diff_temp, lag_diff_temp2), and the ratio of precipitation to temperature (pre_temp, lag_pre_temp2). The ratio of precipitation to temperature offers valuable insights into the equilibrium between temperature and water availability, a critical determinant in evaluating potential stress on vegetation. A disparity in temperature and precipitation may indicate water stress, which has negative impacts on plant development and overall ecosystem health. Furthermore, the combination of increased precipitation and lower temperatures creates an environment conducive to the growth, spread, and establishment of pathogenic microorganisms, particularly fungi. This phenomenon contributes to an increased prevalence of diseases among rubber trees. It is concluded that the temporal lag effects of climate variables and the vegetation index are significant in the prediction of forthcoming vegetation dynamics.

The results show the random-forest-based classification was highly effective in modeling the relationship between NDVI and climate variables. We found a strong correlation between NDVI and climate variables. This correlation indicates that these variables can be used to determine the vegetation cover and to assess phenological changes in response to disease outbreaks. The analysis indicated that the outbreak of Pestalotiopsis disease and the weather pattern have been confirmed to influence the pattern and timing of defoliation and refoliation of rubber trees (wintering period). We conclude that the analysis of the phenological pattern and the proposed prediction model can be used to forecast the early wintering period of rubber trees to prevent the spread of disease and minimize latex yield loss.