## 1 Introduction
Despite these advancements, challenges remain, particularly regarding data quality, feature selection, and model interpretability, highlighting the need for continued methodological refinement in RTA analysis. Although many studies focus on identifying accident risk factors, a growing number explicitly aim to forecast the future number of accidents. For instance, some works adopt ARIMA or exponential smoothing to predict monthly or daily crash frequencies [6; 19], whereas others leverage advanced data-driven methods--such as neural networks or hybrid models--to capture nonlinear trends in accident occurrence [23]. These forecasting approaches provide crucial insights for policy planning and resource allocation, yet they remain susceptible to the inherent volatility of traffic systems, underscoring the importance of ongoing research to enhance both accuracy and interpretability. Effective analysis of the impact of weather conditions on road accidents requires accurate information about the weather at the time of each recorded incident. However, historical data are typically collected only from a sparse network of meteorological stations, restricting direct observational capabilities. To overcome this issue, a range of spatial interpolation approaches has been explored in meteorology and hydrology, specifically designed to generate spatially continuous fields from limited point observations. Several studies have investigated simpler interpolation approaches such as Thiessen polygons, inverse distance weighting (IDW), and areal averages. [PERSON] et al. concluded that IDW offered a balance of simplicity and accuracy, but only when station density was relatively high [24]. Similarly, [PERSON] et al. underscored the significance of station density and local variability, arguing that precipitation is a rapidly changing, localized phenomenon [25]. [PERSON] et al. further highlighted the importance of high temporal resolution--often beyond the capacity of operational radars--particularly in urban areas [26]. Kriging has emerged as one of the most robust geostatistical methods for interpolating meteorological variables, thanks to its ability to exploit the spatial autocorrelation of observations [27]. In its simplest (ordinary) form, kriging estimates unknown values by modeling the variogram of the variable of interest [28]. Given the complexity of real-world environments, many studies have expanded this method by combining it with other data-driven or physical modeling approaches to refine spatial estimates of meteorological variables. Multivariate kriging methods, such as kriging with external drift (KED) or co-kriging, incorporate auxiliary variables like digital elevation models [29; 30], land cover [31], or solar irradiance [32] to correct for local climate effects and enhance predictive accuracy. Researchers also investigate temporal kriging approaches for filling gaps in time-series data [33]. The primary advantages of kriging lie in its ability to quantify estimation uncertainty and adapt to localized spatial structures, making it especially suitable for capturing fine-scale variations. However, it can be computationally demanding, and despite its robustness, its performance is also closely tied to the quality and density of the station network [24]. In light of these challenges, the present study aims to construct a comprehensive dataset that amalgamates detailed accident records with corresponding weather data and sociodemographic indicators, offering a promising avenue for investigation, facilitating the identification of statistically significant relationships among these variables and enabling the application of advanced statistical methodologies to rigorously quantify uncertainty and capture regional heterogeneity. By assembling and visualizing these interconnected data streams, this study seeks to illuminate potential links between environmental conditions and accident dynamics.