MalariaIndicator.Rmd

---
title: "Malaria Suitability Indicator"
author: "Earth Sciences department, Barcelona Supercomputing Center (BSC)"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::knitr}
  %\VignetteIndexEntry{Agricultural Indicators}
  %\usepackage[utf8]{inputenc}
---

<style>
body {
text-align: justify
}
</style>

```{r echo=FALSE}
knitr::opts_chunk$set(warnings=FALSE, message=FALSE)
```

Climate Suitability for Malaria Transmission
-----------------------------


## Introduction

Climate sensitive infectious diseases are a subject of major concern for public 
health around the world in the context of global change. Malaria is an infectious 
disease caused by parasites of the genus _Plasmodium_ and transmitted by _Anopheles_ 
mosquitoes (Boualam et al., 2021). Eradication efforts managed to achieve complete 
elimination of endemic  malaria in Europe since 1974, although sporadic transmission 
events are frequently reported in travelers coming from other endemic areas. Studies 
evaluating the factors that keep malaria out of the continent found that high socio-economic 
conditions and increased life expectancy are key (Fischer et al., 2020). On the 
other hand, climatic conditions are becoming more suitable for malaria 
transmission (Zhao et al., 2016). 

The [Lancet Countdown in Europe collaboration]("https://www.lancetcountdown.org/europe/") 
was launched in 2021 with the purpose of using a variety of indicators to monitor 
trends in the impact of climate on the dynamic of different health processes. As 
part of thisinitiative, the Climate Suitability for Malaria Transmission indicator 
was adapted from the Global Lancet Countdown Collaboration (Romanello et al., 2021) 
for tracking the number of months suitable of malaria transmission per year in the 
continent between 1950 and 2021.This indicator is a threshold based model that 
overlaps the climatic and environmental requirements of _Anopheles_ mosquitoes and 
_Plasmodium vivax_, which was the parasiteof endemic circulation in the continent 
until the mid 1970s (Boualam et al., 2021).

Although the `CSIndicators` was initially developed for providing a set of tools
useful in the context of agricultural applications, it is also possible to expand
its applications to the computation of other type of indicators. In this vignette,
the Climate Suitability  for Malaria indicator will be computed using  the `TotalTimeExceedingThreshold` 
function. The output will be an multidimensional array containing the number of 
months per year suitable for malaria transmission.


## Load libraries

In addition to the `CSIndicators` package, the `CSTools` and `s2dv` packages will
be used for data extraction and visualization.

```{r eval=FALSE}
library(CSTools)
library(CSIndicators)
library(s2dv)
```


```{r echo=FALSE}
#Temporary links to functions
library(CSIndicators)
source("https://earth.bsc.es/gitlab/nperez/Flor/-/raw/master/CST_Load_devel_from_s2dv.R")
```


## 1. Data
### 1.1 Area of interest
The number of months suitable for malaria transmission between 2010 and 2020 will
be computed for the Iberian Peninsula (40º14´24´´N 4º14´21´´W). This area comprises
the territories of Portugal, Spain, France, Andorra and Gibraltar.

```{r echo=FALSE}
knitr::include_url("https://es.wikipedia.org/wiki/Pen%C3%ADnsula_ib%C3%A9rica#/media/Archivo:Espa%C3%B1a_y_Portugal.jpg")
```

### 1.2 Data extraction
Malaria indicator is computed with climate and environmental data. Climate variables
are extracted from the [ERA5-Land reanalysis dataset]("https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=overview") (Muñoz-Sabater et al., 2021), which span between 1950
and present day, at a 9km resolution. Necessary climate variables are 2 meter 
temperature (`tas`), 2 meter dew point temperature (`tdps`) and total precipitation
 (`prlr`). Environmental information is extracted from the [CORINE Land Cover inventory]("https://land.copernicus.eu/pan-european/corine-land-cover"), 
which classifies land into 44 classes at a 100 m resolution. These data are produced 
by the Copernicus service and extracted and reprocessed by the **Earth Sciences department** 
from the **Barcelona Supercomputing Center**.

The function `CST_Load` from the `CSTools` package allows the user to access the
processed data and delivers the data in the form of an `s2dv_cube` object. These
multidimensional objects store data and metadata in several elements. For further 
information about these objects visit [Data retrieval and storage](https://cran.r-project.org/package=CSTools/vignettes/Data_Considerations.html). 

`CST_Load` has its requirements for extracting data. First, the path to the files 
is introduced using *whitecards* delimited by dollar signs. Second, the temporal 
extent and resolution of the output con be provided with `sdates`, which takes a 
vector of all dates of interest. The spatial extent is inputted within `lonmax`, 
`lonmin`, `latmax` and `latmin.`

```{r}
# Path to file locations using whitecards
path_ERA5_CDS <- list(path='/esarchive/recon/ecmwf/era5land/$STORE_FREQ$_mean/$VAR_NAME$_f1h/$VAR_NAME$_$YEAR$$MONTH$.nc')

# Temporal extent and resolution:
# We will work with all months from 2010 to 2020
year_in <- 2020 # initial year
year_fi <- 2021 # last year
month_in <- 1 # first month
month_fi <- 12 # last month

sdates <- paste0(year_in:year_fi, '0101')

# Extract data
vars <- c('tas', 'tdps', 'prlr')

out <- NULL
for(var in vars){
  out[[var]] <- CST_Load_s2dv(var = var,
                           exp = NULL, # We only require observed data
                           obs = list(path_ERA5_CDS),
                           sdates = sdates,
                           lonmax = 355, lonmin = 350,
                           latmax = 45, latmin = 43,
                           storefreq = 'monthly',
                           leadtimemin = month_in, 
                           leadtimemax = month_fi,
                           output = "lonlat")  
}

# library(zeallot)
# c(tas, tdps, prlr) %<-% CST_Load_s2dv(var = c('tas', 'tdps', 'prlr'),
#                                       exp = NULL, # We only require observed data
#                                       obs = list(path_ERA5_CDS),
#                                       sdates = sdates,
#                                       lonmax = 355, lonmin = 350,
#                                       latmax = 45, latmin = 43,
#                                       storefreq = 'monthly',
#                                       leadtimemin = month_in, 
#                                       leadtimemax = month_fi,
#                                       output = "lonlat")  
```

The call above creates a list with the three datasets. Each `s2dv_cube` object has
the following dimensions:
```{r}
dim(out$tas$data)
```
And it is possible to have a look at the data:
```{r}
# A summary of the entire dataset
summary(out$tas$data)
```

```{r}
# Sneak peek of a sample of the temperature data in Jan 2010
out$tas$data[1,1,1,1,,][15:20, 15:20]
```

## 2. Data transformation
### 2.1 Temperature and dew point temperature
This variables will be used to compute relative humidity, which requires them to 
be in Degrees Celsius (C):
```{r}
for(var in c('tas', 'tdps')){
  out[[var]]$data <- out[[var]]$data - 273.15
}

summary(out$tas$data)
```

### 2.3 Relative humidity
After transforming temperature ($T$) and dew point temperature ($T_d$) units from Kelvin (K) to
Degrees Celsius (C), it is possible to compute the relative humidity using 
August-Roche-Magnus equation (Alduchov 1996):

$$
RH=\frac{100*exp(\frac{aT_d}{b+T_d})}{exp(\frac{aT}{c+T})}
$$
Where _a_ and _b_ are the coefficients 17.625 and 243.04, respectively.

```{r}
# Create a template with the same dimensions as the datasets
out$rh <- out$tas

# Compute relative humidity (RH)
out$rh$data <- 100*(exp((17.625*out$tdps$data)/(243.04+out$tdps$data)) / 
                      exp((17.625*out$tas$data) / (243.04+out$tas$data)))

summary(out$rh$data)
```

### 2.2 Precipitation
Precipitation is extracted from the hourly reanalaysis dataset and processes as
monthly meters per second. For the computation of the malaria indicator, it is
necessary to transform it into accumulated mm per month:
```{r}
# To compute the accumulated precipitation per month it is necessary to multiply by:
# - total number of seconds in one hour: 3600
# - total number of hours in a day: 24
# - average number of days in a month: 30.41
# - change factor from meter to mm: 1000
out$prlr$data <- out$prlr$data*3600*24*30.41*1000

summary(out$prlr$data)
```

### 2.4 Land use
Add chunk

## 3. Compute indicator
The suitability for malaria transmission indicator estimated the number of months
suitable for transmission of _P. vivax_, calculated from the following thresholds:

* Accumulated monthly precipitation above 80 mm (Lyon et al., 2017),
* Monthly temperature between 14.5ºC and 33ºC (Grover-Kopec et al., 2006), and
* Monthly relative humidity above 60% (Grover-Kopec et al., 2006)

Additionally, _Anopheles_ suitability is influenced by the surrounding conditions, so that the
following land cover classes are considered highly suitable (Benali et al., 2014):

* Rice fields,
* Permanently irrigated croplands, and
* Sport and leisure facilities

```{r}
# Create a template with the same dimensions as the datasets
malariaSuit <- out$tas

# Determine suitability
malariaSuit$data <- ifelse(out$tas$data >= 14.5 & 
                             out$tas$data <= 33 & 
                             out$rh$data >= 60 & 
                             out$prlr$data >= 80, 1, 0)

# Compute number of number of months per year that are suitable
# (i.e. suitability = 1)
malariaInd <- CSIndicators::CST_TotalTimeExceedingThreshold(malariaSuit, 
                                                            threshold=0.5)

# Sneak peek of the data
malariaInd$data[1,1,1,,][15:20, 15:20]
```

## 4. Visualization
### 4.1. Map
Using the function `PlotEquiMap()` from the `s2dv` package, it is possible to 
have a visual inspection of the number of suitable months in January 2010.

```{r echo=FALSE}
# Until the dependency with s2dverification in CSTools is updated, it is necessary 
# to make these adjustments
source('https://earth.bsc.es/gitlab/es/s2dv/-/raw/master/R/PlotEquiMap.R')
.FilterUserGraphicArgs <- s2dv:::.FilterUserGraphicArgs
.KnownLonNames <- s2dv:::.KnownLonNames
.KnownLatNames <- s2dv:::.KnownLatNames
ColorBar <- s2dv::ColorBar
clim.palette <- s2dv::clim.palette
.IsColor<- s2dv:::.IsColor
```

```{r fig.width=8, fig.height=4}
PlotLayout(PlotEquiMap, c('lat', 'lon'),
           var=malariaInd$data[1,1,,,],
           nrow=1, ncol=2,
           lon=malariaInd$lon,
           lat=malariaInd$lat,
           filled.continents=FALSE,
           brks=seq(0, 12, by=1),
           toptitle="Number of months per year suitable for malaria in the Iberian Peninsula",
           titles=c("2010", "2011"),
           title_scale=0.5,
           coast_width=2,
           filled.oceans=TRUE,
           country.borders=TRUE,
           intylat=1,
           intxlon=1)

```

### 4.2 Time series
Add chunk

## 5. Summary statistics
Add chunk

## 6. Conclusions
Add chunk

## References

* Alduchov O. Improved Magnus Form Approximation of Saturation Vapor Pressure. _J Appl Meteorl Climatol. 1996; 35(4),601-609._ [URL]("https://journals.ametsoc.org/view/journals/apme/35/4/1520-0450_1996_035_0601_imfaos_2_0_co_2.xml")
* Benali A, Nunes JP, Freitas FB, et al. Satellite-derived estimation of environmental suitability for malaria vector development in Portugal. _Remote Sens Environ. 2014;145:116-130._ DOI:10.1016/J.RSE.2014.01.014
* Boualam MA, Pradines B, Drancourt M, Barbieri R. Malaria in Europe: a historical perspective. _Front Med. 2021; 8:876._ DOI:10.3389/FMED.2021.691095/BIBTEX
* Fischer L, Gültekin N, Kaelin M, Fehr J, Schlagenhauf P. Rising temperature and its impact on receptivity to malaria transmission in Europe: a systematic review. _Travel Med Infect Dis. 2020;36:101815._ DOI:10.1016/J.TMAID.2020.101815
* Grover-Kopec EK, Blumenthal MB, Ceccato P, Dinku T, Omumbo JA, Connor SJ. Web-based climate information resources for malaria control in Africa. _Malar J. 2006;5(1):1-9._ DOI:10.1186/1475-2875-5-38/FIGURES/5
* Lyon B, Dinku T, Raman A, Thomson MC. Temperature suitability for malaria climbing the Ethiopian Highlands. _Environ Res Lett. 2017;12(6):064015._ DOI:10.1088/1748-9326/AA64E6
* Muñoz-Sabater J, Dutra E, Agustí-Panareda A, et al. ERA5-Land: A state-of-the-art global reanalysis dataset for land applications. _Earth Syst Sci Data. 2021;13(9):4349-4383._ DOI:10.5194/ESSD-13-4349-2021
* Romanello M, McGushin A, Di Napoli C, et al. The 2021 report of the Lancet Countdown on health and climate change: code red for a healthy future. _Lancet. 2021;398(10311):1619-1662._ DOI:https://doi.org/10.1016/S0140-6736(21)01787-6
* Zhao X, Smith DL, Tatem AJ. Exploring the spatiotemporal drivers of malaria elimination in Europe. _Malar J. 2016;15(1):1-13._ DOI:10.1186/S12936-016-1175-Z/TABLES/5