Week 4: Multi-Decision Criteria Analysis

1 Suitability Mapping: Part 1

1.1 Introduction

Now, with the proliferation of open spatial dataset, risk models derived from environmental, sociodemographic and topological factors are becoming increasingly available for open research. Such models have broadly shown to be useful in delineating geographical areas of risk or suitability for a certain outcomes. Today, we are going to explore this qualitatively using a knowledge-driven approach.

1.1.1 Learning outcomes

To provide an introductory overview to the applicability of knowledge-driven methods, in particular, we are going to learn the Multi-Criteria Decision Approach (MCDA) which is a method that use decision rules from existing knowledge to identify areas potential suitability for an outcome. It is especially useful in data-sparse situations, or when for the first time exploring the potential geographical limits of certain outcome.

For instance, using modest number of raster layers such as population density, urbanisation, approximating to street segments, house prices and deprivation; it is possible to combine such information so to determine regions for which crime events such as burglaries are likely to occur, or suitable in that matter. This approach has been widely used in a number of disciplines over the last decades, and has gained prominence in public health related fields such as vector-borne disease prevention, and disaster sciences such as landslides. We will learn how to apply these methods to the two context.

1.1.2 Datasets & setting up the work directory

Before you begin do make sure to download all data by clicking here. Create a folder on called “Week 4” within your “GEOG0114” folder stored in the desktop of your personal computer. Make sure to extract all data from the zip folder and store it into “Week 4” folder. Open a new R script and set the work directory to Week 4’s folder.

For Windows, the work directory will be:

setwd("C:/Users/AccountName/Desktop/GEOG0114/Week 4")

For MAC, the work directory will be:

setwd("/Users/AccountName/Desktop/GEOG0114/Week 4")

1.1.3 Loading and installing packages

We will need to load the following packages:

• sf: Simple Features
• tmap: Thematic Mapping
• sp: Package for providing classes for spatial data (points, lines, polygons and grids)

The above packages sf, tmap & sp should have been installed in the previous session(s). We will need to install a new package:

• raster: Raster/gridded data analysis and manipulation.
• spatialEco: Spatial Analysis and Modelling Utilities package, which provides the user further functions for dealing with raster data.

install.packages("spatialEco")
install.packages("raster")

# Load the packages with library()
library("sf")
library("tmap")
library("raster")
library("sp")
library("spatialEco")

1.1.4 Loading datasets

We will be dealing with raster data only for this exercise. These will be a series of climate, environmental and other social-anthropogenic gridded data from Kenya. We are going to combine them in order to determine areas that are suitable for disease transmission of Lymphatic Filariasis, a neglected tropical disease, caused by a range of mosquitoes that spread the parasitic worms called the Wuchereria bancrofti through a mosquito bite.

There is evidence that indicates that precipitation, temperature, population density, aridity, dryness, land surface elevation and vegetation are risk factors for mosquito infestation, and in turn, greatly influencing the way mosquitoes breed and feeding behaviour within human dwellings thus maintaining disease transmission. We will use the following techniques for mapping areas of suitability:

• Binary Classification & Ranking
• Suitability Mapping using Weighted Overlay Analysis
• MCDA Mapping

Let us begin loading the following list of raster files, each is a variable of interest:

• Raster: Mean Temperature in Kenya (1000m resolution) named Kenya Mean Teperature.tif
• Raster: Normalized Difference Vegetation Index in Kenya (1000m resolution) named Kenya NDVI.tif
• Raster: Precipitation in Kenya (1000m resolution) named Kenya Precipitation.tif
• Raster: Population Density in Kenya (1000m resolution) named Kenya Population Density.tif
• Raster: Land Surface Elevation in Kenya (1000m resolution) named Kenya Elevation.tif
• Raster: Levels of Dryness in Kenya (1000m resolution) named Kenya Aridity Index.tif

temp <- raster("Kenya Mean Teperature.tif")
nvdi <- raster("Kenya NDVI.tif")
prec <- raster("Kenya Precipitation.tif")
popl <- raster("Kenya Population Density.tif")
elev <- raster("Kenya Elevation.tif")
arid <- raster("Kenya Aridity Index.tif")

Load the shapefile nation and state borders for Kenya:

• Shape file: Kenya’s National border named Kenya_Border_3857.shp
• Shape file: Kenya’s State border named Kenya_States_3857.shp

# load the shapefiles
kenya_border <- st_read("Kenya_Border_3857.shp")
kenya_states <- st_read("Kenya_States_3857.shp")

Note

All shape file and raster data were in projected to the CRS: Spherical mercator 3857

You can inspect each raster to know its dimension, extent, resolution and minimum and maximum values. Since are going to stack all the rasters together, you definitely want the: dimension, extent and resolution to be the same. If there’s a slight difference the stack won’t work.

# for instance temp and aridity
temp
arid

To visualize raster data - you can write the following code:

# For instance take the mean temperature for Kenya
# Spectral colours are useful for diverging scales "-brewer.spectral" is Rd-Or-Yl-Gr-Bu. "-brewer.spectral" reverses the order
tm_shape(temp) + 
    tm_raster(col.scale = tm_scale_continuous(values = "-brewer.spectral"),
        col.legend = tm_legend(title = "Mean Temperature", frame = FALSE)) +
tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
tm_layout(frame = FALSE)

1.2 Suitability mapping using Simple Overlays

1.2.1 Setting up binary classification on raster

One thing to know about this approach - the mappings are purely knowledge-driven and not analytical. For instance, according to previous studies (see example of papers: 1 & 2), we have learnt a bit about the limits or thresholds for these predictors from or below which we can consider an area being suitable for the occurrence of Lymphatic Filariasis (LF).

We have summarized these limits or thresholds in a table as follows:

Risk factor	Threshold for suitability
Temperature	> 15.0 (degree Celsius)
Precipitation	> 350 (mm of rainfall)
Aridity	> 0.20 (i.e., semi-arid (> 0.5) dry-sub-humid)
Elevation	< 1200m (low land)
Population Density	> 0 (Inhabitants must exist source for blood meal)
NDVI	> 0.5 (mild levels of vegetation)

We should use the aforementioned thresholds to produce binary or Boolean maps using the above criteria. This means that the pixel values of each raster layer will be equal to 0 indicating that its an unsuitable condition for LF transmission, and 1 for suitable conditions for LF transmission.

Let us reclassify each layer according to the above criteria, starting with temperature:

# reclassify temperature as a binary or Boolean layer
temp
# lowest value = 1.2
# highest value = 29.6
# reclassify anything below 15 as 0, and the rest above 15 as 1

# reclassify the values into two groups 
# all values <= 15 change to 0
# all values > 15 change to 1
temp_cl <- c(
    -Inf, 15, 0, 
    15, Inf, 1
    )

# convert into a matrix format
temp_cl_mat <- matrix(temp_cl, ncol = 3, byrow = TRUE)
# see matrix
temp_cl_mat
# apply matrix to reclassify() function to categorize the raster accordingly
temp_recl <- reclassify(temp, temp_cl_mat)

When you reclassify the raster for temp to temp_recl. This is what the output should look like:

tm_shape(temp_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("grey", "#F1948A"), labels = c("Unsuitable (<15.0)", "Suitable (15 & above)")),
        col.legend = tm_legend(title = "Temperature", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

Lets repeat the same steps for elevation to reclassify it’s layer according to the given criteria:

# reclassify elevation as a binary or Boolean layer
elev
# lowest value = -11m (below sea level)
# highest value = 4865m (above sea level)
# reclassify anything below 1200m as 1, and the rest above 1200m as 0

# all values <= 1200.0 change to 1
# all values > 1200.0 change to 0
elev_cl <- c(
    -Inf, 1200, 1, 
    1200, Inf, 0
    )

# convert into a matrix format
elev_cl_mat <- matrix(elev_cl, ncol = 3, byrow = TRUE) 
# see matrix
elev_cl_mat
# apply matrix to reclassify() function to categorize the raster accordingly
elev_recl <- reclassify(elev, elev_cl_mat)

The elevation output should look something like:

tm_shape(elev_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("grey", "orange"), labels = c("Unsuitable (>1200m)", "Suitable (1200m & below)")),
        col.legend = tm_legend(title = "Elevation (m)", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

Try reclassifying the remaining raster arid, prec, popl and ndvi. The solutions are provided in the code chunks below.

Click here

Solutions

# nvdi
nvdi

nvdi_cl <- c(
    -Inf, 0.5, 0, 
    0.5, Inf, 1
    )

nvdi_cl_mat <- matrix(nvdi_cl, ncol = 3, byrow = TRUE)
nvdi_cl_mat
nvdi_recl <- reclassify(nvdi, nvdi_cl_mat)

tm_shape(nvdi_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("grey", "green"), labels = c("Unsuitable (0.5 & Below)", "Suitable (> 0.5)")),
        col.legend = tm_legend(title = "NDVI (Vegetation)", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

# prec
prec

prec_cl <- c(
    -Inf, 350, 0, 
    350, Inf, 1
    )

prec_cl_mat <- matrix(prec_cl, ncol = 3, byrow = TRUE) 
prec_cl_mat

prec_recl <- reclassify(prec, prec_cl_mat)

tm_shape(prec_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("grey", "skyblue"), labels = c("Unsuitable (350mm & Below)", "Suitable (> 350mm)")),
        col.legend = tm_legend(title = "Precipitation (mm)", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

# popl
popl

popl_cl <- c(
    -Inf, 0, 0, 
    0, Inf, 1
    )

popl_cl_mat <- matrix(popl_cl, ncol = 3, byrow = TRUE)
popl_cl_mat
popl_recl <- reclassify(popl , popl_cl_mat)

tm_shape(popl_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("grey", "orange"), labels = c("Unsuitable (0 people)", "Suitable (at least 1 person)")),
        col.legend = tm_legend(title = "Population density", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

# arid
arid

arid_cl <- c(
    -Inf, 0.20, 0, 
    0.20, Inf, 1
    )

arid_cl_mat <- matrix(arid_cl, ncol = 3, byrow = TRUE)
arid_cl_mat
arid_recl <- reclassify(arid, arid_cl_mat)

tm_shape(arid_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("grey", "orange"), labels = c("Unsuitable (0.2 & below)", "Suitable (> 0.2)")),
        col.legend = tm_legend(title = "Aridity (Dryness)", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

1.2.2 Finding areas of suitability by multiplying the binary rasters

We can identify the environmentally suitable areas for occurrence of LF by simply multiplying the binary maps. Therefore, only the cells or areas with the value of 1 will be kept in the output raster layer. You can do this by writing the following formula:

Suitable_LF_Binary <- temp_recl*nvdi_recl*prec_recl*elev_recl*popl_recl*arid_recl

Visualizing the output:

tm_shape(Suitable_LF_Binary) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("#f0f0f0", "red"), labels = c("Zone: Not Suitable", "Zone: Highly Suitable")),
        col.legend = tm_legend(title = "Suitability Map (Binary)", frame = FALSE)) +
tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_text("NAME_1", size = "AREA") +
tm_layout(frame = FALSE)

Red region corresponds to areas that are highly suitable for LF. South East of Kenya (i.e., Kwale, Kilifi and Lamu) are environmental suitable based on binary classification

1.2.3 Finding areas of suitability by summing the binary rasters

We can identify the environmentally suitable areas for occurrence of LF by simply summing the binary maps generate ratings for each pixel to show the intensity of suitability for LF. A pixel with a summed value would mean the following:

• 0 would mean that none of the six factors are present at pixel.
• 1 would mean at least one out of the six factors are present at that pixel.
• 2 would mean at least two out of the six risk factors are present at that pixel (and so on…)
• 6 is the maximum meaning that all factors are present there. Therefore, pixel is rated as the highest levels of suitability for LF.

To sum up the raster, we would first need to use the stack() function before using the calc() for the summation.

RasterStack <- stack(temp_recl, nvdi_recl, prec_recl, elev_recl, popl_recl, arid_recl)
Suitable_LF_Summed <- calc(RasterStack, sum)
# check for minimum and maximum
Suitable_LF_Summed@data@min
Suitable_LF_Summed@data@max
# minimum = 2
# maximum = 6

Visualizing the output:

tm_shape(Suitable_LF_Summed) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("#fff5f0", "#FADBD8", "#F5B7B1", "#F1948A", "#E74C3C"), 
        labels=c("Low (2)", "Modest (3)", "Medium (4)", "High (5)", "Highest (6)")),
        col.legend = tm_legend(title = "Suitability Map (Summation)", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_text("NAME_1", size = "AREA") +
    tm_layout(frame = FALSE)

Note

Using this approach allows the user to see of the fly the intensity for LF suitability. States such as Mandera, Wajir, Isiolo and Turkana have areas that have largely low or modest ratings for LF suitability. Whereas, again the highly suitable areas for LF are the South East states (i.e., Kwale, Kilifi and Lamu).

The rating approach is much better because you can see the intensity for LF suitability. Both approaches highly accessible especially for data and public health program managers in sub-Saharan Africa. These kinds of maps are particularly useful to optimize resources (especially within a low-resource setting) for identifying high-risk areas through knowledge-driven approach (and averting to survey potential unsuitable areas to not wasting limited resources).

This approach is great for descriptive knowledge-based decision analysis; however, it does not take into account the factor that certain variables have more importance than others. For instance, precipitation is significantly more important than NDVI because elevated levels of rainfall has a positive impact on a mosquitoes’ survival lifespan and it provides an reservoir for female mosquitoes to breed in standing water. Higher rainfall within a densely populated environment provides not only a breeding ground, but a location closer to their food source (i.e., human and animal blood). Thus, population density is important an factor above NDVI but not as important as precipitation (i.e., high importance: precipitation > population density > NDVI :less importance).

We can take this into account by applying weights for each variable determined by importance for the MCDA process - this approach is referred to as the Saaty’s Analytical Hierarchy Process (AHP). Let us see how we apply this methodology.

1.3 Analytical Hierarchy Process (AHP)

1.3.1 Standardizing the raster variables to the same scale (from 1 to 10)

We are going to use Saaty’s Analytical Hierarchy Process (AHP). We will need to standardize our raster factors in order to make comparisons and combination possible, as all of them contain different measures: temp (degree Celsius), prec (mm), elev (meters), popl (counts/sqkm); while nvdi and arid are derived indices without any units. Before deriving the weights and applying to the equation that is a linear combination of the above variables to estimate the suitability index for LF, we can standardize them on to the same scale.

Note

We can use the Raster reclassification to standardise all rasters on the scale from 0 to 10.

\[ \text{Standardisation} = \frac{R - R_{\min}}{R_{\max} - R_{\min}} \quad (10) \]

Let’s begin to standardize the first variable temp to get the raster scales from 1 to 10.

temp_standardised_10 <- ((temp - temp@data@min)/(temp@data@max-temp@data@min))*10
plot(temp_standardised_10)

Visualize the output with the scale from 1 to 10:

tm_shape(temp_standardised_10) + 
    tm_raster(col.scale = tm_scale_continuous(values = "-brewer.spectral"),
        col.legend = tm_legend(title = "Temperature [Rescaled]", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

We need to repeat this process of standardisation for the remaining raster grids for prec, popl and ndvi. The solutions are provided in the hidden code chunks below.

Click here

Solutions

# precipitation
prec_standardised_10 <- ((prec - prec@data@min)/(prec@data@max-prec@data@min))*10
plot(prec_standardised_10)

# population
popl_standardised_10 <- ((popl - popl@data@min)/(popl@data@max-popl@data@min))*10
plot(popl_standardised_10)

# normalised vegetation difference index
nvdi_standardised_10 <- ((nvdi - nvdi@data@min)/(nvdi@data@max-nvdi@data@min))*10
plot(nvdi_standardised_10)

For elev, we will treat this differently - of course, we’ll apply the standardisation. But for elevation, the risk of LF decreases with higher values for elevation. Therefore, we need to flip the raster values accordingly to invert the raster to appear as increasing values reflect that elevation is decreasing

# see what the original raster looks like
plot(elev)

# standardise the elevation
elev_standardised_10 <- ((elev - elev@data@min)/(elev@data@max-elev@data@min))*10
plot(elev_standardised_10)

# clone the standardised elevation layer and rename it
flipped_elev_standardised_10 <- elev_standardised_10
# invert the raster simply that taking the maximum value from it and subtract it by its pixel value
values(flipped_elev_standardised_10) <- elev_standardised_10@data@max - values(elev_standardised_10)
# view inverted raster
plot(flipped_elev_standardised_10)

Visualize the inverted output (NOTE: Blue: High elevation, Red: low elevation):

tm_shape(flipped_elev_standardised_10) + 
    tm_raster(col.scale = tm_scale_continuous(values = "-brewer.spectral"),
        col.legend = tm_legend(title = "Inverted Elev (Rescale)", frame = FALSE)) +
    tm_shape(kenya_states) + 
    tm_polygons(fill_alpha = 0, col = "black") +
    tm_layout(frame = FALSE)

Note that for arid, we are going to treat it as a constraint (i.e., binary). Anything below 0.2 threshold are environments considered as hyperarid (extremely dry) areas and thus mosquitoes are non-existent there (- meaning that there is no point for us to include that part of the map with hyperarid environments in the suitability analysis!). We will simply mask out that entire region with any pixel contain that value 0.

Let’s create that masking layer:

# to filter only the values == 1
arid_recl[arid_recl[] !=0 ] = NA
plot(arid_recl)

1.3.2 Determining the weight of each variable using Saaty’s AHP

To estimate the weights, a helpful tool has been developed and provided in an Excel Worksheet. It has been created for you so that you can see how the calculations are carried out step-by-step to derive the weights. You can interact with the cells to see the formulas and explanations are provided at each step.

Open the spreadsheet Guidance - AHP Calculator Tool.xlsx and read to follow through the steps carefully. These are steps summarized here:

• Pairwise comparison & Judgement values: Carry out a pairwise comparison of the factors by constructing a matrix as seen in step 1 in the spreadsheet. Use the criteria and make a “Judgement” on the pairwise variables/factors (it could be based on expert opinion or from literature) asserting a variable’s importance over the other.
• Matrices: Create a square matrix with the judgement values inserted to the matrix accordingly (see step 2)
• Reciprocals: Take the reciprocals in the matrix (see step 3)
• Take column sums: Sum each of the columns accordingly (see step 4)
• Normalization of matrix: For each element in a column of the square matrix, divide it by it’s corresponding column sum. Repeat this step for all other elements in that matrix (see step 5).
• Priority Vectors/Weights: Once completed, sum-up the elements across each row and divide it by the number of variables to obtain the priority vector or weights. As a sanity check, you will know that the calculations are correct if the sum of the weights are equal to 1 (see step 6).
• Validation of whether the judgement values are reasonable: We need to calculate a Consistency Ratio (CR), which is derived from the Consistency Index (CI) divided by the Random Index (RI). For the CI, we must estimate an eigenvalue which is derived from the summed products between the summed column and weights (see step 7 and click on the cell E94 to view the formula). Use that eigenvalue (i.e., Lambda_Max) and estimate the CI (see formula in the lecture notes (slide 25), and click on the cell E95). Next, use the Random Index table (developed by Saaty, 1980) to determine the RI based on the number of factors (in this case, it is 5). Finally, calculate the CR by dividing the CI/RI (see step 7 and click on the cell E111 to view the formula). Note that if the CI < 0.1, the judgement values assigned in step 1 were acceptable. If CR is bigger than 0.1, then the judgement values in the pairwise comparison in step 1 were unreasonable (and thus you will have to repeat the step with different values again until you get an acceptable CR estimate that is below 0.1!). Here, our CR is 0.0351 < 0.1.

1.3.3 Weighted Linear Combination (WLC) and derivation of Suitability Maps based from AHP analysis

Our model uses the Weighted Linear Combination (WLC) approach as the decision rule. The formula to estimate the suitability of LF is as follows:

Suitability (LF) = \(w_1x_1 + w_2x_2 + w_3x_3 + w_4x_4 + w_5x_5 + w_6x_6\)

The \(w_i\) and \(x_i\) represent the weights deprived from AHP analysis and raster variables, respectively. The weights are multiplied to its corresponding raster to obtain a raster with values weighted but ultimately scaled with an upper limit of 10.

This is the derived formula to use in calculating the suitability regions for LF in RStudio:

# use the rescaled columns in the formula
suitablemap_WLC <- 0.372*prec_standardised_10 + 0.356*temp_standardised_10 + 0.159*popl_standardised_10 + 0.077*flipped_elev_standardised_10 + 0.037 * nvdi_standardised_10
# see information
suitablemap_WLC

Important

EXPLANATION

Component	Max Contribution (approx)
0.372 × 10 (precipitation)	3.72
0.356 × 10 (temperature)	3.56
0.159 × 10 (population)	1.59
0.077 × 10 (flipped elevation)	0.77
0.037 × 10 (NDVI)	0.37
Total Theoretical Max	10.01

The observed max ≈ 6.786, after taking the weighted sums, meaning that no pixel achieved the maximum 10 across all factors simultaneous. So, 6.786 is highest achievable weighted suitability given your dataset. Likewise, 2.81112 is the lowest achieveable weighted suitability from that given dataset.

Optional: Re-standardised weighted suitability map to 0-10, add the constrained by aridity layer, and finally apply a criteria to help with the interpretation of risk

Range	Category
0.00 – 0.20	Very Low
0.20 – 0.40	Low
0.40 – 0.60	Moderate
0.60 – 0.80	High
0.80 – 1.00	Very High

Final output:

# re-standardised weighted suitability map, constrained by aridity layer
SWLC_standardised_10 <- ((suitablemap_WLC - suitablemap_WLC@data@min)/(suitablemap_WLC@data@max-suitablemap_WLC@data@min))*10
plot(SWLC_standardised_10)

# create the classification i.e., risk labels
SWLC_cl <- c(
    -Inf,0,0,
  0,2,1,
    2,4,2, 
    4,6,3, 
    6,8,4, 
    8,10,5
    )

SWLC_cl_mat <- matrix(SWLC_cl, ncol = 3, byrow = TRUE)
SWLC_cl_mat
SWLC_recl <- reclassify(SWLC_standardised_10, SWLC_cl_mat)

table(SWLC_recl@data@values)

# code for producing the final output
tm_shape(SWLC_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = c("white", "#74add1", "#e0f3f8", "#ffffbf", "#fee08b", "#f46d43"), 
        labels = c("Not Suitable (0)", "Very Low (1)", "Low (2)", "Moderate (3)", "High (4)", "Very High (5)")),
        col.legend = tm_legend(title = "Suitability Map (MCDA-AHP)", frame = FALSE)) +
tm_shape(arid_recl) + 
    tm_raster(col.scale = tm_scale_categorical(values = "white"), 
        col.legend = tm_legend(show = FALSE)) + 
    tm_shape(kenya_states) + 
        tm_polygons(fill_alpha = 0, col = "black") +
    tm_text("NAME_1", size = 0.8) +
    tm_layout(frame = FALSE)

Note

Areas with high suitability for mosquito occupancy based on the optimal conditions (from selected environmental variables and constrained on aridity) are locations in the far south-east region, small pockets in the central region, and far into the deep-west region of Kenya.

1.4 References (see reading list)

1.4.1 Recommended reading

1. Paper: R.W. Saaty (1987). The Analytical Hierarchy Process - What it is and how it is used, Mathl Modelling, Volume 9(3-5), pages 161 to 176 LINK
2. Technical Document: IDRISI Kilimanjaro: Guide to GIS and Image Processing. Chapter 12: Decision Support - Decision Strategy Analysis, Section on ‘A Typology of Decisions’, pages 151 to 154. DOWNLOAD
3. Paper: A.A. Ali, A. Ahmad (2019). Spatial susceptibility analysis of vector-borne diseases in KMC using geospatial technique and MCDM approach, Modeling Earth Systems and Environment, Issue 5, pages 1135–1159 LINK
4. Paper: X. Zhao et al., (2020). Malaria Risk Map Using Spatial Multi-Criteria Decision Analysis along Yunnan Border During the Pre-elimination Period, Am J Trop Med Hyg, Volume 103 (2), pages 793–809 LINK
5. Paper: B. Ahmed (2015). Landslide susceptibility mapping using multi-criteria evaluation techniques in Chittagong Metropolitan Area, Bangladesh. Landslides, Issue 12, pages 1077–1095 LINK
6. Book: [R Programming] Dorman, M. (2014) Learning R for Geospatial Analysis; Chapter 3: Working with Rasters Click link (Note: Downloadable)

1.4.2 Data Sources

1. The population density 2015 data for Kenya was obtained from Worldpop and resampled at 1km Click Here
2. Raster for annual precipitation was created from bio1_27 and bio1_37, averaged and clipped to Kenya. Obtained for WorldClim Click Here
3. Raster for annual temperature was created from bio1_27 and bio1_37, averaged and clipped to Kenya. Obtained for WorldClim Click Here
4. Aridity Index data was obtained for the Global Aridity and PET Database and clipped to Kenya Click Here
5. Normalized Differenced Vegetation Index was obtained from NASA MODIS Click Here (Registration required).
6. Elevation was obtained from the SRTM CSI CGIAR Project, and cropped to Kenya Click Here