summarize_data.Rmd

---
title: "Summary of the state of the GPS data set"
csl: the-american-naturalist.csl
output:
  html_document:
    theme: cerulean
    toc: yes
  pdf_document:
    toc: yes
<!-- bibliography: references.bib -->
editor_options: 
  chunk_output_type: console
---

<!--
IMAGES:
Insert them with: ![alt text](image.png)
You can also resize them if needed: convert image.png -resize 50% image.png
If you want to center the image, go through HTML code:
<div style="text-align:center"><img src ="image.png"/></div>

REFERENCES:
For references: Put all the bibTeX references in the file "references.bib"
in the current folder and cite the references as @key or [@key] in the text.
Uncomment the bibliography field in the above header and put a "References"
title wherever you want to display the reference list.
-->

<style type="text/css">
.main-container {
  max-width: 1370px;
  margin-left: auto;
  margin-right: auto;
}
</style>

```{r general options, include = FALSE}
knitr::knit_hooks$set(
  margin = function(before, options, envir) {
    if (before) par(mgp = c(1.5, .5, 0), bty = "n", plt = c(.105, .97, .13, .97))
    else NULL
  },
  prompt = function(before, options, envir) {
    options(prompt = if (options$engine %in% c("sh", "bash")) "$ " else "> ")
  })

knitr::opts_chunk$set(margin = TRUE, prompt = TRUE, comment = "", echo = FALSE,
                      collapse = TRUE, cache = FALSE, autodep = TRUE, message = FALSE,
                      dev.args = list(pointsize = 11), fig.height = 3.5,
                      fig.width = 4.24725, fig.retina = 2, fig.align = "center")

options(width = 137)
```

```{r packages}
library(dplyr)
library(ecomore)
library(lubridate)
library(magrittr)
library(OpenStreetMap)
library(purrr)
library(sf)
```

The CSV file of the **GPS data set** is
[here](https://github.com/ecomore2/gps/blob/master/data/gps.csv). Go
[here](https://raw.githubusercontent.com/ecomore2/gps/master/data/gps.csv)
if you want to copy and paste this CSV file to your computer.

```{r data}
gps <- readr::read_csv("data/gps.csv", col_types = "icdd")
pacs <- readr::read_csv("../pacs/data/pacs.csv",
                        col_types = paste(c("icfnD", rep("c", 5), rep("D", 4), rep("f", 3)), collapse = ""))
```

There are **`r nrow(gps)` coordinates** in this file from
**`r length(unique(gps$id))` unique cases** (i.e.
`r round(100 * length(unique(gps$id)) / nrow(pacs))` % of the total number of
cases, `r nrow(pacs)`, reported in
[PACS](https://ecomore2.github.io/pacs/summarize_data.html)). By year, the GPS
data split like this:

```{r correcting onset dates}
pacs %<>% mutate(onset2 = correct_onset(.))
```

```{r merging gps and pacs}
pacs_gps <- gps %>% 
  mutate(gps = TRUE) %>% 
  select(id, gps) %>% 
  unique() %>% 
  left_join(pacs, ., "id") %>% 
  mutate(gps = ifelse(is.na(gps), FALSE, gps))
```

```{r looking for gps points with no date}
gps_no_dates <- pacs_gps %>% 
  filter(gps, is.na(onset2)) %>% 
  select(id, onset, hospitalization, consultation, sample_collection)
write.csv(gps_no_dates, "problems/gps_no_dates.csv", FALSE, row.names = FALSE)
```

```{r adding year to pacs_gps}
pacs_gps %<>% mutate(year = year(onset2))
```

```{r GPS rate per year}
pacs_gps %>% 
  group_by(year) %>% 
  summarize(n = n(), perc = 100 * mean(gps), gps = sum(gps)) %>% 
  ungroup() %>% 
  select(year, gps, perc)
```

Considering only confirmed cases, it looks like:

```{r adding test variables}
pacs_gps %<>% add_tests()
```

```{r GPS rate per year for confirmed cases}
pacs_gps %>% 
  filter(confirmed) %>% 
  group_by(year) %>% 
  summarize(n = n(), perc = 100 * mean(gps), gps = sum(gps)) %>% 
  ungroup() %>% 
  select(year, gps, perc)
```

The duplicates are

```{r GPS duplicates}
filter(gps, id %in% names(which(table(gps$id) > 1)))
```

There are `r nrow(gps_no_dates)` geocoded cases for which we don't have any
date:

```{r showing gps points with no date}
gps_no_dates
```

The CSV file of these cases is
[here](https://github.com/ecomore2/gps/blob/master/problems/gps_no_dates.csv). Go
[here](https://raw.githubusercontent.com/ecomore2/gps/master/problems/gps_no_dates.csv)
if you want to copy and paste this CSV file to your computer. The split by
source reads

```{r GPS source table}
gps %>% 
  group_by(source) %>% 
  tally() %>% 
  ungroup()
```

The split of data according to the test and the availability of GPS coordinates is:

```{r test GPS table}
addmargins(with(pacs_gps, table(tested, gps)))
```

The cases with GPS data and reported as negative are:

```{r negative with GPS}
pacs_gps %>% 
  filter(tested, ! confirmed, gps) %>% 
  select(id, gps, confirmed)
```

The split of data according to the availability of village data and GPS coordinates is:

```{r table GPS village}
pacs_gps %>% 
  mutate(village_info = !is.na(village)) %$%
  table(village_info, gps) %>% 
  addmargins()
```

```{r loading OpenStreetMap tiles}
upperleft  <- c(18.484562, 101.983290)
lowerright <- c(17.753338, 103.164520)
nb <- 20
if (! file.exists("bing.rds")) {
  bing <- openmap(upperleft, lowerright, type = "bing", minNumTiles = nb)
  saveRDS(bing, "bing.rds")
} else {
  bing <- readRDS("bing.rds")
}
if (! file.exists("osm.rds")) {
  osm <- openmap(upperleft, lowerright, type = "osm", minNumTiles = nb)
  saveRDS(osm, "osm.rds")
} else {
  osm <- readRDS("osm.rds")
}
```

```{r loading GADM polygons}
for(f in c("gadm36_LAO_1_sf.rds", "gadm36_LAO_2_sf.rds"))
  if (! file.exists(f))
    download.file(paste0("https://biogeo.ucdavis.edu/data/gadm3.6/Rsf/", f), f)
lao1 <- readRDS("gadm36_LAO_1_sf.rds")
lao2 <- readRDS("gadm36_LAO_2_sf.rds")
```

```{r points of interest coordinates}
points_of_interest <- list(ipl = c(17.962593, 102.615130),
                           vki = c(17.965726, 102.605642),
                           wai = c(17.975349, 102.568564))
```

```{r Vientiane polygon}
vt <- filter(lao1, NAME_1 == "Vientiane [prefecture]")
```

```{r GPS points}
gps_sf <- gps %>% 
  select(longitude, latitude) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
```

```{r districts polygons}
districts <- filter(lao2, NAME_1 == "Vientiane [prefecture]")
```

```{r POI points}
poi <- points_of_interest %>% 
  lapply(setNames, c("latitude", "longitude")) %>% 
  reduce(bind_rows) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
```

```{r a function to plot GPS on a map}
plot_map <- function(map, border = "black") {
  proj <- map$tiles[[1]]$projection
  plot(map)
  vt %>% 
    st_transform(proj) %>% 
    st_geometry() %>% 
    plot(add = TRUE, border = border)
  gps_sf %>% 
    st_transform(proj) %>% 
    plot(add = TRUE, col = "red")
  districts %>% 
    st_transform(proj) %>% 
    st_geometry() %>% 
    plot(add = TRUE, lty = 3, border = border)
  poi %>% 
    st_transform(proj) %>% 
    plot(add = TRUE, col = "blue", pch = 3)
}
```

Here is a map of the geolocated cases:

```{r OSM map, fig.width = 2 * 4.24725, fig.height = 2 * 3.5}
plot_map(osm)
```

The blue crosses are the Wattay International Airport, the Vayakorn Inn and the
Institut Pasteur du Laos. Same with satellite image background:

```{r BING map, fig.width = 2 * 4.24725, fig.height = 2 * 3.5}
plot_map(bing, "white")
```