Exploring the Relationship between Melioidosis Morbidity Rate and Local Environmental Indicators Using Remotely Sensed Data

1. Introduction

The distribution of melioidosis, a bacterial infection caused by Burkholderia pseudomallei, can be examined in relation to the soil [1,2,3] and water resources [4]. Rainfall, soil water surfaces, and flooding are commonly associated with an increased incidence of melioidosis [5,6]. Melioidosis infections can be caused by the inhalation, skin abrasion, and ingestion of B. pseudomallei [4,7]. B. pseudomallei is found in soil depths of 0–90 cm [8], having an optimal growth temperature of 37°C [9]. It can survive in soil moisture of more than 10%, lasting for over a year at a 20% survival rate [10,11].

Changes in climate and environmental conditions lead to changes in health status, health-related illnesses, and death [12,13]. They also influence the transmission of infectious diseases such as melioidosis [14]. Epidemiology, the study of the distribution and determinants of health status or health-related events in a given population, puts the results of various studies to use for the protection and control of health problems. Monitoring and collecting environmental data over long periods can support the development and planning of disease prevention and surveillance programmes. These require the presentation of spatiotemporal disease occurrence. Remote sensing is used to detect and track the physical characteristics of an area, such as land surface temperature (LST), normalised difference vegetation index (NDVI), and normalised difference water index (NDWI), by measuring reflected and emitted radiation from a distance. To analyse the relationship of these indices with diseases, such as dengue, malaria, and leptospirosis, high-resolution satellite images are assigned pixel values for each index [15,16,17].

Currently, studies on melioidosis disease have only been limited to field remote sensing. Geospatial data, which are data with georeferenced coordinates to positions on the earth’s surface, have been utilised to monitor and investigate the risk area of melioidosis occurrence. Several studies have conducted spatial analyses of melioidosis distribution in the endemic regions of Australia [2,18], Thailand [19,20,21], and Laos [22]. Geostatistical modelling has revealed that the range distance of the spatial autocorrelation in a quantitative B. pseudomallei count was 7.6 m [23], and the range distance between positive B. pseudomallei samples was 90.51 m in a rice field [24]. To elucidate the conditions in which melioidosis infections can be increased in humans, an environmentally optimal bacterium, such as B. pseudomallei positive can be used.

The Google Earth Engine (GEE) is a cloud-based platform that revolutionises environmental analysis by granting access to vast collections of satellite imagery and geospatial data, empowering researchers to study diverse environmental factors such as land cover changes, climate patterns, and ecosystem dynamics. Leveraging its capabilities in time-series analysis, data visualisation, and algorithm development, the GEE facilitates the monitoring of global-scale environmental changes over time. It can be effectively used to study infectious diseases by integrating geospatial data and advanced analytics. For vector-borne diseases such as malaria, dengue fever [25,26], and COVID-19 [27,28], the platform can play a crucial role in understanding disease dynamics and improving public health interventions. In the case of melioidosis, a study revealed the spatial distribution pattern of melioidosis incidence using local Moran’s I and its spatial risk area using indicator interpolation kriging [21]. However, this study only considered a single variable. Furthermore, the limited information in a local area is not sufficient for understanding the characteristics and potential of melioidosis infection in other locations. Therefore, different environmental indicators in a local area should be utilised to reveal their relationships with melioidosis outbreaks.

To this end, this study aimed to investigate the relationship of melioidosis morbidity rate with local environmental indicators, specifically LST, NDVI, NDWI, and rainfall using remote sensing data and geographically weighted Poisson regression (GWPR). The objectives of this study were to: 1) determine the spatiotemporal dependence of melioidosis distribution and identify the monthly hot and cold spots and 2) classify the monthly data of melioidosis morbidity rate and the environmental indicators over a period of 10 y. The results of this study will be beneficial for the spatial monitoring and surveillance of melioidosis outbreaks in local areas.

2. Data and Methods

2.1. Study Area

This study was conducted in the Ubon Ratchathani province in Northeast Thailand. The province has an area of 15,774 km² and covers 25 districts, which are further divided into 219 sub-districts. It is adjacent to the borders of Cambodia and Laos and hosts three significant rivers: Mun, Chi, and Mekong. The study area experiences three seasons: summer, which begins in mid-February and lasts until mid-May; the rainy season, which extends from mid-May to mid-October; and winter, which is influenced by the northeastern monsoon winds, begins in mid-October and lasts until mid-February. The morbidity rate of melioidosis in the area was obtained from the National Disease Surveillance (Report 506), revealing a morbidity rate of 29.18 per 100,000 people in 2019 and 25.71 per 100,000 people in 2020 [29,30].

2.2. Conceptual Framework

This was an applied study of the spatiotemporal dynamics of provincial melioidosis outbreaks. It shows a map of sub-district boundaries using long-term data over a 10-y period, showing monthly outbreak periods. Spatial autocorrelation was used to determine the spatial distribution pattern, which showed clusters, randomness, and dispersion. The local spatial outlier and cluster revealed the location of the melioidosis hotspots in the study area. We required the understanding of local spatial indicators to track and determine the ones that influence melioidosis. This study used big data from satellite images at multiple time periods and scales through the GEE platform. Data were accessed by writing JavaScript commands and filtering the date, and then reducing the map using median values following the sub-district area. Spatial data were extracted from each sub-district and exported as a CSV file. The data were then linked to the spatial indicators (LST, NDVI, NDWI, and rainfall) and melioidosis according to the location of the sub-district. GWPR analysis was used to determine the relationship between spatial indicators and melioidosis morbidity rate in the local area for each month, as shown in Figure 1. The software used for analysing and visualising the spatial model were as follows: 1) spatial autocorrelation (Moran’s I) using Geoda v.1.20.0.22, 2) GWPR using AcrGIS Pro 3.2.0, and 3) map visualisation using R v.4.2.1 with the tmap package [31].

2.3. Melioidosis Data

The 10-y attribute data from 1 January 2013 to 31 December 2022, which included melioidosis case reports, were obtained from the public health office in Ubon Ratchathani. A data frame, which is a two-dimensional data structure aligned in a tabular manner in rows and columns, was created. The columns included the sex, age, occupation, Tambon ID (sub-district ID), and date of melioidosis cases confirmed by a doctor at the hospital. Then, the data were mutated, and the rows that were missing or not related to the study area were excluded. The morbidity rate per 10,000 people was calculated by dividing the number of melioidosis cases by the sub-district and number of people in the sub-district. Subsequently, the attribute data were connected one-to-one to the map of the sub-district boundary using the field Tambon ID. Geographical coordinates were based on the Universal Transverse Mercator (UTM) system zone 48N (EPSG 32648). Human protocols were approved by The Ubon Ratchathani University-Human Ethic Committee (ID# UBU-REC-171/2565).

2.4. Spatial Data

The GEE is a cloud-based system for processors, data storage, satellite remote sensing analysis, and other environmental and climate data products. The big data of high-resolution satellite images are freely accessible online via the JavaScript, Python, and R computer languages. By analysing the data collected from remote sensing technologies, the environmental indicators associated with the occurrence of a disease can be investigated within a specific locality.

In this study, the use of spatiotemporal analyses allows for a comprehensive understanding of the patterns and trends in the distribution of melioidosis over time. Moreover, the integration of spatial data enhanced our ability to explore the temporal dynamics of this disease and its relationship with various environmental indicators, namely LST, NDVI, NDWI, and rainfall. The application of remotely sensed data in this context provides valuable insights into the complex interactions between melioidosis and local environmental indicators that contribute to its spread. Overall, this approach enabled a detailed examination of the spatiotemporal aspects of melioidosis distribution, which can facilitate the identification of potential risk factors and the development of effective prevention and control strategies.

2.4.1. LST

The MOD11A1.006 product of the Moderate Resolution Imaging Spectroradiometer (MODIS) dataset provides daily global information on terrestrial LST and emissivity. Operated on the Terra satellite, this product has a high spatial resolution of 1 km. Emissivity information is particularly valuable for correcting temperature values based on the characteristics of different land surface types. Overall, the MOD11A1.006 product is a valuable resource for understanding the daily LST variations on a global scale.

2.4.2. Vegetation

The NDVI was acquired from the MOD13Q1 V.6 product, providing 16-d surface reflectance data for vegetation in the red (ρRED) and near-infrared (ρNIR) channels at a resolution of 250 m. NDVI is calculated using the formula: NDVI = (ρNIR − ρRED)/(ρNIR + ρRED). The contrast between the RED and NIR responses serves as a sensitive indicator of vegetation abundance, with maximum differences observed over areas with a full canopy and minimal contrast over regions lacking vegetation. In areas with low to medium vegetation density, the contrast results from changes in both the RED and NIR channels, whereas in areas with high vegetation density, the increase in contrast is primarily attributed to changes in the NIR channel as the RED band becomes saturated owing to chlorophyll absorption [32].

2.4.3. Soil Moisture

The NDWI is a satellite-derived metric from the ρNIR and Short-wave Infrared (ρSWIR) channels, designed for the estimation of water content within internal leaf structures. Employing the MOD09GA V.6 data, updated daily at a resolution of 463.313 m by MODIS, and subjected to cloud cover masking, the NDWI is computed using the formula: NDWI = (ρNIR − ρSWIR)/(ρNIR + ρSWIR) [33,34]. The SWIR channel captures alterations associated with vegetation water content and spongy mesophyll structure, whereas the NIR channel responds to leaf internal structure and dry matter content, excluding considerations for water content [34]. Through the amalgamation of NIR with SWIR, the index efficiently eliminated variations induced by leaf structure and dry matter, yielding a more precise evaluation of vegetation water content. The resultant NDWI product, expressed as a dimensionless value within the range of –1 to +1, not only signifies leaf water content, but also provides valuable insights into the vegetation type and cover.

2.4.4. Rainfall

The Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) data is a quasi-global (50S–50N), land-only rainfall dataset characterised by diverse spatiotemporal resolutions. Operating at a resolution of 5,566 m, the data contain amalgamated information from various sources, including ground station measurements and satellite data. An extension of a previously established climatology dataset, CHIRPS, leverages satellite information to fill data gaps in areas that lack ground station data. The dataset was constructed by integrating daily and monthly infrared cold cloud duration (CCD) precipitation estimates using an innovative blending procedure to incorporate the spatial correlation structure of the CCD. Acknowledging its utility on a global scale [35], the CHIRPS has been instrumental in assessing and monitoring precipitation patterns. Notably, CHIRPS has been employed to evaluate rainfall in various regions of Thailand [36,37,38].

2.5. Data Preparation and Pre-Processing

Remotely sensed data were meticulously selected and filtered based on the geographical area and date, employing the capabilities of a JavaScript API. Satellite data collection involved clipping and extension procedures aligned with district boundaries (Figure 2). Monthly median descriptive statistics were applied to reduce the pixel values within each tambon boundary. For rainfall data obtained from the CHIRPS product, the images were resampled to a 1,000 m resolution. Subsequently, data comprising the four environmental indicators were exported into a CSV file. This comprehensive dataset covered 219 administrative boundaries within Ubon Ratchathani. To establish spatial relationships, attribute data were connected using sub-district codes, linking melioidosis morbidity rate to environmental parameters. This method facilitated the creation of spatial data, enabling a nuanced exploration of the interplay between melioidosis and the chosen environmental indicators.

2.6. Spatial Statistics

2.6.1. Spatial Autocorrelation

Spatial autocorrelation, which examines the similarity among observation values at different locations, reveals distribution patterns, such as clustering, dispersion, and random distribution. In our investigation of melioidosis morbidity rate, we utilised the global Moran’s I statistic with queen contiguity-based weights, signifying a shared boundary edge between spatial units. Neighbouring relationships are assigned values of 1 or 0. Moran’s I ranges from −1 to 1, where positive values denote clustering, negative values indicate dispersion, and 0 suggests a random distribution. The statistical significance was set at a Pseudo P < 0.05.

The utilisation of Local Moran’s I in our study allowed for the identification of hotspots, cold spots, and spatial outliers by evaluating neighbouring relationships, providing insights into zones with either high or low morbidity rates and their spatial connections. Cluster and outlier detection, employing four autocorrelation types, was complemented by the Local Indicators of Spatial Association technique, which further scrutinised spatial patterns, emphasizing districts with similar morbidity rates surrounded by similar districts. Positive values assigned to features with high or low values among neighbours, along with the assessment of dissimilarity with neighbouring I values using negative values, contributed to a comprehensive understanding of the spatial dynamics. The incorporation of z-scores and p-values aided in evaluating the null hypothesis for significance and determining the output feature class for spatial dependency. This holistic approach provided a nuanced examination of the spatial distribution of the morbidity rate of melioidosis in various zones.

2.6.2. Global Poisson Regression (GPR)

A GPR model was used to examine the overall relationship between the morbidity rate of melioidosis and the environmental indicators in the study area. Before integrating the variables into the equation for GPR (Equation (1)), multicollinearity among the predictors was analysed to evaluate their independence. This analysis employed the variance inflation factor (VIF); values below 7.5 indicated non-multicollinearity and were consequently included in the model.

$$ln\left(y_i\right)=ln\left({\beta }_0\right)+\sum _{k=1}^p{\beta }_k{X}_{ki}$$

(1)

where ${ln(y}_i)$ is the expected value of melioidosis morbidity rate in the tambon I, ${\beta }_0$ is the global intercept, $X_{ki}$ is the k^th explanatory variable, and ${\beta }_k$ is the parameter estimating the explanatory variables.

2.6.3. Local Poisson Regression

The GWPR is a spatial statistical technique used in the analysis of spatially distributed count data, where the outcome variable represents the number of events or occurrences in a given area. This method is an extension of the traditional Poisson regression, which considers the spatial heterogeneity of the data. The GWPR method was employed to address spatial relationships and instability issues, specifically concerning count and ordinary data. Utilising an analytical framework derived from the GWR model [39], the GWPR extends its application to investigate the correlation between disease occurrence and spatial variability. This model, a form of conditional kernel regression, utilises spatial weighting functions to estimate locally varying coefficients within the Poisson regression parameters. This approach allowed the creation of local surface maps illustrating spatial fluctuations in the relationship between monthly melioidosis morbidity rate and localised distribution across the tambons. Consequently, this methodology facilitated the identification and in-depth exploration of locations that exhibited correlated relationships with the spatial indicators of LST, NDVI, NDWI and rainfall. The GWPR model was formulated using Equation (2).

$$ln\left(y_i\left(u_i, v_i\right)\right)=ln\left({\beta }_0\left(u_i, v_i\right)\right)+\sum _{k=1}^p{\beta }_k\left(u_i, v_i\right)X_{ki}$$

(2)

where Y represents the expected value of melioidosis morbidity rate at the coordinate location, $\left(u_i, v_i\right)$ denotes the two-dimensional coordinates of the i^th tambon, and ${\beta }_0$ and ${\beta }_k$ represent the locally estimated intercept and effect of variable k for location i, respectively.

The coefficients $\widehat{\beta }\left(u_i, v_i\right)$ were calculated using Equation (3) through spatial weighting. They were weighted by distance based on observations in nearby tambons, with data from tambons closer to the point being weighted by fixed and adaptive kernels. The optimal distance influenced the observed location, which was determined by the size of the nearest spatial unit.

$$\widehat{\beta }\left(u_i, v_i\right)={\left(x^{T}w\left(u_i, v_i\right)x\right)}^{-1}{x}^T w\left(u_i, v_i\right)Y$$

(3)

where $w$ is an n-by-n geographical weight matrix of the sub-districts.

$$w_{ij}=f\left(x\right)=\left\{\begin{array}{c}{\left[1-{\left(\frac{‖u_i-v_j‖}{G_i}\right)}^2\right]}^2, ‖u_i-v_j‖<G_i\\ 0 , otherwise\end{array}\right.$$

(4)

where $w_{ij}$ is the spatial geographical weight, $‖u_i-v_j‖$ represents the Euclidean distance between the tambon i and j, and $G_i$ represents the adaptive bandwidth size.

In the fixed kernel type, the data were weighted by a measure of the distance from the calibration location, considering data limitations. In the adaptive kernel type, the number of neighbouring tambons was optimised for each geographical region. Generally, Gaussian and bi-square kernel functions are used to produce a weight scheme when defining spatial units. We opted for a bi-square function with an adaptive kernel to estimate the spatial weights in the GWPR model based on the suggestions from previous research [39,40,41]. The optimal bandwidth was determined by lower AICc values, indicating the best model fit. Therefore, the AICc was also used to evaluate the model performance and compare the goodness-of-fit measurements between GPR and GWPR. Additionally, we considered the percentage of deviance explained as a measure of how well the model fit, indicating whether the explanatory variables can explain the relationship with melioidosis morbidity rate. Subsequently, the residuals of the GWPR regression model were assessed for spatial correlations among adjacent tambons using Moran’s I. If the spatial autocorrelation value was close to 0, the null hypothesis was accepted, indicating a perfectly random spatial pattern with a significance level of 0.05.

3. Results

3.1. Melioidosis Morbidity Rate

Between 2013 and 2022, a comprehensive analysis of melioidosis cases, which totalled 4,871, was conducted. Gender classification revealed 3,262 cases among males and 1,609 among females. Predominantly, those affected were farmers (1,943), followed by hired individuals (258). The monthly morbidity rates exhibited distinctive patterns, with January recording the highest rate, followed by August and September (Figure 3a). A closer examination of the monthly trends in the morbidity rate of melioidosis from January to December revealed that January had the highest morbidity rate, which gradually decreased until April. Conversely, morbidity rates increased from June to August, followed by a gradual decline until December. These findings provide valuable insights into the seasonal dynamics of melioidosis over a specified period.

As shown in Figure 3b, among the annual melioidosis rates between 2013 and 2022, 2022 recorded the highest rate, followed by 2017 and 2021, whereas 2019 marked the lowest rate. Notably, melioidosis rates exhibited a rapid increase from 2013 to 2017, averaging 4.96. Subsequently, a substantial decline occurred from 2017 to 2019, which represented the lowest rate in the decade. However, the morbidity rate experienced a notable resurgence, reaching 37.49 in 2022.

3.2. Spatial Autocorrelation

Table 1 illustrates the global Moran’s I values depicting the spatial autocorrelation pattern of the monthly melioidosis morbidity rate from 2013 to 2022, highlighting the spatial distribution clustering in the study area. January exhibited the highest spatial correlation. In contrast, the correlation gradually decreased from February to June, followed by a general upward trend, reaching 0.257 in December. The Local Moran's I results pinpointed the spatial clusters, indicating the hot and cold spots of melioidosis morbidity rate.

Spatial hotspots emerged predominantly in the northern part of the study area. The most prevalent spatial cluster (high–high) occurred in January and February, encompassing 29 tambons, followed by August with 23 tambons, and March and October with 20 tambons each. Notably, hotspots in May and June were scattered across the central and northern regions (Figure 4).

3.3. GPR Model

The GPR model associated with cover area represented the intercepts and coefficients of explanatory variables at a statistical significance of p ≤ 0.05 (Table S1). All explanatory variables were selected for the regression analysis because all values ignored multicollinearity when the VIF was <5. The results of the parameters in the GPR model are listed in Table S1. The statistically significant (p ≤ 0.05) months related to the morbidity rate of melioidosis were August, July, and January. In August, all explanatory variables affected melioidosis morbidity rate; the intercept and NDVI were negatively significant, whereas LST, NDWI, and rainfall were positively significant. Additionally, LST was not significant in January, whereas the other variables were significant at p ≤ 0.05. In July, rainfall was excluded from the model.

3.4. Local Poisson Regression

The summary of the descriptive statistics (minimum, median, and maximum) of the coefficients were based on the GWPR model results (Table S1). Table 2 shows the median coefficients representing the estimated impact of environmental indicators on the response variable, which is likely related to the occurrence of melioidosis. A negative intercept coefficient suggests a baseline reduction in the expected count. A positive LST coefficient indicates that higher temperatures were associated with the increase in the expected count of the response variable. A positive NDVI coefficient implies that areas with higher vegetation density corresponded to higher expected counts. A negative NDWI coefficient suggests a negative association between the water index and response variable. A positive rainfall coefficient indicates that increased precipitation was linked to higher expected counts. The values of these coefficients provide insight into the spatially varying relationships between environmental indicators and health-related outcomes, offering a nuanced understanding of the geographical complexities involved in the study. In addition, GWPR calibrated the association between melioidosis morbidity rate and LST, NDVI, NDWI, and rainfall for individual tambons.

3.5. Local Percent Deviance

Figure 5 shows the local percentage of the deviance-explained map that represents a potential relationship and fit accuracy of the explanatory variables that drive the changes in melioidosis morbidity rate. The spatial distribution of the percent deviance varies over time and space every month, which was inherently grouped by the Jenks natural breaks classification. The red colour on the map indicates a higher value of local percentage deviance, showing a clear pattern of spatial variation. The maximum local accuracy ranges were 0.600–0.658 in January, September, and December, respectively, with the highest obtained in December. In December, the risk areas were mostly located in the western region (29 tambons). In January, a trend of directional influence from the west to the central region was observed, while that in September was found in the southeastern region. However, for most of the northern region, the local accuracy ranges were 0.352–0.444 in July, and 0.239–0.338 in November.

3.6. Comparison between GPR and GWPR

The goodness-of-fit measures for both the GPR and GWPR models are presented in Table 3, revealing that the GWPR model exhibits a smaller AICc value than the GPR model. Conversely, the GWPR model exhibited a higher percentage of deviance than the GPR model. The global percentage of deviance ranged from 0.231 to 0.526, with the highest value observed in January. The residuals of the GWPR model were not statistically significant at a p ≤ 0.05, indicating a perfectly random spatial pattern. These findings suggest that the GWPR is a suitable method for elucidating the relationship between the morbidity rate of melioidosis and environmental indicators.

4. Discussion

Melioidosis remains an endemic disease in the study area, with cases reported throughout the 10-y study period. Hantrakun et al. [42] found in their study on the clinical epidemiology of 7,126 melioidosis patients in 60 hospitals in Thailand from 2012 to 2015 that the Ubon Ratchathani province exhibits a high incidence rate, designating it as a high-risk area. Figure 3 further confirms the occurrence of melioidosis outbreaks, demonstrating that the monthly and annual occurrences of the disease vary, thereby affecting the susceptibility of individuals to melioidosis each month. Given that B. pseudomallei, the causative agent of melioidosis, is found in both soil and water, individuals in at-risk populations are more susceptible to infection. Therefore, melioidosis is a significant public health concern that necessitates continuous monitoring, surveillance, and prevention efforts to mitigate its incidence in rural areas.

Because reliance solely on case report data for disease monitoring and surveillance may be insufficient, operations must incorporate spatial data to inform decisions and prevention plans. Spatial autocorrelation analysis using global Moran's I revealed a clustering pattern in the monthly distribution of the melioidosis prevalence. Local spatial correlation analysis using hotspot analysis facilitated the identification of areas with high and low patient numbers, distinguishing between low- and high-risk regions. For instance (Figure 4), a study identified 29 tambons (high–high) as hotspots in January and February, wherein closely clustered locations formed a high-risk group. Conversely, tambons initially deemed low-high may transition into high-risk areas because of their proximity to high-risk groups. Notably, the high-risk group was primarily clustered in the northern region, consistent with the findings of Wongbutdee et al. [21], who identified a significantly elevated incidence of melioidosis in northern Ubon Ratchathani during 2016–2020. Clustering of melioidosis cases suggests heightened exposure to B. pseudomallei in these areas [43]. Furthermore, B. pseudomallei has been isolated from the environment near patient residences in Northeast Thailand [44], indicating a correlation between the presence of melioidosis cases and B. pseudomallei in the surrounding environment. Thus, environmental factors likely contribute to the growth or persistence of B. pseudomallei in the soil and water, leading to the occurrence of melioidosis.

This study utilised satellite image data, specifically LST, NDVI, NDWI, and rainfall, to identify the environmental indicators influencing melioidosis occurrence. The analysis was conducted at a local scale, leveraging proximity-based data analysis, which enhanced predictive accuracy. This approach, employed through the GWPR model, outperformed the GPR model, as indicated by smaller AICc values and higher deviances (Table 3). However, the two models serve different purposes. While the GPR model elucidates global-level indicators influencing melioidosis development, the GWPR model highlights local relationships. Analysis of the explanatory variables in August and November, as well as in September and October, demonstrated a strong association between rainfall and melioidosis morbidity rate. Although previous studies have identified this relationship in numerous countries [5,45,46,47,48], it is not significant in Thailand [49]. The GPR model revealed significant associations between rainfall and melioidosis morbidity rate in certain months (e.g. January and August), but not in other months (March, April, May, and June) (Table S1).

The GPR model can inform policies for the prevention and control of melioidosis at a provincial level. However, its effectiveness is limited owing to variations in local environmental conditions such as temperature, humidity, rainfall, and vegetation. Therefore, satellite image data were also employed to facilitate the surface analysis of land cover, leveraging the ability of satellite images captured using the MODIS platform to monitor environmental changes over time. The GEE is a powerful tool for accessing large geospatial datasets, enabling the analysis and visualisation of geospatial image data in time series, which is instrumental in disease outbreak monitoring and our research, even though only a few studies have been conducted on melioidosis.

The GWPR model utilises the distance weighting of neighbouring locations to estimate the values at the points of interest. The adaptive kernel bandwidth yields a better-fitting model than the fixed kernel employed in this study. The kernel size is determined by the number of observations, with the distance adapted to the density of the nearest neighbours, resulting in a non-uniform spatial weighting shape. This assertion is supported by several previous studies [50,51,52]. The results demonstrate the percentage deviance, explaining the potential relationship between environmental indicators and the morbidity rate of melioidosis in each tambon. Consequently, the coefficients of the best-fitting model indicated the presence of non-stationarity, as evidenced by the different spatial patterns in the local coefficients of each independent variable. Notably, several local coefficients of LST, NDWI, and rainfall were negative, whereas the local coefficient of NDVI was positive (Figures S1–S12). Overall, the contribution of LST showed a negative correlation with melioidosis morbidity rate despite the association of B. pseudomallei contamination in soil and water with temperature. B. pseudomallei exhibited the highest growth rate at 37°C, with modest reductions observed at 30°C, 40°C, and 42°C, and a more pronounced delay at 25°C [53]. Additionally, the maximum temperature was associated with an elevated risk of melioidosis [46].

The NDVI product of the MODIS vegetation indices, produced at 16-d intervals, facilitates consistent spatiotemporal comparisons of vegetation canopy greenness, which is a composite property of leaf area, chlorophyll content, and canopy structure. During the summer period (mid-February to mid-May), few areas had vegetation cover, primarily dominated by dense paddy fields in the irrigation zones. However, previous studies have identified B. pseudomallei activity in soil paddy fields during the dry season [24] and in uncultivated lands [23]. In the rainy season, vegetation cover is increased, comprising paddy fields, grasslands, and forests, which B. pseudomallei has been detected in [54,55,56]. Notably, studies have shown no significant differences in B. pseudomallei activity between paddy fields and other land use types [56].

Heavy rainfall influences soil moisture and flooding, creating favourable conditions for the presence of B. pseudomallei in watershed areas [2,57,58]. NDWI and rainfall exhibited varying coefficients each month, aiding in understanding the spatial distribution of melioidosis morbidity rate at the local level. The use of satellite imagery enables rapid data acquisition and coverage of large areas, particularly in regions with differing rainfall patterns across tambons. Evaluation of rainfall from CHIRPS involves comparison with gauge observations before application, given its approximately 5 km resolution and large scale [35]. Consequently, the GWPR model assisted in weighting the parameter values for neighbouring observations to generate location-specific estimates. Our study revealed a negative coefficient for rainfall, with nearly half of the associated melioidosis cases displaying high deviance percentages. Nonetheless, an association between melioidosis and rainfall has been reported in various countries, such as Australia [45,46], Taiwan [47], Malaysia [48], and Singapore [5]. Additionally, Shaharudin et al. [59] reported a 1% detection rate of B. pseudomallei in soil, indicating a potential risk of melioidosis among flood victims in Kelantan, Malaysia.

In the present study, we utilised environmental indicators derived from remotely sensed data to investigate their spatial relationships with the morbidity rate of melioidosis, which differs in spatial heterogeneity in each local area. We employed boxplots to highlight measures of the central tendency of the explanatory variables (Figure 6), identifying outliers in each month and suggesting potential data non-stationarity. This observation indicated an inconsistency in the trend of monthly melioidosis morbidity rate over the decade, with cases occurring every month.

However, this study has several limitations. First, the resolution of the satellite imagery of environmental indicators may have affected the scale or resolution of the study area. Therefore, future considerations should incorporate higher-resolution or optimal-scale input data layers, such as Landsat, SPOT, and Sentinel. Second, our retrospective spatiotemporal design may have led to an underestimation of melioidosis cases within the province. Strengthening our study will involve predicting at-risk areas and forecasting melioidosis cases. Finally, while the GWPR model offers moderate reliability, it remains incomplete, similar to other spatial models. Addressing these challenges would involve validating the model performance through training data generation and testing data for validation and accuracy assessment for further research or implementation in public health. Additionally, we suggest further research that includes explanatory factors such as demographic and socioeconomic data, soil texture analysis, and health survey data. Furthermore, investigating the relationship between environmental indicators and bacterial presence in soil and water as well as analysing spatiotemporal data using time-series models will provide valuable insights for future studies.

5. Conclusions

This study conducted spatiotemporal analyses of melioidosis morbidity rate in an endemic area of Ubon Ratchathani over a 10-y period and categorised them into monthly occurrences. The distribution pattern of melioidosis morbidity rate indicated clustering, with hotspots predominantly observed in the northern region. The GPR model proved suitable for studying relationships at the global level and aiding in monitoring and prevention efforts at the provincial level. Meanwhile, the GWPR model was employed to estimate local geographical relationships by utilising nearby location weights to estimate the points of interest, yielding results superior to those of the GPR model. Our findings elucidate the relationship between environmental indicators and the morbidity rate of local melioidosis as indicated by the local percentage of deviance, which assesses the explanatory power of the model in terms of melioidosis occurrence. Although the coefficients for each factor varied from negative to positive, they effectively explained the relative contributions of the LST, NDVI, NDWI, and rainfall to the morbidity rate of melioidosis in each local area. Furthermore, the GWPR model revealed spatial heterogeneity attributable to differences in land surface characteristics and geography across various areas. Consequently, our findings offer valuable insights to guide the delineation of area boundaries when planning melioidosis monitoring and surveillance efforts.