Spatial-Temporal Assessment of Environmental Factors Related to Dengue Outbreaks in São Paulo, Brazil.

Dengue fever, a disease caused by a vector-borne flavivirus, is endemic to tropical countries, but its occurrence has been reported worldwide. This study aimed to understand important factors contributing to the spatial and temporal patterns of dengue occurrence in São Paulo, the largest municipality of Brazil. The temporal assessment of dengue occurrence covered the 2011-2016 time period and was based on climatological data, such as the El Niño indices and time series statistical tools such as the continuous wavelet transformation. The spatial assessment used Landsat 8 data for years 2014-2016 to estimate land surface temperature and normalized indices for vegetation, urban areas, and leaf water. Results from a cross correlation for the temporal analysis found a relationship between the sea surface temperature anomalies index and the number of reported dengue cases in São Paulo (r = 0.5) with a lag of +29 (weeks) between the climatic event and the response on the dengue incidence. This relationship, initially nonlinear, became linear after correcting for the lag period. For the spatial assessment, the linear stepwise regression model detected a low relationship between dengue incidence and minimum surface temperature (r = 0.357) and no relationship with other environmental parameters. The poor relationship might be due to confounding effects of socioeconomic factors as these seem to influence the spatial dynamics of dengue incidence. More testing is needed to validate these methods in other locations. Nevertheless, we presented possible tools to be used for the improvement of dengue control programs. ©2019. The Authors.

Introduction

Approximately 390 million dengue infections occur globally every year (Bhatt et al., 2013). Dengue outbreaks have been reported in many countries in Asia and South America such as Brazil, Cambodia, Ecuador, Peru, Philippines, Singapore, and Thailand (Halsted, 2007), and it has been also reported in North America (CDC, 2010). In Brazil, dengue outbreaks have been reported since the 1980s, especially during the rainy season. However, in the last decade, the number of cases drastically increased (Aguiar et al., 2015).

Dengue is an acute febrile disease caused by a virus (Phillips, 2008), and it has four strains (Halsted, 2007). The most common carriers are Aedes mosquitoes, especially the Aedes aegypti. These mosquitoes are usually found in tropical urban areas, especially in places where there is rain water accumulation (Focks et al., 1995). In addition to rainwater, climatic characteristics such as air temperature, humidity, and solar radiation encourage Aedes mosquito to grow (Patz et al., 2005). Because of this, the interannual variability of climate correlates to dengue cases in tropical countries. El Niño indices like the El Niño–Southern Oscillation (ENSO), the Southern Oscillation Index (SOI), and the sea surface temperature anomalies (SSTA) also predict dengue outbreaks (Cazales et al., 2005; Gagnon et al., 2001; Fan et al., 2014; Fuller et al., 2009)

Other indicators of the spread, establishment, and persistence of vector‐borne diseases such as dengue are environmental (Messina et al., 2012). Phillips (2008) showed that the expansion of urban areas is a key driver for the increased endemicity of dengue. Since the urban population is increasing worldwide, dengue cases will grow as the trend continues (Messina et al., 2012). Other environmental factors that predict dengue include terrain elevation and vegetation (Fuller et al., 2009; Moreno‐Madriñán et al., 2014; Stanford et al., 2016).

The variety of factors that lead to rises in dengue cases calls for the development of tools that predict dengue based on its geographical distribution in the past (Messina et al., 2012). For other environmental applications, the combinations among the spatial distribution of environment parameters are used to perform risk assessments through statistical or mechanical modeling of environmental trends. A growing alternative to monitor environmental trends, especially in resource‐limited regions, is the use of remote sensing. Hadjimitsis and Clayton (2009) used remote sensing products as predictors for the computation of such models, exploiting advantages such as (1) the synoptic view of satellite or airborne images, which allows the acquisition of data from the entire environment; (2) information from remote and sometimes inaccessible regions; and (3) historical data.

One of the satellite‐derived environmental variables that dengue outbreaks modeling uses most is rainfall, usually from the Tropical Rainfall Measuring Mission satellite (Ashby et al., 2017; Moreno‐Madriñán et al., 2014; Pinto et al., 2011). However, Stanforth et al. (2016) showed that while temperature, elevation, and vegetation cover significantly correlated with dengue outbreaks reported in Magdalena River watershed in Colombia, rainfall did not. Satellite‐derived products such as land surface temperature (LST), digital elevation models, and vegetation indices can be used to compute temperature, elevation, and vegetation cover. However, since the Moderate Resolution Imaging Spectrometer products—usually used for dengue research (Ashby et al., 2017; Buczak et al., 2012; Stanforth et al., 2016)—do not model these factors spatially, they do not model the spatial distribution of dengue cases accurately in the Magdalene River watershed (Stanforth et al., 2016). This calls for a remote sensing approach that can provide a better spatial understanding of the environmental variables.

The goal of this research was to identify the temporal and spatial characteristics of environmental factors affecting the occurrence of dengue in a large urban area, São Paulo, the largest municipality in Brazil. To do this, we evaluated environmental and dengue cases data from the city to determine the temporal and spatial relationships among dengue cases and environmental factors. The evaluations were conducted based on a time series analysis of the data set as well as remote sensing‐derived indices to measure the geographic distribution of environmental factors.

Material and Methods

Study Site

At approximately 11 million inhabitants, São Paulo is one of the most populous cities in the world (Prefeitura do Município de São Paulo, 2017). This population is distributed within 1,509 km2 of area which consists of 32 subprefectures or 96 districts, as presented in Figure 1. São Paulo is Brazil's major financial and economic center, contributing 11.37% of total Brazilian gross domestic product (IBGE, 2012). However, São Paulo has broad economic diversity, with 13.30% of households living on less than two minimum wages (MWs) jobs, 24.39% living on wages amounting to 2–5 MWs, 25.97% on 5–10 MWs, 11.29% on 10–15 MWs, 10.98% on 15–25 MWs, and 14.06% of families living on more than 25 MWs (SEADE, 2000). The use of MWs as an economic diversity indicator is a distinctive feature of the Brazilian economy. This status is based on the fact that MW affects employment through wages, pensions, benefits, inflation, and the public deficit (Lemos, 2004). Moreover, the Brazilian Institute of Geography and Statistic (IBGE in Portuguese) uses MWs to define the social class of the Brazilian population. Therefore, families living on low MWs represent a lower social class, while a higher MWs represents richer families. This condition along with the urban landscape diversity, the climatic context, and the size of the exposed population make São Paulo a good study site for the current research.

Figure 1

Location of São Paulo Municipality within Brazil and its 96 districts.

Dengue Data

Confirmed dengue cases for the municipality of São Paulo were acquired on a monthly basis data set from the State Secretariat of Health for the period between 2011 and August of 2017 (São Paulo State Secretariat of Health, 2017). For each district of São Paulo municipality, only data of reported dengue cases were collected from the City Secretariat of Health for the period between 2011 and August of 2017 (São Paulo Secretariat of Health, 2017).

Climatological Data

Daily temperature and precipitation data from the weather station of Mirante do Santana (23°30′S,46°37′W) were acquired from the National Institute of Meteorology (INMET in Portuguese, 2017) via the Historical Research and Teaching Database. Figure 2 summarizes the data set for the same period of dengue data. It was observed that higher temperature and higher precipitation occurs in the beginning of each year, which is related to the austral summer. Therefore, São Paulo's climate can be classified as humid, subtropical, and oceanic, without a dry season and with hot summers (classified as Cfa based on Köppen's classification method; Alvares et al., 2014).

Figure 2

Daily time series (2011–2016) of precipitation (mm) and temperature (°C) from the weather station located at “Mirante do Santana,” São Paulo, Brazil.

The most common El Niño indices are the air pressure indices, sea surface temperature indices, and longwave radiation indices. Since sea surface temperature data were recognized to be a key player in ENSO (Rasmussen & Carpenter, 1982), the use of this type of index is more common in the scientific literature. To calculate the sea surface temperature, certain regions were defined for in situ measurements and were divided into four areas: Niño1, Niño2 (combined into Niño1 + 2), Niño3, and Niño4. However, it was identified that a combination between areas 3 and 4 would be the most ENSO‐representative (Barnston et al., 1997). Huang et al. (2013) showed that the area‐averaged Niño‐3.4 SSTA index in the tropical Pacific is more consistent with in situ observations in extended reconstructed sea surface temperature, version 3b. Thus, in this study we decided to use the Niño 3.4 SSTA index. The weekly SSTA index data for Niño region 3.4 (the area within 5°N, 5°S, 120°W) were taken from the Climate Prediction Center of the U.S. National Weather Service (2017). The weekly data set was chosen to match the epidemiological weeks.

Satellite Data

Earth Observation Images

Landsat 8 images from the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS) sensors were used to remotely estimate environmental variables to complete the spatial analysis of this study. TIRS images were acquired around 13:00 UTC time (10:00 in São Paulo) and were used at level 1 processing products to compute LST, while OLI atmospherically corrected products—Landsat 8 Surface Reflectance (L8SR)—were used for the computation of environmental indices. Both products were downloaded from the U.S. Geological Survey (USGS) Earth Explorer website, and Table 1 presents the list of images used.

Table 1

Image Identifier, WRS‐2 Path and Row, and Date for the OLI Images Used in This Study

Image identifier	Path/row	Acquisition date
LC82190762014039LGN00	219/76	2014‐02‐08
LC82190772014039LGN00	219/77	2014‐02‐08
LC82190762015010LGN00	219/76	2015‐01‐10
LC82190772015010LGN00	219/77	2015‐01‐10
LC82190762016109LGN00	219/76	2016‐04‐18
LC82190772016109LGN00	219/77	2016‐04‐18

Note. Dates are formatted as year‐month‐day. OLI = Operational Land Imager; WRS = Worldwide Reference System.

Environmental Indices

Environmental indices such as the Normalized Difference Vegetation Index (NDVI; Rouse et al., 1973), Normalized Difference Water Index (NDWI; Gao, 1996) and Normalized Difference Built‐up Index (NDBI; Zha et al., 2003) were estimated from OLI‐SR to characterize the environment using remote sensing images. NDVI is one of the most used vegetation indices especially because of its advantages such as lower influence of atmospheric variations, greater sensitivity to chlorophyll, the reduction of noise by normalization between −1 and +1 (Gutman, 1991; Luo & Li, 2014), and the possibility to monitor seasonal changes (Ponzoni & Shimabukuro, 2009). The NDWI is also known as the leaf area water‐absent index and is used primarily to estimate the water content within vegetation (Gao, 1996). However, recent studies have been using NDWI as an index to identify water bodies (McFeeters, 1996; Moknatian et al., 2017; Xu, 2006). Zha et al. (2003) developed the NDBI index to be used to highlight urbanized and barren soil areas.

Estimation of LST

The retrieval of LST from Landsat 8 was based on the measurement from the TIRS. To estimate LST from the TIRS images, we used the single‐channel method, which has been successfully applied in several studies (Jimenez‐Munoz et al., 2009; Jimenez‐Munoz & Sobrino, 2003; Sobrino et al., 2004). To compute LST for Landsat 8/TIRS, we used equation (1), which Jimenez‐Munoz et al. (2014) proposed: where LST is given in Celsius (°C), ε is the spectral surface emissivity, L is the top‐of‐the‐atmosphere (TOA) radiance, γ is a parameter computed by equation (2), δ is a parameter computed by equation (3), and ψ 1, ψ 2, and ψ 3 are derived from the water vapor (w) content following the procedure described in Yu et al. (2014, equation 6. where T is the at‐sensor brightness temperature, b is 1324 for Band 10, and 1199 K for Band 11.

T is calculated based on the conversion from digital numbers to TOA radiance (equation (4)), and then from TOA radiance to T (equation (5)) using the thermal constants provided in the metadata file of the Landsat 8/TIRS image (USGS, 2016). where L is the TOA spectral radiance, M represents the band‐specific multiplicative rescaling factor, Q cal is the Band 10 image, and A is the band‐specific additive rescaling factor. where K 1 and K 2 stand for the band‐specific thermal conversion constants from the metadata.

The atmospheric functions ψ 1, ψ 2, and ψ 3 were computed from water vapor (w) measurements by the Moderate Resolution Imaging Spectrometer Atmospheric Profiles product (MOD07) and the following procedure proposed by Yu et al. (2014, equation 6):

ε was estimated following the procedures described in Sobrino et al. (2004), an NDVI‐based method for the estimation of ε. This study uses an NDVI threshold method (NDVITHM) for ε estimation, using calibration values presented for Landsat 8/TIRS band 10 (Skoković et al., 2014), as shown in equation (7). where ρ red is the reflectance at the red spectral band and P is the proportion of vegetation (also known as fractional vegetation cover). P can be calculated as presented in equation (8).

Where: NDVIv is a NDVI value from a vegetation pixel and NDVIs is a NDVI value from a soil pixel. In this study we used the values Sobrino and Raissouni (2000) proposed, 0.5 and 0.2, respectively, for NDVIv and NDVIs.

Continuous Wavelet Transformation

Wavelet analysis involves transformation of a data series with a “mother” wavelet to produce “daughter” wavelets (Wiebe & Sturman, 2011). Thus, data are transformed into the frequency domain, which highlights the spectral characteristics of a time series. Johansson et al. (2009) presented an example of application of wavelet transformation for a dengue incidence time series and stated that the transformation was useful to differentiate multiannual patterns that made seasonal variation clear. In another study that applied wavelet analysis to dengue cases, Nagao and Koelle (2008) showed that there is a shift in the frequency dengue cases in Thailand. According to Johansson et al. (2009) there are two advantages of using wavelet technique for assessing the relationship between ENSO and dengue: (1) It allows the separation by time scale, and (2) it provides a way to measure nonstationary association. Thus, the wavelet technique can be considered as a cross correlation among a time signal and a set of wavelets.

In this study we used the continuous wavelet transformation to analyze how the frequency content of the weekly SSAT index data changes over time. The scale for this time series (n = 305) was set to 64. The computation of the wavelet and coefficients was performed using the Matlab wavelet toolbox (MathWorks, 2017). Compared with the Morlet wavelength, Wiebe and Sturman (2011), found the Haar wavelet to be more suitable for this type of analysis; hence, it was chosen for the temporal component of this study.

Temporal and Spatial Analysis

The above data sets were integrated to conduct the temporal and spatial assessments. For the temporal analysis we tried to identify a relationship between dengue cases and climatological variables. Since those are the variables with a high temporal variability, we complemented the relationship by applying the continuous wavelet transformation to the climatological data. To verify the benefit of this process, we applied the Pearson correlation and cross‐correlation analysis to the data with and without the continuous wavelet transformation. The temporal analysis provides information about the most vulnerable period for dengue outbreaks, helps identify coincidents events that could be their most important determinant variables, and can also help estimate the intensity of future dengue outbreaks. Parallely, remote sensing data such as LST and environmental indices and the spatial distribution of dengue cases per district were used for the spatial analysis. This analysis will provide a sense of where most of dengue cases are located as well as the environmental factors that correlate with high rates of dengue.

Results

Temporal Analysis

Figure 3 shows the monthly fluctuation of confirmed dengue cases along the studied years, with most of the dengue cases in São Paulo occurring at the end of the austral summer, during the first half of the year, when precipitation and temperature are the highest (as shown in Figure 2), but variable annual incidence range from as low as 475 cases in 2016 to 43,359 cases in the previous year (2015).

Figure 3

Distribution of confirmed dengue cases (in log scale) per epidemiological months in the municipality of São Paulo between the years of 2011 and 2017.

Figure 4 presents the reported dengue cases with the climatic variables (precipation and temperature). No correlation was detected between the climatic variables (precipitation and temperature) and number of dengue cases. Correlation analysis with temperature resulted in a Spearman's rank coefficient (rs) value of 0.11 (p = 0.33) and a linear Pearson coefficient (r) value of 0.06 (p = 0.59). Similarly, a rs value of 0.09 (p = 0.40) and a r value of 0.10 (p = 0.36) were the correlation analysis found with precipitation.

Figure 4

Time series of precipitation, temperature, and confirmed dengue cases (from 2011 to 2016).

Dengue cases by month in São Paulo city were plotted with the monthly SSTA index (Figure 5a). Figure 5a shows a pattern with peaks of dengue incidence within the first semester of each year (also observed on Figure 3). Increasing trends of SSTA coincided with dengue peaks in years 2011–2015, although a second and lower dengue peak occurred in 2013 that coincided with a SSTA peak. Conversely, after the highest SSTA peak of the study period (right at the beginning of 2016), there were almost imperceptible dengue peaks in years 2016 and 2017. To further explore the assessment of the relationship between weekly SSTA index and dengue cases, the Haar continuous wavelet transformation was applied to the weekly SSTA index time series and depicted in Figure 5b. The use of the Haar wavelet transformation (Figure 5b) suggests a relationship between coefficients (scalogram) and dengue cases (grey line). The scalogram for the weekly SSTA index shows increasing coefficient values tending to be blue with a tendency to green or red color for decreasing SSTA values. Decreasing coefficient values tend to happen at the second semester of the year right before the dengue season of the following year, which occurs at the first semester of the year (austral summer and austral autumn) during increasing coefficient values. However, 2016 and 2017 did not maintain this trend presumably in connection with the fact that the preceding year, 2015 (highest dengue peak of the study period) presented El Niño conditions during the entire year thus coinciding with a steady steep increase in SSTA (blue scalogram) and consequent drastic SSTA decrease (red scalogram) in the following year (2016). This year presented a drop in dengue incidence from 45,359 to 475 confirmed dengue cases, and the following year (2017) was even lower. The two years with the above normal peaks in dengue incidence (2014 and 2015) were predominantly El Niño years (NOAA, 2019) and presented an overall increasing trend for the monthly SSTA index with a predominant lower coefficient value for the Haar continuous wavelet transformation (<0, blue dominating color in the scalogram).

Figure 5

(a) Monthly SSTA index and dengue cases plot, and (b) wavelet transformation of the monthly SSTA index and dengue cases plot. SSTA = sea surface temperature anomalies.

Our observations suggest that the continuous wavelet transformation of the SSTA index has a potential predictive ability over time that could complement the SSTA index. In terms of wavelet transformation the acute peaks in dengue cases appear to occur when coefficient values are 0 or negative (≤0, blue) but preceded by negative or slightly positive coefficient values (≤1, blue or olive green) the preceding year. Notice that a highly positive coefficient value (>2, light green or red) either during the normal dengue season (the first semester of the year) or at some point during the preceding year would result also in a low peak of dengue cases. Thus, lower wavelength transformation values tend to be associated with more dengue incidence. On the assumption that these described patterns for the SSTA curve and the coefficient values of the wavelength transformation would be a constant, it could be expected that 2017 would have few cases. This is because the last part of 2016 shows a strong decrease in the SSTA trend line and the wavelength transformation showed highly positive coefficient values (>2, light green or red). In fact, Figure 5 shows that 2017's dengue season did not present important peaks in confirmed dengue cases in São Paulo municipality. Only 214 cases were confirmed during the first eight months of 2017, while 2014 and 2015 had 30,066 and 45,359 cases, respectively, in the same corresponding period. The year 2016 differed with years 2011–2015 in the drastically higher coefficient value of wavelength transformation and lower monthly dengue incidence.

Other than our suggestion for a possible influence of ENSO measured through SSTA and wavelength transformation, we could not find a clear proposition to explain such a decrease. Although preventive measures sponsored by the public health authorities have been suggested as the cause of the decline in dengue incidence (Santiago, 2016), it is not clear why such supposedly measures did not protect against the huge Zika virus outbreak that occurred during the same year and whose etiologic agent is transmitted by the same vector, likewise why the peak in the Zika virus epidemic took place in 2016 despite the highly positive coefficient value in wavelet transformation. In fact, the decline in dengue incidence after 2015 could have also been influenced by the Zika virus epidemics (Fuller et al., 2017). The temporal coincidence between the two viruses might have cause some sort of competition that would have caused only one of the two to peak at a time. On the other hand, the reason why the Zika epidemic peaked in 2016 (and not in 2015 as in the dengue peak) might have been in part due to the timing of the Zika virus arrival into the naïve population and the time period it took to spread. For instance, we do not know if had the Zika virus arrived earlier its peak might had been even stronger due to even more favorable El Niño conditions and it might have coincided temporally with that of dengue or conversely, if one of the two would have just inhibited the other. In any case, the causes of this dengue decline are still not fully understood (Lopes et al., 2018).

We compared the continuous wavelet transformation scale values and SSTA index to the number of dengue cases via the Pearson correlation coefficient (r) and rs. Figure 6 shows the absolute r and rs values (black and red columns, respectively) between dengue cases and the SSTA index (first two columns) and between dengue cases and each scale from the continuous wavelet transformation (remaining columns). In both relationships the rs values were higher than the r values; however, for the continuos wavelet transformation scales this difference is minimal at higher scales. This indicates that, while SSTA was not linearly related to confirmed dengue cases, the transformed data were. Furthermore, the correlation between dengue cases and the SSTA was weaker (rs = 0.14) than that between dengue cases and the scale of the continuous wavelet transformation (rs = 0.38, at scale 32) when using a nonlinear relationship estimator. For a linear correlation, the transformated data was stronger (r = 0.37, at scale 30) than the SSTA (r = 0.05).

Figure 6

Absolute Pearson's correlation coefficient and absolute Spearman's rank correlation coefficient for the relationship between denge cases and both SSTA index and continuous wavelet transformation scales. SSTA = sea surface temperature anomalies.

Spatial Analysis

Surface Temperature of São Paulo

As expected, the urban area of São Paulo has a much higher LST than the suburbs; it is a major example of an urban heat island (UHI) in Brazil (Barros & Lombardo, 2016; Lombardo, 1985). Across dates (January 2015, February 2014, and April 2016) and different temperature scales, the spatial pattern in the images is the same, with higher temperatures in urban areas and lower temperatures in the suburbs (Figure 7). Barros and Lombardo (2016) recently discussed this spatial distribution of LST, as they identified different intensities of UHI occurring within São Paulo city. The authors also observed the presence of a park cool island that partially offsets the UHI effect in neighborhoods located in the midwest and south of the area that have more space with vegetation cover. Zipper et al. (2016) observed a similar effect in Madison, WI, USA.

Figure 7

Land surface temperature for São Paulo city (a) for 2014, (b) for 2015, and (c) for 2016.

Spectral Indices of São Paulo City

Spectral indices (NDBI, NDVI, and NDWI) were computed for São Paulo municipality to explore the spatial distribution of dengue cases and remote estimated environmental indices. NDBI computation for São Paulo (Figure S1, see in the supporting information) showed a similar spatial pattern as LST results (Figure 7), with higher index values for urban areas and lower for the suburbs and water, respectively. This was expected since NDBI is an index to quantify the amount of building area, and it did not show a significant difference among 2014 (Figure S1a, see in the supporting information), 2015 (Figure S1b, see in the supporting information), and 2016 (Figure S1c, see in the supporting information). The Kruskal‐Wallis test for equal medians for these three data sets showed that there is no significant difference between medians (H = 4.363, p = 0.11). The NDVI showed a significant variation (Kruskal‐Wallis test, H = 6.445, p = 0.04) on the values computed over continental areas—not including lake areas (see Figure S2 in the supporting information). Although the NDVI calculated from 2015 (Figure S2b; see in the supporting information) and that from 2016 (Figure S2b, see in the supporting information) were similar, these were different to the NDVI calculated from 2014 (Figure S2a, see in the supporting information), especially over aquatic systems. These changes on NDVI values could be related to the 2014–2015 drought, which lowered the water level in all of São Paulo municipality reservoirs, causing problems such as eutrophication and algal blooms (which increases the vegetation signature within the water surface). The same could be observed for the results from the computation of the NDWI (Figure S3, see in the supporting information). In these results, for 2014 (Figure S3a, see in the supporting information) and 2015 (Figure S3b, see in the supporting information) the water bodies showed a medium value (close to the values observed for vegetation), while the NDWI values only increased in 2016 (Figure S3c, see in the supporting information). The spatial pattern of NDVI was the opposite of LST and NDBI, with lower values for water bodies and urban areas and higher values for the suburbs, while for the NDWI higher values were observed over water bodies and vegetation (in the suburbs) and lower values were observed in the urban area.

Spatial Analysis

To evaluate the spatial characteristics of dengue outbreaks in São Paulo municipality, a data set of reported cases of dengue per district was used (São Paulo Secretariat of Health, 2017). Unfortunately, there is no data of confirmed cases at the district level; thus, reported cases were plotted per district per 10,000 habitants (Figure 8). The intensity observed in the temporal assessment of dengue cases in São Paulo was also observed in this spatial plot of reported cases, which showed more reported cases in 2015 (Figure 8b) and fewer in 2016 (Figure 8c). For 2014 (Figure 8a), the districts in the northern and northwestern parts of São Paulo municipality had the highest number of reported cases of dengue. This spatial pattern intensified in 2015 (Figure 8b). In 2016, however (Figure 8c), the eastern city districts had the highest number of reported cases.

Figure 8

Reported dengue cases per 10,000 habitants per district of São Paulo city (a) for 2014, (b) for 2015, and (c) for 2016.

Comparing the spatial distribution of reported dengue cases (Figure 8) to the mapping of LST (Figure 7) and that of NDBI (Figure S1, see in the supporting information) in São Paulo revealed that some of the districts with high dengue incidence had lower LST and NDBI than the main urban area. These spatial patterns could have been not expected since the primary vector (A. aegypti) is traditionally associated with warm temperatures, as does the dengue virus incubation (Messina et al., 2012), and with more urbanization, since urbanization is an important factor for the increase of endemicity of dengue (Gubler, 2004; Phillips, 2008; Russell et al., 2009). However, negative correlations with dengue incidence could occur with temperaures higher than the comfort zone of A. aegypti and conditions such as UHI could contribute to those higher temperatures and also explain negative correlations with urbanization. Such apparent discrepancy however, was not the case in all the districts with high dengue incidence as some others presented the expected coincidence of high dengue incidence with high LST and NDBI. Inconsistent patterns were also observed for the NDVI (see Figure S2 in the supporting information) and NDWI (Figure S3, see in the supporting information) as higher values were found for the northern and southern parts of São Paulo municipality. Yet the number of reported cases was higher just in few sections of the northern part of the municipality and lower in the southern part.

Correlation between average values of environmental indices per district and number of dengue cases per district were computed for the three years. The statistical analysis showed that dengue incidence was slightly correlated with LST (r = −0.25, p value < 0.001) but not with NDBI (r = −0.07, p value = 0.18), NDVI (r = −0.09, p value < 0.09) and NDWI (r = −0.075, p value < 0.20).

Based on these results presented, none of the selected spatial environmental factors seemed to clearly influence the spatial distribution of dengue cases. Therefore, to better evaluate the environmental variables, a forward stepwise regression was applied for the four remotely sensed estimated environmental variables (LST, NDBI, NDVI, and NDWI). For each variable, statistical estimators (minimum, maximum, average, and standard deviation) were calculated, using the variability within the area from each district, for a total of 16 variables. Table 2 summarizes the result of the forward stepwise regression, showing that minimum temperature is the environmental variable with the greatest r value possible (0.357) and the greatest determination coefficient (R 2) possible (0.127) with a standard error estimate of 77.45. Hence, minimum surface temperature was identified as an important variable for modeling the number of dengue cases. All other variables failed the F test and were not selected from the stepwise regression because of their high p values (Table S1, see in the supporting information).

Table 2

Summary of the Stepwise Regression for the 16 Environmental Variables Analyzed for Each District and Coefficients of the Model Showing Minimum Temperature as the Best Variable

Step #	Variables entered	Variables removed	r	R 2	p value	Coefficients
	Constant					115.619
1	Temperature minimum		0.357	0.127	<0.001	−3.254

Discussion

Temporal Assessment

The increase in the SSTA index between 2014 and 2015 affected São Paulo and caused the worst drought in the city since 1930 (Gutiérrez et al., 2014). During this period, the southeastern region of Brazil suffered from a 2‐year‐long drought, and São Paulo city was particularly affected. Escobar (2015) attributed this atmospheric anomaly to an atmospheric pachyderm which blocked the South Atlantic Convergence Zone, causing blockage of the southward flow of the water vapor from the Amazon Rain Forest from the northern region to the southeastern region of the country. Since the water vapor could not migrate to the southeastern region, precipitation levels in the peak of the summer were half the average in 2014 to 2015.

Concurrently, the number of confirmed dengue cases drastically increased in 2014 and 2015 (Figures 3, 4, 5) despite the lower than normal precipitation during this period. Yet, the normal intra‐annual increase in dengue incidence tends to coincide with that of precipitation, which in Sao Paulo occurs during the first semester of the year (Figures 2 and 3). A suggestion to explain similar inconsistency in the watershed of the Magdalena River in Colombia has been proposed by Eastin et al. (2014), Stanforth et al. (2016), and Ashby et al. (2017) who found dengue incidence to be importantly associated with temperature but not with precipitation in interanual trends. Interestingly, incidence is higher during the rainy months within the year, which are also the warmest, both in the Magdalena River watershed and in Sao Paulo (Figure 2). The authors suggest that temperature may be the important driver since human‐filled water containers make mosquito reproduction less dependent on precipitation, which is just associated without causation. Similar discussion has been proposed for other places (Moore et al., 1978; Pontes et al., 2000; Scott et al., 2000). In fact, Sao Paulo experienced slightly higher than normal temperature during the dengue outbreaks of 2014 and 2015. It is important to note that the end of 2014 and the entire 2015 presented El Niño conditions (NOAA, 2019).

In Colombia, higher incidence of dengue and malaria have been reported during El Niño years, which are characterized in this country by below normal precipitation and above normal temperature (Eastin et al., 2014; Gagnon et al., 2001; Kovats, 2000; Poveda et al., 2000). Similarly, higher incidences during drier years have been observed in Thailand (Ungchusak & Kunasol, 1988). Furthermore, Wu et al. (2018) did not find sufficient evidence to support a significant effect of precipitation on dengue fever. Beyond water storage, Cazales et al. (2005) observed that during drought periods, fast‐flowing rivers can recede into a series of stagnant pools, which are ideal for mosquito breeding. Thus, the positive increase in the SSTA index (Figure 5a) or the lower coefficient values (Figure 5b) appear to correlate with the higher number of dengue cases in São Paulo municipality, even during drought conditions. This same pattern was observed by Fan et al. (2014); however, they used the SOI to correlate with the dengue incidence in the Guangdong Province, China. Because of the difference in index as well as geographic location, the authors observed that negative values of SOI were related to warming caused by El Niño. These results suggest that climate variability will affect the potential for dengue epidemics, especially the increase of temperature and consequently drought conditions, which many regions experience in El Niño years.

It is important to highlight that the large dengue outbreaks of 2014 and 2015 occurred despite the implementation of the Brazilian National Plan for Dengue Control and the São Paulo Municipal Dengue Control Program, which started in 2002 and 2005, respectively (Taliberti & Zucchi, 2010). These programs were based on traditional interventions that rely on community visibility alone but not aided by temporal and spatial forecasting methods that could identify when and where control efforts should be made. The consideration of the SSTA index and the continuous wavelet transformation promises to be a helpful tool to observe the relationship between SSTA index (or El Niño effect) and reported dengue cases and thus be a potential pathway of research toward achieving an effective dengue incidence forcasting tool for public health management and policy makers.

These results align with other studies that assess the ENSO‐dengue relationships through other means and in other geographical areas, by identifying a lag between an El Niño event and a peak in dengue incidence (Eastin et al., 2014; Tipayamongkholgul et al., 2009). However, the effect of El Niño to the climate in São Paulo is not yet well understood. Grimm et al. (2000) did not find a significant effect of El Niño and La Niña cycles in the state of São Paulo. Contrarily, Rogers (1988) observed that in the vicinity of São Paulo the ENSO was responsible for the dry period between October and March. In addition, Aceituno (1988) observed that the Southern Oscillation and rainfall in the southeast of Brazil are negatively correlated, which corroborates Rogers' (1988) observations. Nevertheless, the opposite was also observed when El Niño years were related to the increase in rainfall in the period between November and January. All these mixed observations underline the lack of understanding of ENSO effects on the climate of São Paulo, thus highlighting the need for a better characterization of such effects in southeastern Brazil.

Spatial Assessment

Based on the visual interpretation of remotely estimated environmental variables, this study suggests that urban and rural districts have lower dengue incidences than suburban districts (located between urban areas and rural areas). In the suburban districts, a variance in the spatialized LST (Figure 7), NDBI (Figure S1), NDVI (Figure S2), and NDWI (Figure S3) was observed based on the standard deviation values from each district. Thus, the number of dengue occurrences could be related to the variability of environmental parameters. This could also help explain the spatial difference in the number of dengue cases between the southern and northern parts of the municipality. In the southern (lower standard deviation) part of São Paulo environmental factors varied less and the number of dengue cases was lower (Figure 8).

The overall low correlation between environmental variables and the number of reported dengue cases might be due to a complex interaction of counfonding effects of human hosts and their environment, which may have not been considered in this study. Although the annual temperature range in our study site roughly fluctuated between 12 and 32 °C (Figure 2), which includes the optimal survival for the vector, around 21 °C (Brady et al., 2013), it exceeded the lower limit for midrange temperatures of approximately 15–35 °C (Brady et al., 2013) and the optimal diurnal temperature range of 18–32 °C for mosquito survival and virus transmission (Eastin et al., 2014). Temperatures lower than 18 °C happened during midyear right after the dengue season, which strongly suggests temperature as a key factor determining the usual dengue season during the first half of the year.

The amount of daily hours with extreme air temperatures outside the optimal 18–32 °C range is being increasingly viewed as a potential environmental factor in recent ecological dengue studies (Eastin et al., 2014; Lambrechts et al., 2011). The 1,509 km2 of our study site could be enough to detect important spatial differences in the relationship between dengue incidence and diurnal temperature range among its 96 composing districts (between urban and suburban). Furthermore, the cited studies used air temperature, which exhibits lower daytime values and a much smaller diurnal range than the satellite‐derived LST used in this study. This adds to the uncertainty resulting from interchangeably using either one of these variables in urban areas despite their commonly well‐stablished relationship (Schwarz et al., 2012). Yet, a likely more important confounding factor interfering with the spatial assessment may be socioeconomic influences (Banu et al., 2011; Hayden et al., 2010). In fact, socioeconomic factors have been suggested to be important determinants in the incidence of diseases caused by viruses transmitted by A. aegypti (Moreno‐Madriñán & Turell, 2017; Reiter, 2001; Vasconcelos et al., 1998, 1999). A socioeconomic approach was previously used to derive a risk assessment for dengue in Cali, Colombia (Hagenlocher et al., 2013).

To explore the possibility that the poor relationship between dengue incidence and environmental variables might have been caused by socioeconomics, we conducted a complementary analysis where the socioeconomic classes were related to the number of reported dengue cases (Table 3). Dengue cases were distributed according to the reported residences's social‐economic class (percentage of families living within an X number of MWs) per district measured by the State System of Data Analysis (Fundação Sistema Estadual de Análise de Dados, SEADE, in Portuguese). Those districts with more families having more than 25 MWs have less dengue cases. This finding was expected since families with higher wages are more likely to live in wealthier districts, which may have better mitigation programs, sanitation, and living infrastructure. On the other hand, for families living on less than two MWs, the values of the r were close to 0 during the three analyzed years. This unexpected absence of relationship might be explained by confounding factors such as underreporting by people from the lowest socioeconomic class, where there is greater proportion of informal jobs, and thus less incentive to report sickness as compared to employees in the formal economy where workers need to present an excuse to their employers in order to get sick leave. Families living on 2–15 MWs (middle class), especially those in the middle range of 5–10 MWs, showed the most positive relationship with dengue incidence (r = >0.27, p < 0.01). Some differences to families living in richer neighborhoods may be based on the better vector control programs and living conditions of those neighborhoods, which reduce the opportunity of contact between humans and vector (Moreno‐Madriñán & Turell, 2017), and more common underreporting in poorer neighborhoods (Eastin et al., 2014).

Table 3

The Correlation Between Percentage of Residence in Each Economic Class (Based on MW) per District and Reported Dengue Cases

Year	Less than 2 MW	2 to 5 MW	5 to 10 MW	10 to 15 MW	15 to 25 MW	More than 25 MW
2014	−0.004 (p = 0.96)	0.075 (p = 0.46)	0.095 (p = 0.35)	0.000 (p = 0.99)	0.017 (p = 0.87)	−0.081 (p = 0.43)
2015	0.127 (p = 0.21)	0.238 (p < 0.02)	0.270 (p < 0.01)	0.038 (p = 0.71)	−0.149 (p = 0.14)	−0.251 (p < 0.02)
2016	−0.130 (p = 0.20)	0.170 (p = 0.10)	0.131 (p = 0.20)	0.015 (p = 0.88)	−0.131 (p = 0.20)	−0.165 (p = 0.10)

These results are in agreement with those by Vasconcelos et al. (1998, 1999). Both studies used a seroepidemiological analysis on the population of Fortaleza, Brazil, and Island of São Luiz, respectively. Vasconcelos et al. (1998) observed that people living with one to five MWs presented a higher dengue incidence as compared with people living with less than one MW or more than five MWs. Such study also provided an evaluation of educational levels and dengue contamination. It showed that people with middle and high school degrees were the most contaminated (45.2% and 49.6% of the population, respectively). On the other hand, people with no degree, with primary school degree, and with college degree have 39.6%, 34.3% and 38.9% of the population contaminated respectively. Although none of these studies provided an explanation for this apparently illogical relationship, our speculation of underreporting and informal employments in lower socioeconomic levels seems to be a plausible contributing factor. Therefore, understanding the social‐economic dynamics can be the starting point for future research to figure out the optimum population class for vectors control programs to target.

In general, it is understood that low socioeconomic‐related conditions favor the transmission of dengue and other mosquito‐borne diseases (Hagenlocher et al., 2013; Moreno‐Madriñán & Turell, 2017; Ramos et al., 2008; Reiter et al., 2003). Among the common influencing pathways for this relationship are the higher indoor temperatures due to the absence of air conditioning, more in‐and‐out access for mosquitoes due to open nonscreened windows, and sociocultural practices such as water storage as a result of the lack of continuous supply of piped water (Chareonviriyaphap et al., 2003; Eisen et al., 2014). In addition, lower socioeconomic neighborhoods have greater housing density, less access to vector control, and more likelihood of abandoned containers potentialy filled with water (Eastin et al., 2014). Among other possible risk facors are the less opportunity for early diagnosis and treatment, along with more likelihood for people to be outside as a result of less enhancements within households (Moreno‐Madriñán & Turell, 2018).

Although the spatial assessment did not provide an accurate spatial model, it could highlight some trends that should be investigated in the future, especially by policy makers in São Paulo municipality. Moreover, this same concluded observation could also be valid for other diseases which have the same vector, such as yellow fever, Zika (Wong et al., 2013), and chikungunya (CDC, 2016). It is important to emphasize not only that the results presented in this research are valid as a study case for São Paulo, Brazil, but that the same methodology and assessment could also apply to other tropical geographical localities.

Conclusions

In this work, we suggest the usefulness of considering the SSTA index trendlines and/or the continuous wavelet transformation with an SSTA time series to characterize the intensity and temporal distribution of dengue outbreaks in São Paulo municipality, Brazil. Although the effect of ENSO in São Paulo is not well understood and we do not imply a robust relationship with wavelet analysis for prediction, this technique could help estimate the intensity and magnitude of dengue incidences for the next season and be a starting point to develop a more advanced prediction tool. An interseasonal model using this technicque could be developed to initiate proactive public awareness, education campaigns, and planning of budgets and resources for prevention and mitigation. Thus, the use of this technique could be an important tool for policy makers as they budget for vector control programs. However, more validation studies in different geographic locations are needed for the standardization of this method. The spatial assessment of this study was not able to detect important influence from the remotely estimated environmental parameters considered except for minimum temperature for which just a low influence was detected (r = 0.357). Nevertheless, the MWs analysis implies that socioeconomic factors should also be considered to optimize a potential forecasting tool for locating areas at a higher risk for the occurrence of dengue outbreaks. The relevance of the presented study is the contribution to potential development of new strategies for the management of vector control programs. Although the coupling of temporal and spatial analysis was not conducted simultaneously in this study because of the low availability of cloud‐free images over São Paulo municipality, this could be a powerful management tool. The use of higher temporal resolution images could improve the frequency of spatial assessment, and it could be tied to the temporal analysis. Finally, although the present study focused on dengue outbreaks, these analyses may be adopted to investigate other diseases that share the same vector.

Conflict of Interest

The authors declare no actual or perceived conflicts of interest.

Supporting Information S1

Click here for additional data file.