Identiﬁcation of Risk Areas of Dengue Transmission in Culiacan, Mexico

: 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 signiﬁcant 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 signiﬁcant; 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 identiﬁed in October. In conclusion, it was possible to identify which of the study variables were signiﬁcant; in addition, monthly maps of risk areas of dengue transmission for Culiacan were obtained. Each geographical area had its own characteristics that inﬂuenced, in one way or another, the incidence of dengue, highlighting that the strategies for control of dengue must be speciﬁc to each region.


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 ISPRS Int. J. Geo-Inf. 2023, 12,221 3 of 26 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. Figure 1 shows the data management process for generating the spatial layers used to identify risk areas.

Materials and Methods
ISPRS Int. J. Geo-Inf. 2023, 12, x FOR PEER REVIEW 3 of 25 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. 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.

Sources of Information
Information sources were as follows: 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.

Sources of Information
Information sources were as follows:

Study Variables
The study variables were as follows: Percentage of inhabited private dwellings with some level of overcrowding; • Density of eggs per positive ovitrap; • Number and density of probable dengue cases.

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. 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 •

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: where Var is the processed image of the LST, NDVI or NDMI, i is the month (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(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 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 shows the information for the TIRS bands [25].

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: where Var is the processed image of the LST, NDVI or NDMI, i is the month (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(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 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 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: where LST is the land surface temperature in degrees Kelvin (K), C 0 to C 6 are the values of the coefficient for SW [26], TB 10 and TB 11 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: where Pλ is the spectral radiance value of the TOA (watts/(m 2 × 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: where T is the apparent brightness temperature in Kelvin (K); Lλ is the reflectance on the roof of the TOA atmosphere (watts/(m 2 × 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: 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 = NDV I − NDV I min NDV I max − NDV I min (6) where NDVI max is dense vegetation and NDVI min is the bare ground. The LSE is calculated with the following formula: The values for the emissivity of the soil and vegetation of each band are found in the Table 3. The concluding calculation corresponds to the application of the SW algorithm and the calculation of the average value of the LSE (ε): and the difference in LSE values (∆ε): Finally, the temperature is converted from degrees Kelvin to degrees Celsius with the following formula: Normalized Difference Moisture Index (NDMI) The normalized difference moisture index was calculated with the following formula: 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]: where z i corresponds to the values that are known, n corresponds to the number of values that entered the set search radius, d P i 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.

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]: 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).

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].

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: 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 1 , . . . , 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: where B 0 , 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].

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 shows the possible values of the scale ranging from 1 to 9. 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 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): 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: 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: The CR is accepted as long as it does not exceed the values of the Table 5. Table 5. Maximum percentages of the CR consistency ratio.

Matrix Size (n) Consistency Rate
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: (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: where Y is the dependent variable, C n are the coefficients or weights, and X n are the independent variables.

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. 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.

Land Surface Temperature
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: where Y is the dependent variable, Cn are the coefficients or weights, and Xn are the independent variables.

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. 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.   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.   The green areas of the urban zone were represented by an NDVI, ranging from 0. to 0.75, which concurs with other studies [38]. Of the entire study area, in April, the gree 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-u area or bare soil, in which 67.13% report no dengue cases. Figure 6b shows that in Octobe where 15.43% of the AGEB have urban tree vegetation, 55.36% reported cases of dengu 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. 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.

Normalized Difference Moisture Index
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.  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.  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.  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.  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

Dengue Cases
In the municipality of Culiacan in the period 2010-2020, there were 10,450 proba cases of dengue reported. Figure 11 shows that in the first quarter, a decrease in case noted, and the month with the lowest incidence of cases is April at 1.7%. The increase cases begins in August and rises until October with 30.4% of total cases. In Novemb cases decrease until December. 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. 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.

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.

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 Oct 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 this layer, the case density map for April shows that in 89.53% of the AGEB, there we reported dengue cases. Figure 12b shows that in October, 89.26% of the AGEB rep cases of dengue.
(a) (b) Figure 12. Cases density maps with dengue cases for April and October.

Spatial Variables
There is one raster image for each variable (altitude, population density and per age of homes with some level of overcrowding), which is three images in total. Thes iables have neither a spatial, nor a temporal, variability.

Altitude
Despite the fact that the altitude for this specific area does not have a high variab since it is below 152 m, it is still important to include it as it is relevant to other area may influence the dengue disease in some way. Figure 13 provides the altitude map with the dengue cases for October. 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. 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.

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.

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.

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.

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.

Housing Level of Overcrowding
In 87.88% of the AGEB, there is some level of overcrowding, and 7.81% had rep cases of dengue at some point in the study period. Of the 12.12% of the AGEB that d 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 these AGEB, 99.19% presented cases of dengue. All of the 23.97% of AGEB where crowding 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 mon October, which is the month with the most cases. On the map, it can be seen that the with the highest percentage of homes with some level of overcrowding coincide wi areas where there are the highest number of reported cases.

Population Density
In 19.83% of the AGEB, the population density was less than, or equal to, 2500 i

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.

Housing Level of Overcrowding
In 87.88% of the AGEB, there is some level of overcrowding, and 7.81% had repo cases of dengue at some point in the study period. Of the 12.12% of the AGEB that do 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% these AGEB, 99.19% presented cases of dengue. All of the 23.97% of AGEB where o crowding 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 mon October, which is the month with the most cases. On the map, it can be seen that the a with the highest percentage of homes with some level of overcrowding coincide with areas where there are the highest number of reported cases.

Population Density
In 19.83% of the AGEB, the population density was less than, or equal to, 2500 in itants per square kilometer. Of these AGEB, 33.33% did not present any reported den cases. In 25.90% of the AGEB, the population density was greater than 10,000 inhabi

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.   Table 6 shows a summary of basic statistics for the monthly averages of t NDVI, NDMI, eggs per positive ovitrap and dengue cases. The maximum values for the NDMI and NDVI are 0.4 and 0.8, respectively, an between August and September. The maximum LST value is in June, and the mini in January. The average number of eggs per positive ovitrap in Culiacan is 2.2 time from August to October than the average for the rest of the year. The average nu  Table 6 shows a summary of basic statistics for the monthly averages of the LST, NDVI, NDMI, eggs per positive ovitrap and dengue cases.

Statistical Analysis
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; 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). 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). 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.  The output of the model is shown in Table 8. 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, 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 output of the model is shown in Table 8. The scatterplot did not meet the assumption of homoscedasticity, and a pattern wa observed in the data (Figure 18a). Similarly, the assumption of normality was not met, a the points deviate from the diagonal, and there are outliers reported (Figure 18b). The existence of asymmetric variables with outliers meant that a robust regressio was performed and, given the absence of homoscedasticity, the bootstrap method of sim ulation by resampling was used to estimate the standard error of the estimated parameter 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, NDV 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).

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. and percentage of dwellings with some level of overcrowding were not significant (Table  9).

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.

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. 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.

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.

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

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].

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.