## Materials and Methods
### Raster Data and Interpolation  
Annual average values of RH, WDR, WSP, and TEMP were interpolated using the Inverse Distance Weighting (IDW) method [36; 37]; the resulting maps were standardized using a pixel resolution of 500 m. A smooth surface was needed along with the preservation of the original average values. According to the Red de Meteorologia y Radiacion Solar (REDMET), the stations are strategically distributed across the basin; therefore, the distribution of points and the distance between them were expected to be even (geographical location of the REDMET stations are provided in Supplementary Materials Figure S1). Inverse distance weighting (IDW) assumes that each measured point has a local influence that diminishes with distance. It gives greater weights to points closest to the predicted location, and weights decline as a function of distance [36; 37]. This principle fits perfectly with this kind of data since what was sought was not the identification of a trend, but a continuous hypothetical surface based on the accumulated average values throughout a year. An interpolated surface was generated specifically for PM\({}_{2.5}\) to create the dependent variable map. We created a simple IDW interpolation just to show the raw values of PM\({}_{2.5}\) graphically. Nonetheless, an ordinary kriging interpolation method was used for analysis. Kriging has been tested several times to measure its reliability as a predictor and can accurately predict unknown values of a limited sample when contrasted to a set of known values, especially when the samples are randomly distributed [36; 38; 39; 40]. Kriging assumes stationary data, as in \"the condition that the probability of distribution of the variable is constant in time and/or space, meaning that the same function of distribution probability could be expected anywhere and anytime\" [41]. Due to the uniform distribution of the particulate matter, it was presumed to be stationary since the sources of pollution are located in fixed places, namely, roads, factories, houses, and public places. Therefore, the function probability of these variables withholds its properties [42]. For example, it is more likely for the surrounding area of a factory to be more contaminated than a place far away from it, which is the underlying approach of autocorrelation. Accordingly, an interpolated map was produced using the 37 PM\({}_{2.5}\) samples; in order to generate a better understanding of the spatial variability in PM\({}_{2.5}\), if a particular point had two different measurements (day or night), it was decided to use the one with the highest value of PM\({}_{2.5}\). It has been noticed that high values of ambient PM\({}_{2.5}\) concentrations do not necessarily reflect peaks of exposure [26; 43; 44]. We wanted to capture these high values, measured per personal monitor, to consider a hypothetical and possible acute toxic effect for further focused research. The nature of PM\({}_{2.5}\) implies the use of anthropogenic and topographic data to evaluate the relationship between the particles and the location of human activities [45]. It was assumed that the higher the RSL, the higher the values of PM\({}_{2.5}\) in the surrounding areas would be due to the way traffic behaves [45]. Therefore, road line data were used to approximately reflect this variation and produce a continuous surface for the RSL [4]. The lines of the roads were transformed into points separated every 250 m, which were then interpolated using IDW with a resolution of 500 m. Hence, buffered zones of the average speed value were used rather than fixed values. In the case of population density (PD), the 2005 census was used [46] in addition to the electoral sections (Geoelectoral Census Statistics) Estadisticas Censales a Escala Geoelectoral [46], since these were the most accurate sources available that fulfilled our purpose. The centroid of each electoral section was extracted and given the value of each polygon; then, an IDW interpolation was performed.