Identification of Risk Areas of Dengue Transmission in Culiacan, Mexico

Abstract

Dengue is a public health problem in more than 100 countries around the world and in virtually the entire region of the Americas, including Mexico. Mosquitoes of the genus Aedes aegypti transmit dengue; its reproduction requires certain geographical, epidemiological, demographic and socioeconomic conditions. Detailed information on socioeconomic, epidemiological and entomological data is available, but detailed meteorological information is not. The objective of this study was to identify the areas of risk of dengue transmission for each month of the year based on environmental, social, entomological and epidemiological information from 2010 to 2020, in Culiacan, Mexico. LST, NDVI and NDMI were calculated from Landsat 8 satellite images with remote sensing techniques. Additional variables were human population density and overcrowding; mosquito egg density from positive ovitraps; and probable cases of dengue. A descriptive analysis of the study variables and a multiple linear regression analysis were performed to determine the significant variables. In addition, a multicriteria spatial analysis was applied through the AHP technique to identify areas at risk of dengue transmission. The results revealed that the variables NDVI, NDMI and overcrowding were not significant; however, the LST, population density, egg density per positive ovitrap and probable cases were. The highest population in the transmission risk areas was in November, and the highest transmission area was identified in October. In conclusion, it was possible to identify which of the study variables were significant; in addition, monthly maps of risk areas of dengue transmission for Culiacan were obtained. Each geographical area had its own characteristics that influenced, in one way or another, the incidence of dengue, highlighting that the strategies for control of dengue must be specific to each region.

1. Introduction

Dengue is a disease transmitted mainly by the Aedes aegypti vector. Four closely related serotypes of the dengue virus are known: DEN-1, DEN-2, DEN-3 and DEN-4. Dengue causes between 50 and 100 million cases annually in more than 100 countries and is present in virtually the entire region of the Americas [1]. In the last three decades, the Pan American Health Organization (PAHO) has registered more than 24 million cases, of which 2.4 million have occurred in Mexico, demonstrating that dengue fever is, currently, a public health problem [2]. There are geographical, epidemiological, demographic and socioeconomic conditions that favor its transmission.

Most vector-borne diseases have been closely associated with tropical and subtropical climatic regions, while dengue, in particular, is found primarily in tropical and subtropical urban and suburban areas [3,4]. The mosquito that transmits dengue is considered domestic, since its behavior causes it to inhabit areas in close proximity to humans [5]. This disease results from the dangerous combination of environmental abandonment, climatic factors and growing poverty rates, all of which make it difficult to eradicate the mosquito. The growth of urban areas has promoted the establishment of large population groups in areas naturally conducive to vector development. The rapid concentration of people in urban areas is not always accompanied by a similar and efficient provision of public services, such as clean drinking water, drainage and garbage collection [6].

The General Directorate of Epidemiology (DGE) oversees epidemiological surveillance in Mexico, the Institute of Diagnosis and Epidemiological Reference (InDRE) is responsible for the diagnosis and reference of diseases, and the National Center of Preventive Programs and Disease Control (CENAPRECE) is responsible for vector control. These three institutions work in unison; the DGE and the InDRE provide data and information to the CENAPRECE, which leads the national vector control program and is responsible for vector control supervision.

Currently, entomological surveillance is carried out through two basic instruments, entomological surveys and ovitrap surveillance. The entomological survey is applied randomly, the interiors of the selected properties are inspected, the number of containers that are potential mosquito breeding sites is recorded, and the containers containing Aedes aegypti eggs are counted [7]. The ovitraps help determine the presence of Aedes aegypti eggs in samples from 15% of the blocks in the prioritized localities [8]. Entomological surveys are carried out sporadically in small samples, while surveillance with ovitraps is carried out weekly, with a uniform spatial or geographical representation, according to the guidelines of the corresponding operational guide.

The DGE registers the confirmed cases of dengue through the National Epidemiological Surveillance System (SINAVE); CENAPRECE has the database of confirmed cases to carry out the pertinent control interventions that are registered in the Comprehensive Vector Monitoring System (SIMV) [9]. The data from both the entomological and epidemiological surveillance are analyzed to generate information that helps identify the transmission risk areas and signals an alert using meteorological data and endemic channels, with the intention of prioritizing vector control actions, according to available resources [10]. It should be noted that the resources for vector control are very expensive; therefore, it is essential to focus on their optimization, especially in endemic countries with limited economies, such as Mexico.

Despite having epidemiological and entomological data with spatial and temporal references, dengue in Mexico continues to be a public health problem. A significant number of studies have shown that temperature, as one of several factors involved in the transmission of dengue, has a direct relationship with the abundance of the vector, and, therefore, with cases of dengue. Temperature also affects virus replication, maturation and the infection period [11]. The proliferation of the mosquito is also linked to humidity, which is linked to the vegetation providing a suitable environment for its development. These two variables are studied through the NDVI (normalized difference vegetation index) and the NDMI (normalized difference moisture index) [12]. The close relationship between these two indices support the evaluation of the functional state of the vegetation as a propitious condition for the development of mosquitoes [13,14]. A determining factor in the development of dengue outbreaks is the density of the human population in a geographical area; the greater the number of people, the easier it is for the mosquito to transmit the virus [15]. Factors such as urbanization, lack of basic services and variations in temperature, humidity and rainfall can increase the risk of epidemics globally as several insects, especially Aedes aegypti, will adapt to local environmental conditions [16]. Climate change can also be another important factor in the incidence of dengue, as evidence shows that the mosquito has now been detected in areas where the altitude is above 1700 m, where it was previously unknown [4,17]. Many studies have shown a strong relationship between dengue and the oviposition index [18,19]; however, other studies have found that the number of dengue cases correlates with the density of eggs, which depends on having favorable conditions for the proliferation of the mosquito [20,21].

The present study aimed to build a multicriteria spatial model that allows the identification of areas at risk of dengue transmission in each month of the year using environmental, social, entomological and epidemiological variables at the microregional level in the municipality of Culiacan, Mexico. To meet this objective, remote sensing techniques were used to generate environmental variables such as LST (land surface temperature), NDVI and NDMI from Landsat 8 satellite images. A multiple linear regression statistical analysis determined the significance of the variables, and a multicriteria spatial analysis identified the risk areas for dengue transmission.

This study may have implications for prevention and vector control actions by providing a detailed approach to identifying risk areas of dengue transmission for each month of the year for Culiacan. Culiacan is the capital of the state of Sinaloa; it is one of the 18 municipalities of the state, and is located on the coast of the Gulf of California and borders Durango in the Sierra Madre Occidental. In addition, Culiacan is the capital of the state of Sinaloa and is a municipality with a constant incidence of dengue.

2. Materials and Methods

Figure 1 shows the data management process for generating the spatial layers used to identify risk areas.

Monthly data for the LST, NDVI, NDMI, egg density per positive ovitrap, and the density of probable cases in the study area and time period were generated. Altitude data, population density data and housing overcrowding were obtained as spatial data only, considering that they are from a single measurement in time. The generation of the spatial layers from these data served two purposes. The first was to obtain the average of each variable with the basic geo statistical area (AGEB) analysis unit in order to perform a multiple regression analysis to obtain the significant variables. The second purpose was to use the weight of each variable that resulted from the statistical analysis in a multicriteria spatial analysis, through which the risk areas were obtained on a monthly basis.

2.1. Sources of Information

Information sources were as follows:

• Information from NASA and the US Geological Survey (USGS) Landsat 8 OLI-TIRS satellite imagery from 2013 to 2020;
• Digital elevation model of INEGI 2010;
• Statistical data on housing conditions from the INEGI Population and Housing Census 2010 and 2020;
• Entomological data from ovitraps of the SIMV used by CENAPRECE for the surveillance and control of vectors for the period 2010–2020;
• Cases of dengue from the SINAVE that regulates the DGE for the 2010–2020 period.

2.2. Study Variables

The study variables were as follows:

• Land surface temperature (LST);
• Normalized difference vegetation index (NDVI);
• Normalized difference moisture index (NDMI);
• Altitude;
• Population density;
• Percentage of inhabited private dwellings with some level of overcrowding;
• Density of eggs per positive ovitrap;
• Number and density of probable dengue cases.

2.3. Layer Generation Techniques and Information Analysis

The layer generation techniques were as follows:

• Satellite images of the study area were processed to obtain the LST, NDVI and NDMI through ArcGIS Desktop 10.6.1 software;
• Spatial data layers were generated using ArcGIS Desktop 10.6.1 software through thematic maps of altitude and the percentage of homes with some level of overcrowding;
• Spatial layers of population density, egg density and probable case density were generated through interpolation methods in ArcGIS Desktop 10.6.1 software;
• Information was processed through the SQL Server 2019 database manager to generate tabulated data of the entomological and epidemiological information by month for the entire study period;
• Spatial layers of the variables of interest were generated in raster format in ArcGIS Desktop 10.6.1 software;
• The statistical analysis was performed annually by AGEB in IBM SPSS Statistics version 20 and in RStudio version 2022.07.1;
• The spatial analysis was performed on a monthly basis in ArcGIS Desktop 10.6.1.

2.4. Study Area and Temporality

Figure 2 shows the study area.The municipality of Culiacan occupies 10.96% of the total area of Sinaloa and has 1483 localities with a population of 1,003,530, of which 48.93% are men and 51.07% are women. In the period from 2010 to 2020, Culiacan reported more than 10,000 probable cases of dengue, representing 37.2% of the total cases in the state [22]. There is a diversity of climates over the municipality: very warm and warm–dry (37.40%), very warm and warm semi–dry (31.96%), warm sub–humid with medium humidity (27.98%), warm sub–humid with low humidity (1.49%), semi–warm sub–humid with medium humidity (1.13%), and semi–warm sub–humid with low humidity (0.04%). Temperatures range from 18 to 26 °C and the annual precipitation ranges from 400 to 1200 mm. Culiacan is located between parallels 24°02′ and 25°17′ north latitude; the meridians 106°52′ and 107°50′ west longitude; and varies in altitude between 0 and 1800 m [23].

2.5. Spatial and Temporal Variables

2.5.1. Remote Perception

A total of 132 Landsat 8 satellite images from the period 2013–2020 were generated, and the LST, NDVI and NDMI were obtaining. The monthly average of all the images of the study period was calculated using the following formula:

$$\overline{x}Va{r}_i=\frac{{{\displaystyle \sum }}_{j=1}^{n}Va{r}_{ij}}{n}$$

(1)

where Var is the processed image of the LST, NDVI or NDMI, i is the month (1–12), and n is the number of images obtained in the month of the study period.

The Landsat 8 satellite has two sensors: the operational terrestrial imager (OLI) and the infrared thermal sensor (TIRS) [24].

Table 1. Spectral ranges of the OLI sensor bands.

Band	Description	Wavelength	Resolution
1	Coastal spray	0.43–0.45 µm	30 m
2	Blue	0.45–0.51 µm	30 m
3	Green	0.53–0.59 µm	30 m
4	Red	0.64–0.67 µm	30 m
5	Near Infrared (NIR)	0.85–0.88 µm	30 m
6	Short Wave Infrared 1 (SWIR1)	1.57–1.65 µm	30 m
7	Short Wave Infrared 1 (SWIR1)	2.11–2.29 µm	30 m
8	Panchromatic (PAN)	0.50–0.68 µm	15 m
9	Cirrus	1.36–1.38 µm	30 m

details the information for the OLI sensor bands.

The TIRS bands are acquired at a resolution of 100 m, but are resampled at 30 m.

Table 2. Spectral ranges of the TIRS sensor bands.

Band	Description	Wavelength	Resolution
10	Thermal Infrared Sensor 1 (TIRS 1)	10.6–11.19 µm	100 m
11	Thermal Infrared Sensor 2 (TIRS 2)	11.5–12.51 µm	100 m

shows the information for the TIRS bands [25].

Land Surface Temperature (LST)

The split window (SW) method was used to calculate the LST with the following formula:

LST = TB₁₀ + C₁ (TB₁₀ − TB₁₁) + C₂ (TB₁₀ − TB₁₁)² + C₀ (C₃ + C₄W) (1 − ε) + (C₅ + C₆ W) ∆ε

(2)

where LST is the land surface temperature in degrees Kelvin (K), C₀ to C₆ are the values of the coefficient for SW [26], TB₁₀ and TB₁₁ are the brightness temperatures of bands 10 and 11 (K), ε is the mean value of the land surface emissivity (LSE) of the TIRS sensor bands, W is the water vapor content in the atmosphere, and ∆ε is the difference in LSE [26].

This method requires the conversion to reflectance in the ceiling of the atmosphere (TOA) and the brightness temperatures of bands 10 and 11. The conversion to reflectance in the TOA is calculated with the following formula:

Pλ = ML Qcal + AL

(3)

where Pλ is the spectral radiance value of the TOA (watts/(m² × srad × μm)), M_L is the band minus the specific scaling multiplicative factor obtained from the metadata, A_L is the band minus the specific escalation additive factor obtained from the metadata, and Q_cal is the standard product quantified and calibrated by pixel values (DN). This value refers to each of the bands in the image.

To perform the conversion at brightness temperature on the satellite, the following formula was used:

$$T=\frac{K2}{\ln \left(\frac{K1}{L\lambda }+1\right)}$$

(4)

where T is the apparent brightness temperature in Kelvin (K); Lλ is the reflectance on the roof of the TOA atmosphere (watts/(m² × srad × μm)); K1 is the conversion constant $K1$ specific to each band, which is the thermal constant is supplied in the metadata; and K2 is the specific K2 conversion constant for each band, which is the thermal constant is supplied in the metadata [25].

The satellite images of the red and infrared bands were corrected for the effect of electromagnetic energy in the water particles suspended in the atmosphere, using the dark subtract method. The correction is applied assuming a reflectance percentage of 1% in dark areas [27].

Normalized Difference Vegetation Index (NDVI)

Once the images have been corrected, the calculation of the normalized difference vegetation index must be made from the following formula:

$$NDVI=\left(\frac{\left(NIR-R\right)}{\left(NIR+R\right)}\right)$$

(5)

where NIR is the infrared band pixel value (NIR) and R is the pixel value of the red band (R) [28].

The fractional vegetation cover (FVC) is an index that estimates the fraction of the area occupied by vegetation, and is obtained from the NDVI. This index is required to find the land surface emissivity (LSE) values:

$$Pv=\left[\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}\right]$$

(6)

where NDVI_max is dense vegetation and NDVI_min is the bare ground.

The LSE is calculated with the following formula:

LSE = ε_s × (1 − FVC) + ε_V × FVC

(7)

The values for the emissivity of the soil and vegetation of each band are found in the

Table 3. Emissivity for TIRS bands 10 and 11.

Emissivity	Band 10	Band 11
Land	0.971	0.977
Vegetation	0.987	0.989

.

The concluding calculation corresponds to the application of the SW algorithm and the calculation of the average value of the LSE (ε):

$$\Delta =\frac{LS{E}_{10}+LS{E}_{11}}{2}$$

(8)

and the difference in LSE values (∆ε):

∆ ε = LSE₁₀ − LSE₁₁

(9)

Finally, the temperature is converted from degrees Kelvin to degrees Celsius with the following formula:

0 Kelvin = −273.15 Celsius

(10)

Normalized Difference Moisture Index (NDMI)

The normalized difference moisture index was calculated with the following formula:

$$NDMI=\left(\frac{\left(NIR-SWIR 1\right)}{\left(NIR+SWIR 1\right)}\right)$$

(11)

where NIR are the pixel values of the near infrared band (NIR), and SWIR 1 are the pixel values of the short wave infrared band 1 (SWIR) [29].

For the generation of the spatial layers of the calculation of the LST, NDVI and NDMI, a model was generated in ArcGIS model builder.

Egg Density Per Positive Ovitrap

A positive ovitrap density layer was devised by the Aedes aegypti egg density interpolated for each month from the average number of eggs from a total of 2421 ovitraps during the study period. This used the interpolation spatial analysis method of inverse distance weighting (IDW). This method determines cell values through a linearly weighted combination of a sample point set and assumes that the variable being mapped decreases its influence at a greater distance from its sample location [30]:

$$z_p=\frac{{{\displaystyle \sum }}_{i=n}^n\left(\frac{z_i}{d_i^P}\right)}{{{\displaystyle \sum }}_{i=1}^n\left(\frac{1}{d_i^P}\right)}$$

(12)

where z_i corresponds to the values that are known, n corresponds to the number of values that entered the set search radius, $d_i^P$ corresponds to the distance between each point and the site of interest, and P corresponds to power.

Dengue Case Density

An interpolated probable case density layer of geo-referenced probable cases was made for each month of the study period for a total of 10,450 reported cases. This used the spatial analysis method called IDW.

2.6. Spatial and Temporal Variables

2.6.1. Altitude

The altitude was generated from the digital elevation model (MDE); this model is a visual and mathematical representation of the height values with respect to mean sea level [31]. An altitude thematic map was made and subsequently converted to the raster format.

2.6.2. Population Density

The population density was obtained through the spatial analysis method called kernel; as a result, the population density per square kilometer is obtained. A magnitude per unit area is calculated from a spatial layer of points generating a smooth surface for each point [32]:

$$\begin{gathered} Density=\frac{1}{{\left(radius\right)}^2}{\displaystyle \sum }_{i=1}^n\left[\frac{3}{\pi }·po{p}_i{\left(1-{\left(\frac{dis{t}_i}{radius}\right)}^2\right)}^2\right] \\ \mathrm{For}\mathit{dist}\mathrm{i} < \mathit{radius} \end{gathered}$$

(13)

where i = 1, …, n are the entry points. The points in the sum are included if they are within the radius distance of the location (x, y); pop_i is the population field value of the point i, which is an optional parameter; and dist_i is the distance between the point i and the location (x, y).

2.6.3. Housing Level of Overcrowding

A spatial layer of the percentage of overcrowding at the AGEB level was generated; this is the maximum level of aggregation that is available. Overcrowding is defined by the percentage of homes that have 2.5 or more people per bedroom [33].

2.7. Statistical Analysis

Multiple linear regression was used to analyze a mathematical association, through a linear equation, between the value of a dependent variable and the value of one or more independent variables. Multiple regression analysis examined the influence on the response or dependent variable, corresponding to the number of probable cases to detect areas of dengue transmission, through the following spatial and temporal variables:

• Land surface temperature;
• Normalized difference moisture index;
• Normalized difference vegetation index;
• Density of eggs per ovitrap;
• Density of probable cases;
• Altitude;
• Average number of inhabited homes with some level of overcrowding;
• Population density per km².

The statistical analysis used the information from all the variables as reported at the AGEB level. Multiple linear regression analysis refers to the extension of a simple linear regression model with more than one independent variable—when determining the behavior of the variable Y in a regression hyperplane from an optimal combination of predictor, or independent variables referring to X₁, …, X_n.

There are different aspects to consider in a multiple linear regression model, such as (a) the regression equation, (b) the validity and fit of the model, and (c) the analysis of assumptions.

The general regression equation is represented by:

Y = B₀ + B₁ X₁ + B₂ X₂ + … + B_n X_n + E_i

(14)

where B₀, the intercept, and B_i represent the partial regression coefficients for each independent variable and are interpreted as the change that is generated in the dependent variable, Y, for each unit that a predictor varies, B_i, with the rest of the constant predictors. The estimation of the regression parameters is carried out such that the sum of the squares of the errors or residuals is minimal and optimizes the prediction [34].

2.8. Spatial Analysis

A multi-criteria spatial analysis was carried out on a monthly basis through the hierarchical analytical process (AHP). This method comprises a set of techniques that allow for the evaluation of various choice alternatives in the light of multiple criteria (study variables) and priorities (weights). From the spatial point of view, the alternatives are observation units or portions of the territory that are evaluated based on their geographical characteristics. This method uses the weighted linear combination [35].

The AHP uses paired comparison matrices using a fundamental scale.

Table 4. Fundamental pairwise comparison scale.

Value	Definition	Comments
1	Equal importance	Criterion A is just as important as Criterion B
3	Moderate importance	Experience and judgment slightly favor criterion A over B
5	Big importance	Experience and judgment strongly favor criterion A over B
7	Very great importance	Criterion A is much more important than Criterion B
9	Extreme importance	The greater importance of criterion A over B is beyond doubt.
2, 4, 6 and 8	Intermediate values between the previous ones, when it is necessary to clarify

shows the possible values of the scale ranging from 1 to 9.

This matrix complies with the properties of reciprocity, homogeneity and consistency (the matrix must not contain contradictions in the assessment performed) [36].

Consistency is obtained through the consistency index (CI):

$$CI=\frac{{\lambda }_{max}-n}{n-1}$$

(15)

where λ_max is the maximum eigenvalue and n is the dimension of the decision matrix.

A multi-criteria evaluation comprises a set of techniques that allow for the evaluation of various choice alternatives in light of multiple criteria and priorities.

The random consistency index (RCI) is calculated with the following formula:

$$RCI=\frac{1.98\times \left(n-2\right)}{n}$$

(16)

where n is the dimension of the decision matrix.

Once the consistency index and the random consistency index have been obtained, the consistency ratio (CR) is obtained with the following formula:

$$CR=\frac{CI}{RCI}$$

(17)

The CR is accepted as long as it does not exceed the values of the

Table 5. Maximum percentages of the CR consistency ratio.

Matrix Size (n)	Consistency Rate
3	5%
4	9%
≥5	10%

.

Once the consistency has been verified, the weights are obtained, which represent the relative importance of each criterion or the priorities of the different alternatives with respect to a certain criterion. This process requires the normalization of the criteria for using the eigenvalue method, where the following equation must be solved:

$$A ·w={\lambda }_{max} ·w$$

(18)

where A represents the comparison matrix, w, represents the eigenvector or preference vector, and λ_max the eigenvalue [37].

To conclude with the AHP, the following equation is applied:

Y = C₁ × X₁ + C₂ × X₂ + C₃ × … + C_n × X_n

(19)

where Y is the dependent variable, C_n are the coefficients or weights, and X_n are the independent variables.

3. Results

3.1. Spatial and Temporal Variables

A total of 12 raster images were obtained for each variable (LST, NDVI and NDMI), resulting in a total of 36 images. On average, 12 images of each variable were processed per month of the 8 years that comprise the study period.

A total of 12 raster images were obtained for the variable density of eggs and density of cases, resulting in a total of 24 images.

The Supplementary Materials show the complete series of monthly maps for these variables and demonstrate that the incidence of cases was lowest in April and highest in October with 1.65 and 28.78%, respectively.

3.1.1. Land Surface Temperature

Figure 3 shows that the range of the average LST oscillates between 28 and 49 °C during the year; January has the lowest average temperature, while June has the highest average temperature. There is a trend of continuous increase during the first half-year from 28 to 49 °C, while in the second half there is a gradual decrease, with a marked seasonality between the months of August, September and October at around 40 °C, subsequently, the temperature declines from November to January.

In the majority (96.69%) of the AGEB, the temperature ranges from 40 to 45 °C in April (Figure 4a), and there were 33.62% reported cases during this time. In October (Figure 4b), 59.50% of the AGEB have temperatures ranging from 35 to 40 °C, and there were 76.85% reported cases of dengue during this time.

3.1.2. Normalized Difference Vegetation Index

Figure 5 shows that the vegetation index decreases slightly in the first half of the year. In June there is a slight increase in vegetation, and this continues until September, after which it declines slightly from October to December.

The green areas of the urban zone were represented by an NDVI, ranging from 0.35 to 0.75, which concurs with other studies [38]. Of the entire study area, in April, the green areas represented 4.14,% and in October, they reached 26.75%.

Figure 6a shows the NDVI map for April, where 99.72% of the AGEB have a built-up area or bare soil, in which 67.13% report no dengue cases. Figure 6b shows that in October, where 15.43% of the AGEB have urban tree vegetation, 55.36% reported cases of dengue.

3.1.3. Normalized Difference Moisture Index

Figure 7 shows the annual behavior of the NDMI. In the first half of the year there is a slight decrease in humidity, and the data are very close to the average, indicating that there is little variability in the humidity index in these months. From July, an increase in the humidity index is noted, and it reaches its maximum point in the rainy season in August and September. In October, the humidity index decreases until December and returns to near average values.

Of the entire study area, 9.48% has a high humidity in April, and this increases to 22.76% in October.

Figure 8a shows the NDMI map for April when 94.77% of the AGEB have zero humidity and 66.57% report no dengue cases. Figure 8b shows October when 7.99% of the AGEB have a built-up area or bare soil and 48.28% report cases of dengue.

3.1.4. Egg Density Per Positive Ovitrap

Figure 9 shows the seasonality of the density of Aedes aegypti eggs and an increase with the rainy season. February has an average density of 26 eggs, whereas September shows the highest average with 82 eggs. From December to February, there is a low variability between monthly averages and the density of eggs as 50% of the data ranges from 25 to 36 eggs. This increases in August and continues to rise in September, with an average of 83 eggs. From October, a decline occurs, resulting in an average of 30 eggs in December.

Of the entire study area, 95.59% of the AGEB have positive ovitraps for April and October. On average there are 17 ovitraps per month with 435 records of each ovitrap. The months with the highest egg densities are August, September and October with 14.95, 16.63 and 14.97%, respectively.

A layer of positive ovitrap egg density was obtained for the analyses carried out; derived from this layer in Figure 10a, you can see the map of eggs by positive ovitraps for April, where 90.08% of the AGEB present an egg density less than 50. Figure 10b shows that in October, 96.42% of the AGEB are above 50 eggs per positive ovitrap.

3.1.5. Dengue Cases

In the municipality of Culiacan in the period 2010–2020, there were 10,450 probable cases of dengue reported. Figure 11 shows that in the first quarter, a decrease in cases is noted, and the month with the lowest incidence of cases is April at 1.7%. The increase in cases begins in August and rises until October with 30.4% of total cases. In November, cases decrease until December.

In April, 67.22% of the AGEB did not report cases of dengue, whereas in October, there were 81.82% of the AGEB that reported cases of dengue.

A case density layer was obtained for the analyses, and in Figure 12a, derived from this layer, the case density map for April shows that in 89.53% of the AGEB, there were no reported dengue cases. Figure 12b shows that in October, 89.26% of the AGEB reported cases of dengue.

3.2. Spatial Variables

There is one raster image for each variable (altitude, population density and percentage of homes with some level of overcrowding), which is three images in total. These variables have neither a spatial, nor a temporal, variability.

3.2.1. Altitude

Despite the fact that the altitude for this specific area does not have a high variability, since it is below 152 m, it is still important to include it as it is relevant to other areas and may influence the dengue disease in some way.

Figure 13 provides the altitude map with the dengue cases for October.

3.2.2. Housing Level of Overcrowding

In 87.88% of the AGEB, there is some level of overcrowding, and 7.81% had reported cases of dengue at some point in the study period. Of the 12.12% of the AGEB that do not have any level of overcrowding, 65.91% had no cases reported in the study period.

Where more AGEB are recorded, the overcrowding ranges between 40 and 50%. Of these AGEB, 99.19% presented cases of dengue. All of the 23.97% of AGEB where overcrowding is between 60 and 80% presented cases at some point in the study period.

Figure 14 shows the map by level of overcrowding and the cases for the month of October, which is the month with the most cases. On the map, it can be seen that the areas with the highest percentage of homes with some level of overcrowding coincide with the areas where there are the highest number of reported cases.

3.2.3. Population Density

In 19.83% of the AGEB, the population density was less than, or equal to, 2500 inhabitants per square kilometer. Of these AGEB, 33.33% did not present any reported dengue cases. In 25.90% of the AGEB, the population density was greater than 10,000 inhabitants per square kilometer, of which 97.87% had reported cases of dengue at some point in the study period.

Figure 15 presents the population density map with the dengue cases for October, showing that the areas with the highest population density coincide with the areas where there are more registered cases of dengue.

3.3. Statistical Analysis

Table 6. Descriptive statistics of the study variables.

Month	LST				NDVI				NDMI				Eggs by Positive Ovitrap				Dengue Cases
	Min.	Median	Mean	Max.	Min.	Median	Mean	Max.	Min.	Median	Mean	Max.	Min.	Median	Mean	Max.	Min.	Median	Mean	Max.
January	25.0	27.8	27.8	30.9	0.1	0.2	0.2	0.5	−0.1	0.1	0.1	0.1	7.9	27.5	29.2	77.8	0.0	2.0	2.5	16.0
February	26.7	31.2	31.2	34.8	0.1	0.2	0.2	0.4	−0.1	0.0	0.0	0.1	9.8	26.0	26.4	74.6	0.0	1.0	1.2	10.0
March	35.5	39.1	39.0	43.1	0.1	0.2	0.2	0.4	−0.1	0.0	0.0	0.1	9.4	30.9	33.1	90.0	0.0	0.0	0.8	5.0
April	38.9	42.4	42.5	46.6	0.1	0.2	0.2	0.4	−0.2	0.0	0.0	0.1	12.5	38.6	38.2	102.0	0.0	0.0	0.5	5.0
May	43.5	47.5	47.6	51.8	0.1	0.2	0.2	0.4	−0.2	0.0	0.0	0.1	13.9	33.8	36.1	84.2	0.0	0.0	0.6	5.0
June	45.0	48.8	48.9	53.8	0.1	0.2	0.2	0.4	−0.1	0.0	0.0	0.1	15.5	34.2	35.6	70.1	0.0	0.0	0.5	5.0
July	36.7	43.8	43.5	47.0	0.1	0.2	0.2	0.6	0.0	0.1	0.1	0.2	20.0	38.7	39.3	74.8	0.0	0.0	0.5	6.0
August	31.9	39.6	39.4	43.5	0.1	0.3	0.3	0.8	0.0	0.1	0.1	0.4	32.0	68.9	72.1	138.6	0.0	0.0	1.0	11.0
September	32.3	39.9	39.5	44.0	0.1	0.3	0.3	0.8	0.0	0.1	0.1	0.4	30.4	81.5	83.8	160.0	0.0	2.0	2.8	18.0
October	34.7	39.9	39.7	42.1	0.1	0.2	0.3	0.7	0.0	0.1	0.1	0.2	28.4	67.7	71.3	137.7	0.0	6.0	8.3	46.0
November	31.7	34.8	34.8	37.3	0.1	0.2	0.2	0.5	0.0	0.1	0.1	0.2	14.7	43.5	45.7	85.4	0.0	5.0	6.9	61.0
December	26.0	28.7	28.8	31.8	0.1	0.2	0.2	0.5	0.0	0.1	0.1	0.2	9.7	28.7	29.8	74.3	0.0	2.0	3.3	25.0

shows a summary of basic statistics for the monthly averages of the LST, NDVI, NDMI, eggs per positive ovitrap and dengue cases.

The maximum values for the NDMI and NDVI are 0.4 and 0.8, respectively, and occur between August and September. The maximum LST value is in June, and the minimum is in January. The average number of eggs per positive ovitrap in Culiacan is 2.2 times higher from August to October than the average for the rest of the year. The average number of probable cases of dengue in Culiacan is 5.5 times higher between October and November, compared to the rest of the year.

Spearman’s bivariate correlations were calculated for multicollinearity analysis. All VIFs were less than 10 as presented in

Table 7. Descriptive statistics of the study variables.

Model	Collinearity Statistics
	Tolerance	VIF
(Constant)
Altitude	0.744	1.354
LST	0.744	1.344
NDMI	0.303	3.295
NDVI	0.343	2.919
Dwellings with some level of overcrowding	0.599	1.670
Population density	0.647	1.545
Cases density	0.813	1.230
Eggs density	0.577	1.733

; this showed that there are no multicollinearity problems.

In the scatterplot of standardized forecasts on the x-axis and standardized residuals on the y-axis, no band is observed in the shape of the points. The assumption of equality of variances was not met; when a certain pattern was observed in the distribution of the data, there was a tendency towards a cone shape (Figure 16a). The Durbin–Watson statistic was reported as 1787, indicating that the residuals are independent. The assumption of normality was not met according to the graph of accumulated proportions of the expected and observed variable (Figure 16b).

As the assumption of linearity was not met, following the recommendation of Foley, the variable, cases, was transformed, for which, when presenting values equal to zero, the Log10 transformation (Y + 0.001) was used [39]. By transforming the variable Y (cases) using the formula Log10(Y + 0.001), partial regression graphs were obtained where a slight to moderate linear trend was shown, except for the case density variable, for which this variable, Y, was also additionally transformed as Log10 (density of cases + 0.001); the partial graphs for this model generally showed a linear trend (Figure 17).

The output of the model is shown in

Table 8. Statistical analysis results.

Model	R	R2	Corrected R2	Standard Error of the Estimate	Change Statistics					Durbin–Watson
					Change in R2	Change in F	gl1	gl2	Significance Change in F
1	0.893 a	0.797	0.797	0.182744495	0.797	1173.088	8	2386	0.000	1.927

^a. Predictor variables: (constant), log10 (density of cases + 0.001), population density, NDVI, altitude, egg density, dwellings with some level of overcrowding, LST, NDMI. Dependent variable: log10 (Cases + 0.001).

.

The scatterplot did not meet the assumption of homoscedasticity, and a pattern was observed in the data (Figure 18a). Similarly, the assumption of normality was not met, as the points deviate from the diagonal, and there are outliers reported (Figure 18b).

The existence of asymmetric variables with outliers meant that a robust regression was performed and, given the absence of homoscedasticity, the bootstrap method of simulation by resampling was used to estimate the standard error of the estimated parameters in order to account for the equality of variances [40]. The final model remained as Y = log10 (cases + 0.001), and the independent variables as altitude, LST, population density, log10(density of cases + 0.001), and egg density, given that the variables NDMI, NDVI, and percentage of dwellings with some level of overcrowding were not significant (

Table 9. Robust regression coefficients estimated using the bootstrap method.

Variable	Estimator	p-Value		Interpretation
Intercept	−1.00594	0.001
Altitude	−0.00209	0.007	*	−0.209
LST	0.00819	0.014	*	0.823
NDMI	0.28496	0.583		32.971
NDVI	0.19690	0.597		21.762
Dwellings with some level of overcrowding	−0.00068	0.590		−0.068
Population density	0.00002	0.000	*	0.002
Log10 (dengue cases + 0.001)	1.91386	0.000	*	1.914
Egg density	−0.00497	0.001	*	−0.496

* Significant variables.

).

3.4. Spatial Analysis

With the interpretation values of the significant variables obtained through the statistical analysis, the fundamental Saaty scale was applied to define the weights of each criterion and carry out the hierarchical analytical process [36]. The inconsistency index was 0.8636, for which the matrix is valid, the weights are shown in Figure 19.

Once all the spatial data layers were normalized, the weighted superposition method was applied to carry out the multi-criteria analysis and determine the areas at risk of dengue transmission.

3.5. Risk Areas

The dengue transmission risk areas obtained through the spatial model were classified into five strata according to the quintile method: 1, high; 2, medium high; 3, medium; 4, medium low; and 5, low. These areas will be explained in terms of the population found in the risk areas of the medium-high and high strata.

The AGEB that present risk areas are less than 10% of the total AGEB in the municipality. Figure 20 shows the risk areas for the months of August through to January, as follows: August, 0.06%; September, 3.91%; October, 52.91%; November, 41.27%; December, 1.48%; and January, 0.38%.

The total population in the risk areas was divided as follows: August, 12.53%; September, 12.53%; October, 26.69%; November, 31.03%; December, 19.66%; and January, 1.19%. The total percentage of cases that were registered in the risk areas were 0.26% in August, 3.79% in September, 51.92% in October, 36.76% in November, 6.55% in December and 0.73% in January.

4. Discussion

Many different techniques have been used to study the problem of dengue in Mexico and in the world; however, it remains an unresolved problem. As a multi-factorial public health problem, the identification of the factors that influence the incidence of dengue depends on the characteristics of each specific geographic area or region.

The use of remote sensing techniques applied to public health problems has been neglected and lacks substantive research. A systematic review in Colombia demonstrated this showing that although the use of these techniques in the investigation of vector-borne diseases has increased in recent years, their application in official control programs has been limited [41]. Remote sensing has not been exploited in Mexico; however, there have been attempts, such as those by the Ministry of Agriculture and Rural Development (SAGARPA), to use NDVI for phytosanitary epidemiological surveillance [42].

We found no evidence of the application of remote sensing techniques to guide vector control actions in Mexico and believe that this study could be the first of its kind.

Previous studies show that there is a strong relationship between rainfall and dengue transmission, in addition to temperature, humidity and land use. For example, Pathirana et al. found that dengue outbreaks were closely associated with post-rain periods, around two weeks after a major rainfall [43].

In our study area, the normalized difference vegetation and humidity indices were not significant. The behavior shown by the vegetation index is similar to that of the humidity index; however, the vegetation index presents more spatial variation. Nevertheless, the NDMI readings show that the humidity increased two months prior to an increase in dengue cases, that is, first, the humidity increases as a result of the rains that begin in June, followed by an increase in cases. This phenomenon has been identified in other studies, such as in Medellín, Colombia, which found that precipitation was positively correlated with the incidence of dengue, with a lag of 20 weeks [44]. The same was reported in another study in Anzoátegui, Venezuela, which found that, on average, the highest incidence of dengue occurred two months after the peak of precipitation [45].

The overcrowding variable, or the percentage of homes with some level of overcrowding and population density, was not significant for Culiacan. On the other hand, the population density variable, although it contributes very little to the risk areas, was significant at 3%. This coincides with the study carried out in Anzoátegui, Venezuela, where an analysis of the associated risk factors in cases of dengue infection showed no association with socioeconomic variables, with the exception of an association with the low level of the population’s knowledge about the disease [46,47]. The inclusion of the variable of houses with piped water, which can generate garbage dumps for mosquitoes, also showed no correlation with the level of dengue cases, despite the fact that a study carried out by Peña et al. found that the protection of water supply sources in homes was only fair to poor, and there were areas of landfills and micro-landfills [48].

In the present study, we found that the altitude variable was significant and contributed 6% to the generation of risk areas, which is an important consideration, despite the fact that no variability was found in any specific study area and that there are areas where this variable had a great impact on the proliferation of the dengue-transmitting mosquito. In addition to this, it is noted that over time the mosquito has learned to adapt to new conditions, as demonstrated by the epidemic outbreak in the state of Querétaro, where it was confirmed that Aedes aegypti had adapted to higher altitudes [49].

The variable of eggs per positive ovitrap contributed 18% to our risk areas. This variable shows a seasonal behavior that is aligned with the rainy season, which is notable in the study area, although this is not the case in other areas. One study identified that where there is constant rain throughout the year, producing constant humidity, there is no association of rainfall with the vector because there is no fluctuation, or defined seasonal pattern [50].

In the present study, the LST was used as the temperature variable and was identified as one of the most important variables for determining risk areas, with a 29% contribution. According to the literature, the proliferation of the dengue-transmitting mosquito occurs in temperature ranges from 15 to 40 °C [5]. In Culiacan, the average temperature is within this range, with the exception of the dry season in April to July where the average can reach up to 48.9 °C. When the temperature rises the humidity decreases and vice versa. Marques et al. reports that temperature influences the Aedes mosquito, from its development to its relationship with the virus, which makes it the most important climatic variable [51]. The marked seasonality of the rainy season in Culiacan coincides with an increase in dengue cases. This is verified by other studies, such as Barrera et al., who found that the Aedes aegypti population was driven by climate and human activities, and that peak mosquito density preceded the peak incidence of dengue during the rainy season [52].

At 33%, the case density variable made the highest contribution to the risk areas. This is undoubtedly a key determinant in the transmission of dengue, since the cases of infections, together with the presence of the vector, are potentially conducive to the appearance of dengue outbreaks. The probable cases that appear in areas with the vector are an important reason for the surveillance of this disease. This is known in Spain, where the virus does not currently circulate, but the regular arrival of infected people from endemic countries presents a constant risk. In addition, multiple studies in Mexico, and globally, include probable or confirmed cases as another study variable [53].

The determination of risk areas at the micro-regional level supports the fight against dengue. Each place has particular characteristics, and the actions of the authorities accord with the regional conditions. This is illustrated by another study in the city of Manta, Ecuador, which analyzed the variables associated with dengue transmission. This study considered the relationships with space, having the neighborhoods as the unit of analysis, and identified neighborhoods with the highest risk of transmission, to design focused strategies that considered the entomological and demographic data of each neighborhood [54,55].

5. Conclusions

In this study, LST, NDVI and NDMI information was generated from Landsat 8 satellite images through remote sensing techniques. These variables, together with altitude, population density, overcrowding, mosquito egg density due to positive ovitraps, and probable cases of dengue, were analyzed to generate transmission risk areas. The NDVI, NDMI, and overcrowding were not significant for the incidence of cases in Culiacan; however, all the other variables were. November showed the highest population living in the transmission risk areas; however, the transmission risk area itself was highest was October, thus reflecting the concentration of the population in specific areas.

One of the complexities of dengue is that it varies in each specific geographical area and it is not possible to apply the same control strategies or surveillance throughout the entire territory, thus it is necessary to know the conditions of each region, or area, in detail. The inclusion of remote sensing techniques in the dengue surveillance and control program identifies significant opportunities in the prevention and control of dengue and a potential solution to this public health problem.