Remote Sensing and Multi-Criteria Evaluation for Malaria Risk Mapping to Support Indoor Residual Spraying Prioritization in the Central Highlands of Madagascar

Abstract

The National Malaria Control Program (NMCP) in Madagascar classifies Malagasy districts into two malaria situations: districts in the pre-elimination phase and districts in the control phase. Indoor residual spraying (IRS) is identified as the main intervention means to control malaria in the Central Highlands. However, it involves an important logistical mobilization and thus necessitates prioritization of interventions according to the magnitude of malaria risks. Our objectives were to map the malaria transmission risk and to develop a tool to support the Malagasy Ministry of Public Health (MoH) for selective IRS implementation. For the 2014–2016 period, different sources of remotely sensed data were used to update land cover information and substitute in situ climatic data. Spatial modeling was performed based on multi-criteria evaluation (MCE) to assess malaria risk. Models were mainly based on environment and climate. Three annual malaria risk maps were obtained for 2014, 2015, and 2016. Annual parasite incidence data were used to validate the results. In 2016, the validation of the model using a receiver operating characteristic (ROC) curve showed an accuracy of 0.736; 95% CI [0.669–0.803]. A free plugin for QGIS software was made available for NMCP decision makers to prioritize areas for IRS. An annual update of the model provides the basic information for decision making before each IRS campaign. In Madagascar and beyond, the availability of the free plugin for open-source software facilitates the transfer to the MoH and allows further application to other problems and contexts.

1. Introduction

Malaria is an infectious disease transmitted by female Anopheles mosquitoes. It remains a major global public health concern. In 2018, worldwide malaria cases were estimated at about 228 million, 93% of which were reported in the African Region [1]. The Roll Back Malaria (RBM) partnership has defined a global framework to end malaria [2].

As the fourth-highest cause of both morbidity in health centers and mortality in hospitals, malaria remains a public health issue in Madagascar. Until 2013, the confirmed number of simple cases of malaria in health facilities annually represented 6.5% of all outpatient consultations. This proportion was 5.6% in 2016 [3]. In Madagascar, the prevalence varied greatly according to the epidemiological facies of malaria: from 1% in the Central Highlands (CHL), to 5% in the sub-desert facies, and to 9% in the Tropical and Equatorial facies, where transmission is high and perennial [4]. Unlike eastern and western coastal areas of Madagascar, the CHL are characterized by seasonal and unstable malaria transmission. Historically, malaria in the CHL of Madagascar has been marked by deadly epidemics. During the last major epidemic in the CHL, in 1988, malaria caused an estimated 25,000 deaths [5]. Human in-country mobility is increasing between high malaria transmission and low transmission areas. A sentinel network of fever surveillance shows the persistence of malaria transmission in stable areas [6].

In high malaria prevalence areas, the main malaria control measures implemented by the Malagasy National Malaria Control Program (NMCP) are the use of insecticide-treated nets (ITN) and ProActive Community Treatment (Pro-ACT). In areas of low malaria prevalence, as in the CHL, control differs and relies essentially on indoor residual spraying (IRS) [7,8]. Global funds and other institutions, such as the President Malaria Initiative (PMI)/USAID help to fight malaria in high burden countries [9]. The NMCP benefits from those financial supports [10]. In 2015, the Malagasy Ministry of Public Health (MoH) declared that six out of 31 districts in the CHL are in the pre-elimination phase. These districts are those with a rapid diagnostic test (RDT) positivity rate under 5% in the general population [4]. It becomes crucial to consolidate these achievements and to extend new areas towards the pre-elimination phase.

For Madagascar, the implementation of IRS implies a heavy logistical mobilization and an adapted financial support mostly depending on external funding. To control the malaria epidemic in the central highlands in the late 1980s, five years of full insecticide campaign coverage were implemented. Annual programs of indoor residual spraying of DDT were carried out between December 1993 and January 1998 in most rural areas at altitudes of 1000–1500 m. This strategy helped to stop this deadly epidemic [11,12]. Since 1999, rotational and selective interventions were carried out [12,13] based on the incidence of malaria observed during the year preceding the transmission season. Some IRSs have been carried out in response to epidemics. For IRS campaigns in the CHL, the NMCP advocates an alternation of two years of total coverage and two years of selective IRS. This strategy aims to reduce IRS coverage and instead enhance selective interventions. In a context of limited resources, optimizing malaria control targets becomes essential for countries like Madagascar. Given that allocated resources will not cover costs for total coverage of IRS, the aim is to prioritize areas with the highest malaria transmission risk.

Several methods have been used to map malaria risk. GIS-based spatial modeling techniques were adopted by Rincón-Romero and Londoño for malaria risk mapping in Columbia and by Ferrao et al. in Mozambique [14,15]. Others promote spatial statistical approaches to map malaria risk. For example, Kleinschmidt et al. combined a logistic regression model and geostatistical analysis to produce a malaria risk map in Mali [16]. Yankson et al. used model-based geostatistical methods to map malaria risk in children under five in Ghana [17]. Hanaf-Bojd et al., Wondim et al., and Ali and Ahmad used the multi-criteria evaluation (MCE) method to map malaria risk in the South of Iran, Ethiopia, and West Bengal, respectively [18,19,20].

As for most, if not all, vectors, environmental and climatic factors play a critical role in the Anopheles life cycle and in the maintenance of malaria transmission. The environmental temperature impacts the internal temperature of mosquitoes [21]. The environmental conditions also considerably impact breeding sites. Wetlands, for example, are crucial for the development of the larvae of Anopheles [22,23]. The changes of environment and climate can also induce changes in vector behavior [24,25].

To overcome the lack of entomological data on malaria vector distribution at the scale of the country, a suitable habitat for Anopheles can be characterized by using an environmental proxy [26,27,28,29]. The malaria patterns depend on location and availability of breeding sites that can be recovered using spatial information [30,31].

Several risk factors of malaria have been identified in the literature, mainly through a statistical modeling approach [32,33]. Environment and climate data, such as distance from water bodies, temperature, elevation, drainage density, rainfall, and land use/land cover, are widely used because they are closely linked to malaria and its vectors. Other methods, such as the use of MCE, can be developed when epidemiological and/or entomological data are too scarce to develop statistical models [34,35,36].

MCE is a common method for assessing and aggregating many criteria [34,37,38,39]. In order to form a single index of evaluation and to provide necessary information for decision makers, MCE basically consists of combining information from several criteria. It is considered a semi-quantitative method using a participative approach where stakeholders bring their knowledge and expertise [40].

The combination of the geographical information system (GIS) approach and MCE began in the 1990s, with full integration by Eastman [41]. This combination is widely used as a decision support tool [34,42,43]. In Madagascar, where field malaria data is limited, Rakotomanana et al. used an MCE approach to evaluate malaria risk in six zones of the CHL [22]. Their study is based on environmental factors such as temperature, elevation, and rice fields as model input data. Landsat 7, Spot 4, and radar images were used to map rice fields. Yet, inhabited areas, which represent the vulnerable population, and precipitation were not taken into consideration. Furthermore, as the method can only be carried out by GIS specialists, the map has not been updated since 2007.

The main objective of this study is to map malaria risk in all CHL communes, taking into account the vulnerable population, and then validate the results by comparing them to epidemiological data. A second objective is to develop a “user-friendly” tool, enabling regular updating of risk maps, to support Malagasy MoH for selective IRS implementation.

2. Materials and Methods

2.1. Study Area

The CHL are located in the central plateau of Madagascar. They encompass 97,000 km² and include 31 districts (Figure 1), representing 27% of the country’s districts. The west side presents a smooth slope, whereas the east presents an abrupt slope towards the coastal regions. The altitude of the CHL ranges from 200 to 2700 m. The average annual temperature is 20 °C and the climate is subdivided into two seasons: rainy season (October to May) and dry season (June to September). In the CHL, the average annual precipitation is between 800 and 1500 mm. The CHL have the highest human density in Madagascar. Almost half of Malagasy people live in the CHL. In 2014, the population in the CHL was estimated at 9,200,000 inhabitants, with an annual growth rate of 2.8%. More than 80% of Malagasy people live in rural areas. The main activities of the population are agriculture and livestock farming. Rice is the main food and rice fields are highly abundant in Madagascar, especially in the CHL.

2.2. Land Cover Update

A total of 18 images from the Landsat 8 sensor, with 30 m spatial resolution, were used to update the land cover map from 2014 to 2016. Images were downloaded from the United States Geological Survey (USGS) earth explorer website (https://earthexplorer.usgs.gov/). The object-based image analysis (OBIA) classification approach was used. Unlike the standard per pixel approach, the OBIA approach does not treat the pixels one by one, but in context by grouping pixels within objects based on their spectral value, size, shape, and context [44]. Firstly, an object segmentation of images was performed in order to group all pixels having the same radiometric characteristics. Then, the segmented images were classified into nine classes through a membership rule-based method where each segment attribute had conditions to respect. Specifically, each land cover class was classified according to pre-defined rules. The implementation of these rules was based on thresholds (

Table 1. Classification parameters for land cover.

Classes	Parameters
Rice field	GLCM contrast ≤ 70; LWM ≤ 75; Mean layer ≤ 85
Water body	NDWI > 0.37; Sum of reflectance (Σ (b2, b3, b3)) < 204
Hydrographic network	Brightness > 65; LWM < 56; NDVI < 0
Wet cultivation	−0.15 < NDVI < −0.01
Other	Everything that is not assigned to the above four classes

Note: GLCM (Gray Level Co-Occurrence Matrix), LWM (Land and Water Masks) = b5/(b3 + 0.0001) × 100, NDVI (Normalized Difference Vegetation Index), b2: blue band, b3: green band, b4: red band, b5 = near infrared band (Spectral band).

).

As Anopheles breed in various areas, four land cover classes were identified as areas for mosquito development that were considered potential breeding sites: rice field, water body, hydrographic network, wet cultivation [45].

A post-classification phase was performed to assess the accuracy of the classification results. The ground truth data were obtained from field missions in five districts of the CHL during rainy season. Three confusion matrices were generated in order to represent the misclassification, as well as classification accuracy for 2014, 2015, and 2016. Cohen’s Kappa index was calculated for each confusion matrix to test for the concordance of the classification.

2.3. Spatial Multi Criteria Evaluation

Three yearly malaria risk models were established following the workflow shown in Figure 2 using six criteria involved in malaria risk.

2.3.1. Criteria Identification and Differentiation

The choice of the criteria was based on a literature review and seven expert opinions (epidemiologists, entomologist, a doctor, geomaticians, and a modeler).

Inhabited zone: Malaria only occurs in inhabited areas as humans are the only reservoir of malaria. Village locations, in vector format (point type), were extracted from the database of the Malagasy Cartographic Institute named “Foiben-Taosarintanin’i Madagasikara” (FTM) in 2000. In order to obtain updated village locations, these were combined with data gathered from OpenStreetMap (https://www.openstreetmap.org/) and data obtained from 2013 to 2016 field missions. A 1 km buffer, representing inhabited areas, was created from the previously combined data and then rasterized to a 30 m spatial resolution.

Population density: The last national population census campaigns took place in 1993 and 2018, but data from the latter is not yet accessible. Given the lack of updated data, the population density was extracted from the WorldPop project database (www.worldpop.org.uk). WorldPop is an open-access library providing population data with a spatial resolution of 100 m. The 2010 population density was used for the 2014 model, while data of 2015 were used for 2015 and 2016 models. Since data with better resolution do not exist at the present time, a resampling to 30 m was performed.

Elevation: Elevation is important in malaria risk assessment because it is directly related to temperature, and impacts the life cycle of the mosquito vectors [46,47,48] and the development of Plasmodium species, which are the parasites responsible for malaria disease. Elevation data were obtained from digital terrain models provided by the Shuttle Radar Topography Mission (SRTM) campaign held in 2000. The SRTM data for CHL, with a spatial resolution of 30 m, was downloaded from the USGS Earth explorer website (https://earthexplorer.usgs.gov/). Two tile sets of SRTM images (SRTM_46_16 and SRTM_46_17) were used to cover all CHL.

Temperature: Temperature influences mosquito development and longevity [49,50]. The hot and rainy season lasts from October to May in Madagascar. During this period, the population densities of Anopheles and malaria cases are high. Mean temperatures for this timeframe, from 2014 to 2016, were computed. The USGS, through the International Research Institute (IRI) library (http://iridl.ldeo.columbia.edu/), provides temperature data every eight days with a spatial resolution of 1 km. This source platform allowed extracting a 30 m spatial resolution through its expert mode tool.

Precipitation: Rainfall influences the availability of mosquito breeding sites; it may predict vector abundance [51,52]. Rain increases the chance of larval habitat availability; however, excessive rainfall can cause leaching of larval breeding sites. Daily precipitation with a spatial resolution of 12 km data was obtained from the IRI website (http://iridl.ldeo.columbia.edu/) provided by the United States National Oceanic and Atmospheric Administration (US NOAA) agency. These were GeoTIFF raster images, in which values are expressed in millimeters of rain. Using the same expert mode feature as for temperature, a resolution of 30 m was extracted. Total precipitation during the October to May period was calculated for each studied year by adding daily values during the eight months.

Distance to wetland: Extracted from the updated land cover for 2014 to 2016, at a resolution of 30 m, distance to wetland was used as an environmental determinant of malaria transmission risk. Agricultural parcels, standing water, and streams are considered as potential mosquito breeding sites [53]. A maximum distance of 5 km from wetlands was formed to take into consideration mosquito blood meal seeking and resting areas [54]. It was considered that the malaria risk decreased progressively as distance increased within the buffer.

2.3.2. Factor Standardization

Criteria were categorized and divided into constraints and factors. Constraints are defined as masks that consider whether a zone will be part of the calculation [22]. The process uses a Boolean function by classifying the criteria into suitable or not suitable areas (

Table 2. Constraint descriptions and parameters.

Criteria	Description
Inhabited zone	BF: recoded to 1 for inhabited areas that are potentially at risk, and to 0 for uninhabited area that are not at risk
Elevation	BF: recoded to 0 for elevation <1000 m (permanent risk) and >1500 m (no risk); recoded to 1 for elevation between 1000 m and 1500 m
Population density	BF: recoded to 1 in areas with d <800 pop/km² and to 0 in areas with d ≥800 pop/km²

Note: BF (Boolean function), d (Population Density).

). Eastman described factors as criteria that define a certain degree of aptitude for all regions [55]. In this study, factors are expressed by continuous values and act in a progressive way on the aptitude following fuzzy logic functions (

Table 3. Threshold for the fuzzy logic function used in the spatial MCE models.

Criteria	Function	Control Points
		a	b	c	d
Population density	Decreasing sigmoid	-	-	400 ≤ d < 800	d ≥ 800
Elevation	Increasing sigmoid	500	1500	-	-
Distance to wetland	Decreasing sigmoid	-	-	1000	5000
Precipitation	Symmetric sigmoid	0	80	1000	2000
Temperature	Symmetric sigmoid	18	28	32	35

Note: d is in population per km². Elevation and distance to wetland are in meters (m). Precipitation is in millimeters (mm). Temperature is in degrees Celsius (°C).

).

Each factor was integrated into a GIS and mapped. Factors were quantitative variables (e.g., temperature, elevation) and measured in various units (e.g., degrees Celsius, meters). Standardization was performed to bring the factors on a continuous scale of aptitude, ranging from least aptitude (0) to most aptitude (255). Standardization enables the comparison and combination of factors for each pixel in the study area. The fuzzy logic functions were used to represent the transition between the totally unsuitable and the fully suitable to better represent real life settings [56]. Here, factors were defined using sigmoid membership functions (increasing, decreasing, or symmetrical). Functions and parameters are presented in Table 3.

2.3.3. Factor Weighting

The second step of the MCE consisted in weighting factors, i.e., assigning weights to each standardized factor according to their degree of relevance. To determine the weight of each factor, an analytical hierarchy process (AHP) was implemented [57]. This method consists of a pair-wise comparison of the factors [22,57]. Weights assigned to each factor are determined by experts’ contribution using Saaty’s continuous rating scale [57]. Experts were asked to compare factors two by two and to determine their importance. As per this scale, factors were rated as: 9 for extremely more important, 7 for very highly more important, 5 for highly more important, 3 for moderately more important or 1 for equally important than the second factor. Inversely, less important factors were assigned a score of: 1/3 for moderately less important, 1/5 for much less important, 1/7 for very much less important or 1/9 for extremely less important. This procedure generates a table of the pair-wise comparison of all factors. In this table, the sum of the weights of all factors equals 1. The pair-wise comparison used in this study is detailed in Supplementary File (Table S1). The consistency of the comparison is then assessed by computing the consistency ratio (CR). This corresponds to the probability that ratings were randomly generated. The consistency of the rating is considered acceptable when CR <0.1 [58].

2.3.4. Criteria Aggregation

A weighted linear combination (WLC) aggregation of criteria was the last step of the MCE. It consisted of multiplying the weights from the weighting of each standardized factor and then summing them. To exclude unsuitable areas, i.e., areas where there was no risk of malaria, this sum was multiplied with Boolean constraints. The aggregation of criteria per Equation (1) resulted in a gradient risk map of malaria, per pixel of 30 m.

$$\mathrm{S}={{\displaystyle \sum }}_{i=1}^n{w}_i{X}_i{{\displaystyle \prod }}_{j=1}^m{C}_j,$$

(1)

where $\mathrm{S}$ is the aptitude for an event, $W_i$ is the weight of factor $i$, $X_i$ is the factor $i$, $C_j$ it the constraint j, n is the number of factors, and m is the number of constraints.

Malaria risk models were constructed for the years 2014, 2015, and 2016.

As the NMCP uses communes as the scale for IRS intervention in Madagascar, risk gradient per pixel was therefore aggregated by commune.

Aggregation of pixel risk values at commune level was carried out following three steps. Firstly, to determine the probability distribution of gradients per pixel, using the maximum likelihood method, six candidate distributions (normal, log-normal, geometric, gamma, Poisson, and Weibull) were tested in 100 randomly selected communes. Secondly, outlier values were detected according to the respective distribution selected per commune and corrected with the confidence interval of the crude mean to derive an adjusted mean. Finally, the optimal class number to perform the Jenks natural breaks algorithm of commune-level risk gradients was determined through the Gaussian mixture model [59,60].

2.3.5. Validation of the Model Framework

Model output was compared to data on Annual Parasite Incidence (API) provided by the NMCP. These data record the incidence rate of malaria per commune and per year. API validation data was only available for 208 communes and for 2015.

The area under the curve (AUC) of the receiver operating characteristic (ROC) curve was used to assess the model performance. The perfect model shows a value of 1.0 in the ROC curve [61]. We used the 2015 API data to validate the output of the model for the 2016 IRS campaign. The 2015 API in each commune, corresponding to observed values, was recoded as a binary variable to compute the AUC ROC: API <1‰ was considered as low incidence and API ≥1‰ as high incidence. Sensitivity and specificity were computed using the best cut-off of the classified malaria risk map (Figure 4).

2.4. Uncertainty Analyses

Uncertainty analyses were performed to verify the stability of the model. Indeed, the parameters used in the development of the model mainly depended on the choice of experts, as well as on the literature. These analyses consisted of varying the weights of the various standardized factors of the model, and checking the impact of these variations on the standard deviation of model output. Ten variations of the initial weights corresponding to 10 different expert choices (test at ±5%, ±10%, ±15%, ±20%, and ±25%) were simulated. The adjusted weights were calculated using Equation (2):

$${\mathrm{w}}_{\mathrm{i}}=\left(1-w_m\right)\ast \frac{w_{i0}}{1-w_{m0}}, 1\le \mathrm{i}\le \mathrm{n}, \mathrm{i}\ne \mathrm{m},$$

(2)

where ${\mathrm{w}}_{\mathrm{i}0}$ and ${\mathrm{w}}_{\mathrm{m}0}$ are, respectively, the weights of the ith risk factor and, in the base model, of the main changing risk [62,63,64].

2.5. Plugin Development

In order to provide a dynamic, simple, and free decision support tool for the NMCP, an extension of the QGIS 2.x software, called “MCE for Public Health”, was developed. The Python programming language was used because it is the most suitable language for QGIS. To make the tool more efficient, all calculations were performed with the application programming interface and Geospatial Data Abstraction Library (GDAL) of QGIS. Qt Designer Software was used to design the graphical user interface (GUI) [65]. It was designed to be ergonomic, easy to use, and to prevent user typing errors [65].

2.6. Change Detection

In order to evaluate the density of change in the different zone surfaces (with a 30 m step), in terms of area at risk, a change detection process was carried out. It was applied to the 2014, 2015, and 2016 spatial models. According to Singh, change detection can be defined as a process identifying state changes of an object and/or phenomenon on different dates [66]. For this study, image differencing was adopted. It consists of evaluating the difference, pixel by pixel, between two images at different dates. This method is thoroughly explained by Singh [66]. The risk gradient, per class, was calculated in order to appreciate the tendency of each risk class from one year to another. Through this method, zones subjected to changes had different statistics than those that had not changed. A risk gradient with a declining surface shows negative pixel values, while gradients that have an increased surface have positive values.

3. Results

3.1. Land Cover Updating

The application of the OBIA classification method, through the membership rules, allowed classifying the land cover in the CHL. The classification showed that the rice fields and wet cultivation classes represented 17% in 2014, 27% in 2015, and 11% in 2016. The water body and hydrographic networks classes covered 6%, 6%, and 5%, respectively, in 2014, 2015, and 2016. The remainder is covered by the “other” class, which includes bare soils, built-up areas, and clouds.

Based on the strength of agreement established by Landis and Kock, the results of the land cover classifications for 2014, 2015, and 2016 were statistically satisfactory (

Table 4. Cohen’s Kappa index for 2014, 2015, and 2016 land cover classifications.

	Kappa Index
Land cover classification (2014)	0.7741
Land cover classification (2015)	0.7438
Land cover classification (2016)	0.8320

) [67].

3.2. Malaria Risk Model

According to their respective experiences, experts working on the malaria program parameters in relation to the transmission of malaria attributed the weights of each standardized factor (

Table 5. Factor weights.

Factors	Weights
Population density	0.4990
Distance to wetland	0.1824
Temperature	0.1698
Elevation	0.0910
Precipitation	0.0577

). The CR for the pair-wise comparison matrix (Table S1) showed a value of 0.06, which indicates that the conducted comparison was consistent. Population density was considered to be the most important factor in malaria risk, followed in decreasing order of importance by distance from wetlands, temperature, the altitude, and precipitation.

The MCE, through the application of the WLC, allowed having a malaria risk map with a gradient ranging from 0 to 255 for each pixel of 30 m spatial resolution. Areas with a value of 0 had the lowest risk, while areas with a suitability of 255 had the highest risk. Figure 3 presents the different maps of malaria transmission risk for 2014, 2015, and 2016, with a 30 m pixel scale. Figure 4 shows the 2016 spatial model map, adjusted to the commune scale. Three distributions were found to best fit the commune-level risk (Weibull, log-normal, and normal). According to the adjusted mean and available malaria incidence data for validation, three optimal classes were selected to categorize the risk values. CHL presented 60%, 32%, and 8% of communes with high (rating of 3), moderate (rating of 2), and low malaria (rating of 1) risk, respectively. The central and eastern parts of the CHL were dominated by communes at a high risk of malaria. Communes with low risk were located on the western side of the CHL and the other communes had a moderate malaria risk.

Giving the lack of API validation data, only the malaria risk model of 2016 was validated. Validation showed an acceptable fitting with an AUC of 0.736 (95% confidence intervals (CI): [0.669–0.803]) (Figure 5) [68].

Table 6. Results of ROC analysis.

	2015 Annual Parasite Incidence
2016 Model Output	High	Low
High	TP: 111	FP: 49
Low	FN: 11	TN: 37
Accuracy: 0.712 Sensitivity: 0.910 Specificity: 0.430

Note: TP (True Positive), TN (True Negative), FP (False Positive), FN (False Negative), Accuracy = (TN + TP)/(TN + TP + FN + FP), Sensitivity = TP/(TP + FN), Specificity = TN/(TN + FP).

shows the results of ROC analysis. Using a rating of 2 as a threshold, a sensitivity of 0.91 (11 positive cases were missed) and a specificity of 0.43 (49 negative cases labeled as high risk) were obtained. By applying this cut-off, 91% of the communes were identified as high risk by the model.

3.3. Uncertainty Analysis

The CHL uncertainty map (Figure 6) showed that the maximum value of standard deviation was significantly below 0.1, indicating a robust model [63,64]. Calculated standard deviations revealed values around the mean, demonstrating that the risk model of malaria remains stable even when the weights assigned to each standardized factor are changed by the stakeholders.

3.4. “MCE for Public Health” Plugin

The plugin is an interactive tool, adapted for the QGIS 2.x version: a dynamic, free, and open source semi-automatic tool. It groups the main steps for MCE evaluation: reclassification of constraints, normalization, computing weight for each factor using a pair-wise comparison matrix, and aggregation of factors. The plugin runs for a maximum of 10 constraints and 15 factors. It can support four raster formats of input and output data: rst (idrisi format), gtiff (GeoTIFF format), img (erdas imagine format), and jpeg. The Full MCE plugin file is available on GitHub (https://github.com/SaGEOTeam/FullMCE). The plugin is presented both in French and English. Figure 7 shows the main interface of the developed plugin.

3.5. Change Detection on Models

Two change detections of malaria risk were carried out to determine if the risk category had changed. The results of the image differencing technique, for 2015–2014 and 2016–2015, were two images with pixel coded in three categories: increase in risk, no change, and decrease in risk (Figure 8).

Based on the three risk models, at pixel scale,

Table 7. Surface affected by changes resulting from the image differencing method.

	Surface 2015–2014 (km²)	Surface 2016–2015 (km²)
decrease	12,198.96	21,435.32
no change	12,334.38	3883.71
increase	10,086.82	9301.13
Total	34,620.16	34,620.16

shows the changes on the risk surface of 2015–2014 and those of 2016–2015. Between 2014 and 2015, 64% of the surfaces in the study area were subject to risk change where 35% affected the transition to a higher risk than in 2014. Between 2015 and 2016, 89% of the surface changed risk with 61% surfaces subject to a decrease in risk. The models have not changed significantly between 2014 and 2015. From 2015 to 2016, the drastic decrease in rainfall and wetland during the malaria transmission season led to a decrease of the high-risk class of the model and a subsequent increase of the moderate class.

4. Discussion

The MCE method had previously been carried out in six areas of the Malagasy CHL [22]. Here, the model is improved by including all areas of the CHL and by taking into account rainfall and population distribution. The comparison of available API data with the risk model enabled an assessment of the accuracy of MCE in this study. Validation of the malaria risk model showed a good result. For the first time, to our knowledge, the development of the “Full MCE” plugin enables non-specialists, such as staff of the Malagasy NMCP, to carry out a complete MCE approach with a free spatial decision tool. This initiative allowed the use of a single software platform to perform all the MCE treatments. In April 2019, NMCP and other stakeholders (Ministry of Public Health, Ministry of Agriculture, universities and research institutes) from 10 countries were trained to use the plugin and the tool was distributed to them. An update with the latest version 3.x of QGIS is currently in progress to allow a wider diffusion.

The availability of Landsat 8 images allowed updating the Malagasy CHL land cover map with a medium spatial resolution of 30 m. The previous land cover map had dated from 2009. Although the OBIA approach is better suited for images at a high spatial resolution, it is also currently used for images with medium spatial resolution such as Landsat images [69,70]. At the scale of the present study (i.e., commune scale), the results obtained from the Landsat 8 images seem appropriate.

The MCE approach is a knowledge-based and pragmatic method adapted to situations in which there is a lack of data. This method is rapid and easy to implement [36]. It is an alternative to a statistical method requiring large datasets and is adaptable and/or transposable to different areas [36] with particular aspects of malaria transmission. AHP is a method which is very commonly used [22,63,64,71] to weigh factors in MCE. It greatly facilitates the work assigned to experts as (i) pair-wise comparisons are easier to carry out for experts than simultaneous comparison of all factors, and (ii) it is simpler to compare factors using qualitative ordinal variables (“extremely more important than”, etc.) rather than defining precise weights. AHP also has the advantage of enabling an easy assessment of the ratings consistency through the CR.

Most studies on spatial modeling of malaria transmission risk [22,72] do not take into account the human reservoir aspect in the models. Previously, only population density was used, while this variable supports all pixels across the entire study area, even uninhabited areas. The notion of reservoir is however important in terms of risk of disease and particularly in public health [73,74]. The particularity of the present model is that it takes into account only the areas where the population is present, which represent a potential human reservoir.

GIS-based decision support was recently used to target IRS control programs in Zambia [75,76]. Chanda et al. used epidemiological and entomological data to determine areas that needed interventions [75]. In Madagascar, the geographical extent coupled with the issue of being landlocked make routine surveillance and data collection system difficult. Despite efforts from the NMCP stakeholders, complete and up-to-date data are usually not available before elaborating the annual IRS target. The present model shows that MCE approaches overcome these difficulties.

However, the present study has some limitations that need to be mentioned. Data collection on climate and demography was difficult in Madagascar, whether free or paid. Indeed, the number of meteorological stations is not sufficient for national use, and their interpolation may have a direct impact on the result. The scarcity of meteorological stations imposed us to use remotely sensed substitute data. The IRI for Climate and Society provided a solution to this issue by integrating climate and environmental data from various sources [77]. Forecast climate data was previously used to predict malaria risks in the African region [28,77,78]. Dambach et al. confirmed the role of land surface temperature and rainfall in mosquito density [79].

Datasets used to compute malaria risk were not always available at similar scales. Population data was only available with a spatial scale of 100 m. Lack of access to high spatial resolution data is frequent in low income countries and considerably hinders the quality of produced risk maps. Furthermore, using datasets with different spatial resolutions is clearly not ideal, and it is recommended to downgrade the resolution of all datasets to the lowest spatial resolution. Here, we took the risk to resample datasets to 30 m but finally produced an operational risk map at the commune level (i.e., with a much coarser resolution than the one used for the datasets). The latest general population censuses in Madagascar were conducted in 1993 and in 2018 but results of the latter are not yet available. When the 2018 Malagasy population census data is released, updated population density maps with finer resolution will be produced and should considerably improve incidence and prevalence estimates, as well as risk maps produced in Madagascar (including the malaria maps derived from this work).

The available API data, required for model validation, did not cover the entire study area. Only about half of the study area had complete data covering the period of malaria transmission. This lack of API data could significantly impact the accuracy of the model validation.

In our study, malaria control data, such as recent IRS campaigns, were not taken into account. It can be hypothesized that they impacted malaria transmission data (API) and thus model validation. Since these intervention data influence malaria suitability, they can easily be considered in future malaria risk model constructions.

Despite these limitations, the tool synthesizes current knowledge on malaria transmission in the Malagasy CHL. Based on the available data, it provided information that was useful to decision makers. An improved version of the map could be made if better data and knowledge were available.

From the perspective of spatial modeling of malaria risk in other regions, particularly in coastal areas, it would be interesting to produce a risk map at the national level. Madagascar is currently divided into five malaria transmission strata (CHL, fringes, East, West, and desert south) [80,81]. These strata represent endemic areas (coasts) and low transmission areas (CHL and fringes). These areas differ according to their transmission intensity, ecozone, and vectors [80]. The development of a model for the other epidemiological facies might require using different variables or different weights. These variables and parameters need to be determined by experts.

This study focuses on mapping malaria risk in CHL, an area with low transmission. Population movement flow is an important factor impacting malaria risk because malaria risk increases in areas receiving people from high transmission zones [81]. In 2017, Girond et al. developed an early warning system. It used new technologies for early detection and forecasting of malaria based on a Malagasy sentinel surveillance system [82]. This early warning system combined with our spatial model would allow a more accurate mapping of malaria risk and better predictions of epidemics.

5. Conclusions

Malaria transmission is limited by environmental and climatic criteria because they directly affect the life cycle of Anopheles and the Plasmodium they transmit. GIS combined with MCE, through their capacity for storage, data management, analysis, modeling, and mapping of spatially referenced data, is a useful tool to apprehend spatial decision issues [83]. This combination has improved the understanding of areas at risk of malaria. The resulting risk map is used for decision making to target priority communes to focus IRS campaigns in Madagascar. The “Full MCE for Public Health” tool, which is dynamic, fast, and easy to use, should be easily appropriated by decision makers to prioritize IRS, and its flexibility makes it easy to use in other contexts: for other diseases and for other countries, to simulate various scenarios. It will be useful to a much wider community of stakeholders involved in risk assessment, especially in areas where data is lacking.