Spatiotemporal trends and machine learning-based prediction of temperature variability during the T. Aman rice-growing season in Bangladesh

Climate change and climate variability are among the major threats to social, economic, and ecological conditions worldwide^1,2. According to the Intergovernmental Panel on Climate Change (IPCC, 2021), global warming is unequivocal, with widespread impacts on ecosystems, agriculture, and livelihoods. Climate variability arises from the constant changes between short-term weather deviations and long-run climate drifts³. Climate change has particularly devastating implications for global food security as it affects rainfall patterns, temperature, and relative humidity for an extended time.

The abovementioned factors are crucial in defining the specific region’s climate, and each of them has a major impact on the agricultural field productivity. Therefore, even the slightest changes in any of the factors can result in a significant decline in crop output^4,5,6. Extreme temperature plays a fundamental role in various domains, such as ecosystem stability, agricultural production, water availability, and sustainable socioeconomic development. Recently, there has been growing concern about temperature anomalies that deviate significantly from the average in a warming world. The last few years have witnessed increasing attention to the differences concerning temperature extremes between regions, as these extremes significantly affect both natural ecosystems and human communities^7,8,9. Rainfall is the most significant hydrological element, and it plays a vital role in relation to climate change. It influences food crops, harvesting, employment, and the environment. Flood and droughts, salinity, and temperature are in addition to rainfall projected to decrease rice and cereal output by eight and twenty-three per cent by 2050 from the base in 1990^10,11.

Bangladesh is one of the ten countries most vulnerable to the impacts of climatic extremes¹². The immediate impacts of high and low extremes in Bangladesh, like those occurring elsewhere globally, largely result from increased daily temperatures, changes in rainfall patterns, and other climatic events^13,14. Considering that Bangladesh is fundamentally an agricultural economy, with about 80% of the populace directly or indirectly dependent on agriculture, the consequences of climate variability and change on economic development and food security are particularly salient to the country¹⁵. Food security in Bangladesh will be threatened by lower agricultural yield, which is linked to climatic factors. Therefore, research is desperately needed to fully understand how much climate change is affecting Bangladesh and to support laws that would force the required steps to lessen its effects on agricultural output¹⁶.

Furthermore, it is commonly known that environmental factors, particularly rainfall, have a considerable influence on agricultural output¹⁷. Bangladesh has been severely impacted by the global and regional shifts in rainfall patterns brought about by climate change. The last 35 years have seen the nation go through unparalleled spells of rain and drought¹⁸. These variations have had unfavorable effects; for example, in 2007 a major flood that affected over 9 million people was caused by high rainfall¹⁷. In addition, natural disasters cost Bangladesh 111,289 metric tons in lost rice output in the 2019–2020 year¹⁵.

Moreover, it is predicted that from 2005 to 2050, Bangladesh’s overall rice output would decrease by 7.4% yearly due to temperature extremes¹⁹. Therefore, it is crucial to evaluate changes in the frequency, severity, and duration of severe temperature occurrences in addition to estimating changes in mean temperature trends²⁰. Still, there is little knowledge on the features of climatic extremes at the regional level. Therefore, it is imperative to monitor spatiotemporal fluctuations regularly and methodically in temperature extremes, especially in areas like Bangladesh^20,21

Like temperature and rainfall, relative humidity is another important element that affects crop productivity indirectly. High relative humidity reduces transpiration rates, which in turn decreases the upward movement of water and nutrients from the roots, ultimately lowering crop output. It also encourages disease and insect infestation. On the other hand, low humidity combined with sunny weather may lead to leaf wilting and in extreme cases, result in plant death due to water deprivation²². Therefore, it is crucial to assess whether there are trends in climatic parameters such as rainfall, temperature, and relative humidity, as well as patterns in their variability, to make accurate predictions.

Previous studies across Asia and beyond have highlighted significant changes in climatic variables and their implications for water resources, agriculture, and society. For example, a review of temperature trends in India stressed the importance of reliable detection methods²³, while spatial and temporal shifts in mean and extreme rainfall and temperature across urban centers in Rajasthan were reported, with direct consequences for urban water management²⁴. Similarly, climate change has already altered water availability and agricultural productivity in China, underlining the vulnerability of food security to climatic extremes²⁵. In Europe, the importance of data homogeneity for deriving valid climatic trends has been demonstrated^26,27, and change-point detection techniques continue to influence climate research²⁸. Non-stationarity in U.S. flood peaks has also been revealed, highlighting risks for long-term hydrological planning²⁹. Collectively, these studies emphasize that the choice of statistical technique, whether regression, non-parametric tests such as the Mann–Kendall (MK)^30,31 or Sen’s slope estimator³², or homogeneity tests^26,27,28, directly shapes how climatic trends are understood. The implications consistently point to increased climate risks for agriculture, water management, and disaster preparedness, thereby justifying the need for detailed region-specific analyses in countries like Bangladesh.

The Sen’s slope estimator is a widely applied non-parametric method for quantifying the magnitude of monotonic trends in climatic and hydrological variables³².In contrast to many statistical approaches that may be sensitive to serial correlation, time-series length, or seasonality, non-parametric techniques such as the Mann–Kendall test, Spearman’s rank correlation, and Sen’s slope estimator are widely applied for climate trend analysis due to their distribution-free nature and robustness. In this study, a novel trend analysis approach was also employed alongside the MK test to enhance the identification of trends^33,34,35. In Bangladesh, limited research has been conducted utilizing innovative trend analysis to detect trends in climate variables.

Even though trend analysis is useful for understanding past and present climate changes, future projections are more useful for policymakers in developing strategies that consider expected changes in climate variables. Physical models require substantial computational effort and resources, and while they are powerful tools for simulating climate processes, their outputs are not free from limitations. Systematic biases, such as sea surface temperature errors in regions like the Eastern North Pacific, highlight the need for careful evaluation and bias-correction before using climate model projections in impact studies. However, there are drawbacks to statistical methods like autoregressive (AR), moving average (MA), autoregressive moving average (ARMA), and autoregressive integrated moving average (ARIMA). For example, MA models use past errors as explanatory variables, while AR models only rely on past values. It is only appropriate to use ARMA models with stationary time series data. On the other hand, the application of artificial intelligence (AI) models, more specifically, machine learning (ML) techniques in climate forecasting has recently drawn interest^36,37. Unlike physical models, ML models perform very well by handling large and complex datasets with ease and without requiring a large amount of data³⁸. These models perform well even when dealing with non-stationary time series data³⁸. In a recent review, backpropagation-based models outperform neural network-based models for rainfall forecasting³⁶. The multilayer Perceptron (MLP) technique is used in our work to forecast meteorological variables like temperature, RH and rainfall. MLP has been widely praised for its reliable and consistent results when used to anticipate various hydrological and meteorological events^{39,40,41,42,43,44,45,46,47,76}. This study is unique for taking a novel approach to combining MLP models with other methods, such as Random Forest (RF), Neural Network, Artificial Neural Network (ANN), and Ensemble algorithm. We can estimate and predict climate variables like temperature, relative humidity, and rainfall with precision because of this thorough investigation.

Bangladesh Delta Plan 2100 (BDP-2100) identifies six climate-hotspot regions, of which the Barind, Coastal, and Haor areas are most critical for T. Aman rice production because of their vulnerability to temperature extremes, salinity intrusion, drought, and flash flooding⁴⁸. These regions contribute significantly to national rice output and represent contrasting agro-ecological conditions where temperature variability directly affects phenology, yield, and food security. Other hotspot regions, such as the Hill Tracts or riverine zones, were excluded as their cropping systems and microclimates differ substantially from T. Aman production areas, making them less comparable for this study’s focus.

Therefore, the current study aims to analyze the spatiotemporal trends of temperature variability during the T. Aman rice-growing season across Bangladesh and develop robust ML models to accurately predict temperature patterns, thereby supporting climate-resilient agricultural planning and decision-making. This study also aims to explore the spatiotemporal trends in monthly and seasonal minimum and maximum temperature in Bangladesh using the Sen’s slope estimator, and MK test. Furthermore, using MLP neural network, RF, Neural network, Artificial Neural Network (ANN) and Ensemble algorithm, the study attempts to forecast monthly and seasonal trends of the temperature in Bangladesh till 2030. The study also aims to evaluate the risk areas associated with changes in climatic patterns in Bangladesh by using Moran’s I. Additionally, techniques for analyzing change rates in data from a time series of climatic factors and variations will be employed, along with spatial mapping methods like Krigging and Inverse Distance Weighting (IDW).

Study area and data source

Bangladesh is in the northeastern region of South Asia, situated between latitudes 20°34′ and 26°38′ N, and longitudes 88°01′ and 92°41′ E. The study area comprised three distinct ecosystems—Barind, Coastal, and Haor—spanning 15 districts across Bangladesh. The Barind ecosystem includes the districts of Bogra, Pabna, Dinajpur, Rajshahi, Rangpur, and Sayedpur. The Coastal ecosystem encompasses Barisal, Bhola, Mongla, Khulna, Patuakhali, Sathkhira, and Khepupara. The Haor ecosystem covers Moulvibazar and Sylhet districts [Fig. 1]. Data was obtained from the BMD for these three ecosystems covering the period from 1961 to 2023. These datasets include Minimum and maximum temperature data from 35 meteorological stations across the three climatic hotspots (Coastal region, Barind region, and haor region) of Bangladesh, according to the Bangladesh Delta Plan (BDP)-2100. It should be noted that BMD has a total of 48 stations, including newly established ones, and stations with more than 15% of missing data were excluded to ensure accuracy.

Data Preprocessing

Missing data estimation

Missing data estimation is a critical procedure in meteorological research, particularly in the context of climate data analysis, where accurate and complete datasets are imperative for various applications. One commonly employed methodology for addressing missing data is the standard ratio approach, complemented by multiple imputations (MI) when confronted with data gaps⁴⁹. The standard ratio approach operates on the premise of leveraging information from neighboring stations to infer missing values. This technique presupposes a degree of correlation between rainfall, temperature, and RH measurements at proximate stations, thereby facilitating the estimation of missing data based on available observations from adjacent areas. The initial step in the standard ratio approach involves the computation of a weighting factor derived from the ratios between rainfall, temperature, and RH measurements at the target station and those at its neighboring counterparts. This weighting factor is expressed as:

Here, ({N}_{s}) denotes the average of the rainfall, temperature, and RH data presently available at the target station, ({N}_{i}) represents the average of the rainfall data at the ith neighboring station, and n signifies the total count of adjacent stations contributing to the computation. This methodology facilitates the estimation of missing data at the target station by applying the derived weighting factor to available observations from nearby stations. However, its efficacy may be constrained in scenarios where neighboring station data are inaccessible or when the correlation between stations is weak.

Outlier analysis

Outlier detection is a pivotal aspect of data analysis, particularly in meteorological studies where anomalous data points can significantly influence results.

Grubbs’ test and Box plot

Grubbs’ test is a widely utilized method for identifying outliers within a dataset. This statistical test is predicated on a two-sided hypothesis, wherein the null hypothesis posits the absence of outliers, while the alternative hypothesis suggests the presence of outliers either at the lowest or highest end of the data distribution⁵⁰.

The formula for Grubbs’ test is represented as follows:

where, x_i indicates the individual observation, (overline{x }) represents the average value of observations and s represents the sample standard deviation.

The critical value at significance level α for a sample of size n is given by:

where ({t}_{frac{alpha }{2n}, n-2}^{2}) is the critical value from the student’s t-distribution with n-2 degrees of freedom.

The decision rule is as follows: if G > G_critical, the most extreme observation is considered an outlier; otherwise, no outlier is detected.

In addition to Grubbs’ test, outliers were also identified using the Box Plot method, which relies on quartile statistics. The distribution is summarized by the lower quartile (Q₁), the median (Q₂), and the upper quartile (Q₃). The interquartile range (IQR) is defined as:

Potential outliers are determined using the following fences:

Observations lying below the lower fence or above the upper fence are classified as mild outliers. More extreme outliers can be identified using the thresholds,

Assessment of autocorrelation

Durbin-Watson statistics

Furthermore, autocorrelation within climatic variables data is assessed through Durbin-Watson autocorrelation statistics⁵¹.

The most used formula for Durbin Watson autocorrelation statistics is given by,

The Durbin-Watson statistics are calculated to ascertain the presence and direction of autocorrelation in the residuals ({text{e}}_{text{t}}) derived from the experimental observations. A Durbin-Watson statistic value of 2 indicates the absence of positive or negative autocorrelation, while values ranging from 0 to 4 are typically observed.

Normality test

Normality tests are indispensable in statistical analysis to assess whether a dataset follows a normal distribution, which is a fundamental assumption for many statistical techniques. Two commonly used methods for testing normality are the Q-Q (Quantile–Quantile) plot⁵². The Q-Q plot is a graphical technique that compares the quantiles of the dataset against those of a theoretical normal distribution. If the data points fall approximately along a diagonal line on the plot, it suggests that the dataset is normally distributed. However, deviations from this line indicate departures from normality. Specifically, if the points deviate systematically, such as bending upwards or downwards, it indicates skewness or kurtosis in the data, respectively.

Testing homogeneity

In climate research, ensuring the homogeneity of data is essential to mitigate the influence of measurement errors, instrument changes, and other environmental factors on the accuracy of trend analyses. In the present study, two standard tests were employed to assess homogeneity before trend analysis: Alexandersson’s standard regular homogeneity test⁵³ and Von Neumann’s ratio test⁵⁴. Alexandersson’s test is a widely recognized method for detecting structural breaks or shifts in data series, while Von Neumann’s ratio test evaluates the ratio of variances between consecutive segments of the dataset to identify significant deviations indicative of non-homogeneity. By implementing these tests, researchers can safeguard against erroneous trend conclusions, thus enhancing the reliability of findings in climate studies.

Standard normal homogeneity test (SNHT)

Standard Normal Homogeneity Test (SNHT)is a method utilized to detect abrupt shifts within time series data of meteorological and hydrological variables⁵³. The test, proposed by Alexandersson, involves computing a statistic (text{T}(text{k})) that compares the means of the first “k” years with the means of the subsequent “n-k” years in the dataset.

where, (k=text{1,2},3,dots dots dots , n).

If a structural break occurs at year “K,” the statistic (text{T}(text{k})) peaks around the vicinity of the break year, providing insight into the timing of the shift. The test statistic ({text{T}}_{0}) is defined based on this concept, where if ({text{T}}_{0}) surpasses the critical value, dependent on the sample size, the null hypothesis is rejected, indicating the presence of a structural break.

Von Neumann’s ratio test

The Von Neumann’s Ratio (VNR) test, proposed by Von Neumann⁵⁴, is a non-parametric statistical test used to assess the randomness of a time series by examining the presence of first-order autocorrelation in the residuals. For a series X_i with mean (overline{X }), the test statistic N is defined as: (text{N}={sum }_{text{i}=1}{left({X}_{i}-{X}_{i+1}right)}^{2}/{{sum }_{i=1}left({X}_{i}-overline{X }right)}^{2})

Method for trend analysis

M–K test

The MK test serves as a robust and widely employed method for evaluating trends in hydrological and meteorological time series data, particularly when the data does not conform to a normal distribution³¹. This nonparametric test effectively discerns monotonic movements, including both expanding and diminishing trends, within the dataset. It is particularly advantageous in scenarios where data may contain missing values or lack adherence to specific statistical distributions. The test statistics are given by,

where, (n) is the number of observations and the sign function is given by,

When the sample size (n) exceeds 10, the MK test statistic (S) is considered to follow an asymptotically normal distribution. In this scenario, the expected value E(S) of the test statistic is equal to 0. To calculate the variance of S, V(S), which provides a measure of dispersion around the mean, the following formula is employed:

The calculation of S involves considerations of ties (T) within the data, facilitating the assessment of variance V(S) and subsequently the determination of the standard deviation Z. By comparing the Z value to critical values at a chosen significance level, typically 95%, researchers can make informed decisions regarding the acceptance or rejection of the null hypothesis, thereby identifying statistically significant trends in the dataset.

Sen’s slope estimator

Sen’s slope model, widely utilized for estimating the magnitude of time series patterns, offers robustness compared to traditional regression slopes, as it remains unaffected by outliers or incomplete data⁵⁵. This method accurately estimates time series data by computing slopes for all pairs of ordinal time points and subsequently deriving the overall slope using the median of these individual slopes. Specifically, for N pairs of data, the slope estimate Q_k is computed using the formula:

Here, X_i and X_j represent data values at times i and j, respectively. The median of these N values of Q_k constitutes Sen’s slope estimator, which is calculated based on whether N is odd or even. A positive Q value indicates an upward (increasing) trend in the time series, while a negative value signifies a downward (decreasing) trend. When N is odd,

If N is even,

Spatial mapping method

In this study, the Inverse Distance Weighting (IDW) method was employed to interpolate the data and generate a spatial rainfall map^56,57.

The IDW interpolation follows the above Equation, where each scatter point is assigned a weight function ({w}_{i}), and ({f}_{i}) represents the specified function values at the scatter points. The total number of scatter points is denoted as ({n}_{i}). The conventional configuration of the weight function is computed using the following Equation,

In this context, ({h}_{i}) represents the distance between the scatter point and the interpolation point. The power parameter p is a positive real number selected freely. Here, ((x,y)) denotes the coordinates of the interpolation point, while (({x}_{i},{y}_{i})) denotes the coordinates of the scatter point.

Method for risk area assessment: Moran’s I

Local Indicators of Spatial Association (LISA) include various localized statistical methods that analyze the relationships between individual locations and their neighbors to identify spatial clustering patterns. In this study, we use the local Moran’s ({I}_{i}), a decomposition of the global Moran’s I, to quantify spatial autocorrelation^58,59. While the global Moran’s I provide an overall measure, the local Moran’s ({I}_{i}) assesses each observation’s relationship with its neighbors, classifying them based on the degree of spatial autocorrelation. The local Moran’s ({I}_{i}) is calculated using a weight matrix ({w}_{ij}) that defines the strength of the relationship between observation (i) and its neighboring observations (j). This standardized measure, like a z-score, identifies whether values are high or low relative to the mean, revealing patterns of clustering or dispersion. Additionally, the bivariate local Moran’s ({I}_{i}) extends this method by comparing different variables across time, thus facilitating the examination of temporal changes in spatial patterns.

For the ith region the local Moran’s I is defined as follows:

Note that the global Moran’s I is proportional to the sum of the local Moran’s I obtained for all region as follows:

Method for forecasting

Probabilistic and deterministic methods such as ARMA, ARIMA, and SARIMA have been widely applied for predicting hydrological and climatic datasets. These models assume linear relationships and are effective for capturing autoregressive and seasonal patterns; however, they may fail to adequately represent the non-linear interactions that are common in hydro-climatic systems. To overcome these challenges and enhance forecasting accuracy, this study explores the application of ML models, which are better equipped to capture non-linear relationships and complex patterns in hydro-climatic systems. Artificial intelligence (AI) models, particularly ANN, have gained popularity for addressing these limitations and solving environmental problems. As a non-linear black-box model, ANN mimics the structure and function of biological neural networks, consisting of interconnected neurons that process information to establish relationships between inputs and outputs. Using a feed-forward MLP algorithm, the model includes three layers: an input layer, one or more hidden layers, and an output layer. The number of neurons in the hidden layer was determined through trial and error. MLP is among the most widely used neural network architectures, particularly effective for classifying stationary patterns. It operates using the backpropagation algorithm, which minimizes errors and adjusts hidden layer parameters. Key features of MLP include non-linear processing elements (such as logistic or hyperbolic tangent functions) and extensive interconnectivity, where every element in one layer connects to all elements in the next layer.

Mathematical representation of the MLP model

Input Layer: The input vector in Fig. 2 is represented as ({varvec{X}}={left[{{varvec{x}}}_{1},{{varvec{x}}}_{2},boldsymbol{ }{{varvec{x}}}_{3},{{varvec{x}}}_{4}right]}^{{varvec{T}}}), where ({x}_{i}) represent the i-th input feature.

First Hidden Layer: The first hidden layer has five neurons. The transformation from the input layer to the first hidden layer can be described as:

where:

• ({{varvec{W}}}^{(1)}) is a 5 × 4 weight matrix for the first hidden layer.

• ({{varvec{b}}}^{(1)}) is a 5 × 1 bias vector for the first hidden layer.

• ({{varvec{Z}}}^{left(1right)}) is the linear combination of inputs before applying the activation function.

The output of the first hidden layer after applying the activation function f(⋅), typically the Rectified Linear Unit (ReLU), is given by:

Second Hidden Layer: The second hidden layer has three neurons, and the transformation from the first hidden layer to the second hidden layer is:

where:

• ({{varvec{W}}}^{(2)}) is a 3 × 5 weight matrix for the second hidden layer.

• ({{varvec{b}}}^{(2)}) is a 5 × 1 bias vector for the second hidden layer.

• ({{varvec{Z}}}^{left(2right)}) is the linear combination of activations from the first hidden layer.

Applying the activation function f(·) to ({{varvec{Z}}}^{left(2right)}):

Output Layer: The output layer consists of a single neuron to predict the continuous target variable. The output is computed as:

• ({{varvec{W}}}^{left(3right)}) is a 1 × 3 weight matrix for the output layer.

• ({{varvec{b}}}^{(3)}) is a scalar bias for the output layer.

• (hat{y}) represents the predicted output.

Several soft computing ML techniques are available for predicting and forecasting meteorological variables. Among these, the MLP based artificial neural network (ANN) has become widely popular among researchers. This popularity is due to the MLP’s computational simplicity and its ability to generate high-quality predictions. For instance, A study conducted in Poland introduced a stacked ML model that combines a MLP and a RF to forecast daily lake surface water temperature using daily air temperature as the input variable⁶⁰. In the field of meteorology, ML approaches have been widely utilized for predicting weather variables. A study on the development of an enhanced bidirectional recurrent neural network (Bi-RNN) comparing the performance of four ML approaches, including Bi-RNN, MLP, and RF⁶¹. Artificial neural networks (ANNs), particularly MLPs, have been employed in numerous forecasting models. For example, a study used an MLP to forecast monthly relative humidity in Delhi, India⁶². Overall, the combination of ML algorithms like MLP, RF, ensemble and neural networks, along with feature selection methods has demonstrated significant promise in predicting meteorological variables and other outcomes based on various input variables. In this study, we used MLP, ANN, RF, support vector machine (SVM) and the ensemble algorithm for predicting minimum and maximum temperature. Additionally, we combined these models to enhance prediction accuracy. The flow of proposed prediction methodology can be found in Fig. 3.

Convolutional neural networks (CNNs)

To identify local patterns and spatial–temporal relationships in temperature data, convolutional neural networks (CNNs) were used. CNNs employ convolutional layers in which hierarchical features are extracted by sliding filters (kernels) across the input. A non-linear activation (e.g., ReLU) and frequently a pooling step for dimensionality reduction come after each convolutional operation.

For a convolutional layer, the feature map at position (i,j) can be expressed as:

where x is the input, w_m,n are kernel weights, b is the bias term, and g(⋅) is the activation function. The extracted feature maps are passed through fully connected layers to generate the final temperature prediction.

Long short-term memory (LSTM)

Long Short-Term Memory (LSTM) networks were utilized for modeling sequential patterns in temperature data. Unlike traditional recurrent neural networks, LSTMs have memory cells and gating mechanisms (input, forget, and output gates) that govern how information is stored, updated, and transmitted between time steps. This makes LSTMs ideal for detecting long-term dependencies in time series data.

In this study, the LSTM was trained with lagged temperature readings (up to three years) and year, month, and day indicators as sequential inputs. Model hyperparameters, such as the number of hidden units, layers, learning rate, and batch size, were optimized for performance. The model was trained on normalized data and validated against an independent test set.

Model performance evaluation metric

To evaluate the performance of the model, four commonly used error metrics—MAE, MSE, RMSE, and MAPE were applied, as described below.

MSE: Mean Square Error (MSE) calculates the average of the squared differences between predicted and actual values⁶³. By squaring the errors, MSE gives more weight to larger discrepancies, which means it is more sensitive to significant errors or outliers. This characteristic makes MSE useful when large errors are particularly undesirable and need to be penalized more. However, since the errors are squared, the resulting value is in squared units of the original data, making it less intuitive to interpret compared to MAE or RMSE.

MAE: Mean Absolute Error (MAE) measures the average magnitude of errors between predicted and actual values, without considering their direction⁶³. It calculates the absolute differences between predictions and observations and then averages them over all data points. Since MAE treats all errors equally, it provides a straightforward and easily interpretable indication of the average prediction error in the same units as the data. However, it does not penalize larger errors more heavily, making it less sensitive to outliers compared to other metrics like MSE or RMSE.

MAPE: Mean Absolute Percentage Error (MAPE) expresses the average absolute errors as a percentage of the actual values, providing a scale-independent metric that facilitates comparison across different datasets⁶³. By showing errors relative to the size of the observations, MAPE offers an intuitive understanding of the prediction accuracy in percentage terms. However, it can produce extremely large or undefined values when actual data points are zero or close to zero, making it sensitive to small denominators and less reliable in such cases.

RMSE: Root Mean Square Error (RMSE) is the square root of the Mean Square Error, offering an error metric that is in the same units as the original data⁶³. Like MSE, it emphasizes larger errors due to the squaring of differences before averaging, but taking the square root brings the metric back to the original scale, making it easier to interpret. RMSE is widely used because of its sensitivity to large errors, making it an effective measure when significant prediction inaccuracies are particularly problematic.

Descriptive statistics of maximum and minimum temperatures in Bangladesh (1961–2023)

The descriptive analysis for maximum temperatures across ecosystems and regions reveals notable variations. Overall, the mean maximum temperature during the T.Aman season was 31.62 °C, with a standard deviation of 2.16 °C, ranging from 17.85 °C to 39.90 °C [Table 1]. Among ecosystems, Barind exhibited the highest mean maximum temperature (31.82 °C), followed by Coastal (31.67 °C) and Haor (31.62 °C) [Table 2]. The distribution was slightly left-skewed in all ecosystems, indicating a tendency for temperatures to cluster closer to the upper limit, as reflected by skewness values ranging between − 0.44 (overall) and − 0.35 (Barind).

For minimum temperatures, the overall mean was 24.02 °C, with a standard deviation of 3.10 °C, spanning a range from 9.62 °C to 29.98 °C [Table 1]. The Coastal ecosystem recorded the highest mean minimum temperature (24.31 °C), followed by Barind (24.06 °C) and Haor (24.02 °C) [Table 2]. The data demonstrated negative skewness across all ecosystems, with values ranging from − 1.42 (Barind) to − 1.53 (Coastal), indicating a concentration of temperatures near the upper end of the range.

Trend analysis of temperature

Maximum temperature

The MK test and Sen’s slope analysis revealed significant positive trends in maximum temperature during the T.Aman season across the Barind, Coastal, and Haor ecosystems. The Barind ecosystem exhibited a moderate upward trend, with a z-value of 23.61 and Kendall’s tau of 0.16, while the Coastal ecosystem showed a similar positive trend with a z-value of 23.88 and Kendall’s tau of 0.16. The Haor ecosystem demonstrated the strongest increasing trend, reflected by a z-value of 25.68 and Kendall’s tau of 0.17. Correspondingly, Sen’s slope estimates quantified the rate of increase as 23.61, 23.88, and 25.68 units per period for the Barind, Coastal, and Haor ecosystems, respectively [Tables 3 and 4].

The MK test and Sen’s slope analysis at the district level revealed significant positive trends in maximum temperatures during the T.Aman season across districts within the Barind, Coastal, and Haor ecosystems. In the Barind ecosystem, the highest trend was observed in Bogra, with a z-value of 24.99 and Kendall’s tau of 0.17, followed by Rajshahi (z-value: 24.63, tau: 0.17) and Pabna (z-value: 24.12, tau: 0.16). Sayedpur also exhibited a notable increasing trend with a z-value of 17.52 and tau of 0.15. However, Rangpur displayed a comparatively weaker trend with a z-value of 10.74 and tau of 0.07.

In the Coastal ecosystem, the district of Khepupara showed the strongest upward trend, with a z-value of 26.85 and Kendall’s tau of 0.19, followed by Pautuakhali (z-value: 24.64, tau: 0.17) and Barisal (z-value: 25.35, tau: 0.17). Bhola and Khulna exhibited moderate trends with z-values of 23.52 and 20.99, respectively. In contrast, Mongla (z-value: 13.24, tau: 0.11) and Sathkhira (z-value: 13.69, tau: 0.09) displayed relatively weaker upward trends.

For the Haor ecosystem, Sylhet recorded the strongest positive trend, with a z-value of 27.30 and Kendall’s tau of 0.18, while Moulvibazar demonstrated a moderate trend with a z-value of 17.78 and tau of 0.12. Sen’s slope analysis across all districts further confirmed the increasing rates of maximum temperature, with the steepest slopes observed in Khepupara (26.85) and Sylhet (27.30) [Supplementary tables S1 and S2].

Minimum temperature

The MK test and Sen’s slope analysis revealed positive trends in minimum temperatures during the T.Aman season across all three ecosystems: Barind, Coastal, and Haor. In the Barind ecosystem, the MK test yielded a z-value of 16.45 and Kendall’s tau of 0.11, indicating a moderate upward trend in minimum temperatures. Similarly, the Coastal ecosystem exhibited a stronger increasing trend with a z-value of 20.14 and tau of 0.13, reflecting a more pronounced rise in minimum temperatures. The Haor ecosystem also showed a positive trend, though slightly weaker than the Coastal ecosystem, with a z-value of 19.34 and tau of 0.13 [Tables 3 and 4].

The trend analysis for district-level minimum temperature during the T.Aman season, as indicated by the MK test and Sen’s slope test, reveals varying rates of change across districts in the Barind, Coastal, and Haor ecosystems. In the Barind ecosystem, the MK test z-values range from 5.62 (Sayedpur) to 14.11 (Bogra), with tau values reflecting moderate to weak upward trends, except for Sayedpur, which shows a weaker trend. In the Coastal ecosystem, the z-values range from 6.01 (Mongla) to 15.06 (Sathkhira), with tau values suggesting a moderate upward trend, particularly in districts like Sathkhira, Bhola, and Barisal, showing stronger increases in minimum temperature. The Haor ecosystem also shows variability, with the z-value for Sylhet reaching 20.17, indicating the strongest upward trend, while Moulovi Bazar presents a moderate trend with a z-value of 11.93. Overall, the trend analysis highlights a general but varied increase in minimum temperatures across these regions, with the Haor ecosystem exhibiting the most pronounced trend, especially in Sylhet [Supplementary tables S1 and S2].

Model performance evaluation for maximum temperature prediction

The evaluation of five ML models SVM, MLP, ANN, RF, and Ensemble across the Barind, Coastal, and Haor ecosystems identified the Multi-Layer Perceptron (MLP) model as the best-performing model in all regions. In the Barind ecosystem, the MLP model achieved the lowest MSE (0.82), MAE (0.69), MAPE (1.92), and RMSE (0.90), outperforming other models. Similarly, for the Haor ecosystem, MLP demonstrated superior accuracy with an MSE of 1.50, MAE of 0.91, MAPE of 2.74, and RMSE of 1.23. In the Coastal ecosystem, the MLP model consistently outperformed others, yielding an MSE of 1.47, MAE of 0.92, MAPE of 2.75, and RMSE of 1.21. While other models, such as SVM, ANN, CNN, LSTM, RF, and Ensemble, provided reasonable predictions, they exhibited higher error rates across all metrics, underscoring the MLP model’s superior capacity to accurately predict maximum temperatures in all three ecosystems during the T.Aman season [Table 5 and Fig. 4].

The evaluation of five ML models SVM, MLP, ANN, RF, and Ensemble across the districts identified the Multi-Layer Perceptron (MLP) model as the best-performing model in all regions. The MLP consistently outperformed other models, showing the lowest errors across ecosystems. For instance, in the Coastal ecosystem, the MLP achieved an MSE of 0.98, an MAE of 0.73, a MAPE of 2.28, and an RMSE of 0.99. Similarly, in the Haor ecosystem, the MLP reported an MSE of 1.17, an MAE of 0.82, an MAPE of 2.64, and an RMSE of 1.08. Other models like SVM and RF generally exhibited higher error values, while ensemble models provided a balance between performance and robustness. These findings emphasize the MLP’s efficacy in capturing local-level variations with higher accuracy and lower prediction errors [Supplementary Tables S3–S17; Figures S1–S15].

Model performance evaluation for minimum temperature prediction

Across all three ecosystems, the Multi-Layer Perceptron (MLP) model emerged as the best-performing model for predicting minimum temperatures. In the Barind ecosystem, MLP had the lowest MSE (0.48), MAE (0.46), MAPE (1.78), and RMSE (0.69), outperforming other models. The Coastal ecosystem showed similar results, with MLP achieving the best performance, with an MSE of 0.44, MAE of 0.48, MAPE of 1.85, and RMSE of 0.66. In the Haor ecosystem, MLP again excelled, yielding the lowest MSE (0.48), MAE (0.47), MAPE (1.85), and RMSE (0.69). Although models such as SVM, ANN, CNN, LSTM, RF, and Ensemble also provided accurate predictions, MLP consistently provided the best overall fit across all ecosystems, demonstrating its superior ability to predict minimum temperatures in the T.Aman season [Table 6 and Fig. 5].

The best model evaluation for district-level predictions across the different ecosystems reveals that the Multi-Layer Perceptron (MLP) consistently outperforms other models in terms of error metrics such as MSE, MAE, MAPE, and RMSE. In most districts, the MLP demonstrates the lowest MSE and RMSE values, indicating better predictive accuracy. For instance, in districts within the Barind and Coastal ecosystems, the MLP model achieves MSE values as low as 0.63 and 0.52, respectively, with corresponding RMSE values of 0.80 and 0.72, suggesting superior model fit and generalization. Additionally, the MLP maintains relatively low MAE and MAPE values across most districts, reinforcing its consistent performance. While models like SVM and Ensemble also show competitive performance, especially in specific districts, MLP generally proves to be the best overall model for district-level prediction of minimum temperatures [Supplementary Tables S3–S17; Figures S1–S15].

tenfold cross validation is performed for each model to check the validity [Table 7]. The multilayer perceptron (MLP) consistently produced the most accurate predictions of minimum temperature. The mean error in Barind was 0.72 with a standard deviation of 0.11, in Haor it was 0.88 with a standard deviation of 0.08, and in the Coastal zone it was 0.58 with a standard deviation of 0.18. For maximum temperature, MLP again performed best in Barind with a mean error of 0.88 and a standard deviation of 0.07, and in the Coastal zone with a mean error of 1.01 and a standard deviation of 0.12. In Haor, however, the support vector machine (SVM) achieved the lowest error with a mean of 1.11 and a standard deviation of 0.33. Other models such as ANN, CNN, LSTM, and random forest generally produced higher errors with larger variability across all zones. Overall, MLP emerged as the most reliable model, with SVM showing an advantage only for maximum temperature in Haor.

Assessment of risk area

Moran’s I analysis revealed distinct spatial patterns in temperature risk levels during the T.Aman season. For maximum temperature, all districts in the Barind region were categorized as severe risk (40 °C to < 42 °C), while the coastal region exhibited moderate risk across all districts (38 °C to < 40 °C). In the Haor region, Sylhet and Sunamganj fell under the moderate risk category (38 °C to < 40 °C), with the remaining districts classified as mild risk [Fig. 6].

For minimum temperature, Dinajpur, Rangpur, and Pabna in the Barind region were classified as very severe risk (temperatures below 4 °C), while other districts in the region showed severe risk. The coastal region uniformly exhibited moderate risk (6 °C to < 8 °C). In the Haor region, Sylhet was categorized as very severe risk, Moulvibazar as severe risk, and the remaining districts as moderate risk [Fig. 7]. All these categories were determined based on the official thresholds of the Bangladesh Meteorological Department (BMD, 2025)⁶⁴.

Agricultural implications

The spatial risk patterns and temperature predictions from the MLP model provide actionable insights for agricultural planning and climate adaptation in Bangladesh. Farmers can use these findings to adjust sowing dates to avoid periods of extreme heat or cold, select rice cultivars with enhanced tolerance to temperature stress, and implement microclimate management practices such as shading, mulching, or agroforestry to buffer crops against thermal extremes. Regional planning authorities can utilize risk maps to target adaptation resources more effectively, including optimized irrigation scheduling, development of shelter belts, and allocation of extension services for climate-smart farming practices. Collectively, these strategies can help mitigate the adverse impacts of temperature variability, stabilize rice yields, and enhance the resilience of farming communities to climate change.

Bangladesh has experienced significant shifts in temperature patterns from 1961 to 2023, with rising temperatures becoming a growing concern, especially during the T.Aman seasons. This analysis highlights regional disparities in climate change impacts on agriculture, water resources, and public health. The findings of the study found that during the T.Aman season, the overall mean maximum temperature is consistently 31.62 °C, while the mean minimum temperature is 24.02 °C, with significant variability observed across ecosystems and districts. The Barind and Coastal ecosystems stand out with relatively higher mean temperatures compared to other ecosystems. Districts such as Rajshahi and Rangpur tend to experience higher maximum temperatures, whereas Barisal and Bhola record lower minimum temperature extremes. Global Climate Model (GCM) projections suggest that winters globally will become warmer⁶⁵, a pattern supported by this study. Historically, January and April were the coldest and hottest months in Bangladesh, respectively. However, recent findings indicate May has become the hottest month in the southwestern region, and January temperatures are gradually declining, likely due to ongoing climate change. Furthermore, this study aligns with previous studies emphasizing the influence of temperature variability on rice productivity, particularly in coastal areas where temperature impacts are more pronounced than rainfall anomalies^66,67.

The Haor ecosystem demonstrated the strongest increasing trends in maximum temperature, while the Barind and Coastal ecosystems exhibited moderate upward trends. District-level analysis highlighted Sylhet in the Haor ecosystem and Khepupara in the Coastal ecosystem as experiencing the steepest increases. For minimum temperatures, the Coastal ecosystem exhibited the most pronounced upward trend, while the Haor ecosystem followed closely. Sylhet recorded the strongest district-level trend for minimum temperatures, emphasizing significant warming in the region. A gradual rise in temperature, particularly during the monsoon and winter seasons, has emerged as a concerning trend. Extreme heat can lead to crop wilting, sterility, and drought, while low temperatures can cause cold injuries and hinder plant growth and development^68,69. The combination of rising temperatures and limited rainfall poses significant risks to agriculture, including drought and crop loss, which can further drive humans^70,71.

The findings align with earlier predictions that temperature increases across Bangladesh are consistent with global warming trends⁷². Rice, a highly temperature-sensitive crop, is particularly vulnerable to these changes⁷³. Contributing factors include population growth, deforestation, and greenhouse gas emissions⁷⁴. Notably, the northern region, traditionally colder in winter due to its proximity to the Himalayan piedmont, may experience higher yields with increasing maximum temperatures, as warmer conditions could improve crop growth.

The MLP model consistently outperformed other ML models (SVM, ANN, RF, Ensemble) in predicting both maximum and minimum temperatures across ecosystems and districts. For maximum temperatures, MLP achieved the lowest errors across all ecosystems, with an MSE of 0.82, MAE of 0.69, MAPE of 1.92%, and RMSE of 0.90 in the Barind ecosystem. In the Coastal ecosystem, it yielded an MSE of 1.47, MAE of 0.92, MAPE of 2.75%, and RMSE of 1.21, while in the Haor ecosystem, the model showed superior accuracy with an MSE of 1.50, MAE of 0.91, MAPE of 2.74%, and RMSE of 1.23. Similarly, for minimum temperatures, MLP consistently outperformed the other models, achieving an MSE of 0.48, MAE of 0.46, MAPE of 1.78%, and RMSE of 0.69 in the Barind ecosystem, an MSE of 0.44, MAE of 0.48, MAPE of 1.85%, and RMSE of 0.66 in the Coastal ecosystem, and an MSE of 0.48, MAE of 0.47, MAPE of 1.85%, and RMSE of 0.69 in the Haor ecosystem. Several studies have highlighted the effectiveness of ML models, particularly Multi-Layer Perceptron (MLP), in predicting rainfall patterns. A study in Bangladesh demonstrated that the MLP model outperformed both statistical and other ML approaches in accurately predicting rainfall trends⁷⁵. Similarly, a study in India evaluated various statistical and ML techniques and found that MLP-based models provided the most accurate predictions for rainfall forecasting⁶². In line with these findings, our study also revealed that the MLP model outperformed other approaches in rainfall prediction. The superior performance of MLP can be attributed to its ability to capture complex, non-linear relationships in temperature data. Unlike traditional statistical models, MLP can efficiently handle large datasets with high variability, making it particularly suitable for temperature prediction in diverse climatic conditions.

The findings of the study revealed significant spatial heterogeneity in temperature risks during the T.Aman season, directly impacting rice production and agricultural sustainability. The Barind region’s categorization under severe and very severe risk highlights its susceptibility to extreme temperatures, which could adversely affect rice yields and exacerbate water stress. Mitigation measures such as introducing heat-tolerant rice varieties and efficient irrigation practices are essential. In contrast, the coastal region’s moderate risk levels may allow for relatively stable rice production, potentially benefiting from temperature moderation due to proximity to large water bodies. The diverse risk patterns in the Haor region, particularly Sylhet’s very severe minimum temperature risk, emphasize the need for localized adaptive strategies, including cold-resistant rice varieties and timing adjustments for planting. These insights are critical for developing region-specific climate adaptation strategies to safeguard rice production and ensure agricultural resilience.

The study focuses on T. Aman rice because of its significant contribution to food security and the alignment of its growth period with peak temperature variability, making the findings particularly relevant. We acknowledge that crops with different phenology, such as Boro, Aus, or wheat, or regions with differing microclimates and management practices, may experience distinct responses to temperature trends. Nonetheless, the methodological framework employed—integrating spatio-temporal trend detection, machine learning-based prediction, and spatial risk mapping, is broadly applicable. With appropriate data, this approach can be extended to other crops, seasons, or regions, providing a valuable tool for guiding national-level agricultural adaptation strategies under climate change.

Limitations

We acknowledge several limitations in our trend detection and spatial analysis. The Mann–Kendall (MK) test assumes independence in the time series; the presence of autocorrelation can inflate Type I error rates, leading to false detection of trends. Additionally, the power of trend detection is influenced by sample length, variance, and the magnitude of the slope, with small slopes or high variance reducing confidence in detected trends. Non-parametric methods like MK and Sen’s slope estimator are effective for identifying monotonic trends but cannot capture nonlinear or non-monotonic climatic fluctuations. Furthermore, spatial interpolation using Inverse Distance Weighting (IDW) has inherent limitations in capturing local heterogeneity, which may affect the accuracy of spatial temperature variability patterns. The performance of machine learning models in this study is inherently dependent on the quality and quantity of the input data. Since the models were trained primarily on historical temperature records, their predictions may not fully account for future climatic anomalies or unprecedented extreme events, which could affect the accuracy and reliability of future temperature forecasts. It is important to note that the spatial interpolation methods employed, such as Inverse Distance Weighting (IDW), may not fully capture fine-scale temperature variations. Consequently, the precision of spatial risk assessments and the identification of localized climate impacts could be affected. Future studies may benefit from integrating higher-resolution datasets or advanced geostatistical techniques to improve spatial accuracy and better inform localized adaptation strategies.

Implications and recommendations

To address the challenges posed by rising temperatures and climate variability during the T.Aman season, targeted policy interventions are essential to safeguard rice production and agricultural sustainability in Bangladesh. Research and development should prioritize the introduction of heat-tolerant and cold-resistant rice varieties suited to the specific climatic conditions of vulnerable regions, such as the Barind and Haor ecosystems. In addition, adjusting planting dates to avoid periods of peak heat can further reduce crop stress and yield losses. Region-specific adaptive strategies, including efficient irrigation systems, water management, and crop calendar adjustments, should be implemented to mitigate the adverse impacts of extreme temperatures and water scarcity. For the Barind region, investments in drought-resistant rice varieties and advanced irrigation practices are crucial, while in coastal areas, promoting saline-tolerant varieties and improved drainage infrastructure can stabilize yields. Expansion of public awareness campaigns and farmer training programs on climate-smart agricultural practices will strengthen capacity at the grassroots level. Furthermore, incorporating advanced machine learning models like MLP into agricultural planning can enhance the accuracy of temperature and climate risk predictions, enabling policymakers to design proactive, region-specific interventions. Collectively, these measures can help maintain rice productivity, support food security, and improve resilience against climate-induced stress during the T.Aman season.

This study reinforces significant trends in maximum and minimum temperatures across Bangladesh during the T.Aman season, with notable variations across ecosystems and districts. The Multi-Layer Perceptron (MLP) model consistently outperformed other machine learning models, demonstrating robust predictive capability. However, the results should be interpreted with caution. Prediction uncertainties exist due to reliance on historical data, and spatial interpolation methods introduce approximations that may affect localized risk assessments. Consequently, the findings are most reliable for the studied regions and seasons, and extrapolation to other crops, seasons, or regions should be undertaken carefully. Despite these limitations, the study highlights the urgent need for climate adaptation strategies and advanced predictive modeling to mitigate the adverse impacts of rising temperatures, providing guidance for policymakers and researchers to enhance agricultural resilience and support sustainable development in Bangladesh.