---
title: "Temperature-Humidity Index @ Kaunas"
output: 
  flexdashboard::flex_dashboard:
    orientation: columns
    source_code: embed
    vertical_layout: fill
---

```{r setup, include=FALSE}

# Install / load required R libraries
if (!"pacman" %in% rownames(installed.packages())) install.packages("pacman")
pacman::p_load(dplyr, tidyr, ggplot2, plotly, dygraphs)
pacman::p_load(httr, data.table, dtplyr, xts, eurostat, pals, pbapply, doParallel)
euCodes <- eurostat::harmonize_country_code(eurostat::eu_countries$code)
pbapply::pboptions(type = "txt", style = 3)
cpuCores <- parallel::detectCores()
cl <- parallel::makeCluster(cpuCores)
doParallel::registerDoParallel(cl)

# Load helper functions
source("function_urlDownload.R")
source("function_xts_heatmap.R")

```

```{r download-new, include=FALSE}

# Download the newest data
urlLocation <- 'https://aqicn.org/data-platform/covid19/report/39374-7694ec07/'
csvFile <- urlDownload(urlLocation, showProgress = TRUE)

```

```{r download-old, include=FALSE}

# Create lists of year and quarter names
yNames <- as.character(2019:2021)
qNames <- paste0("Q", 1:4)

# Create a data frame with the url location and year/quarter combinations
# Combine url location and year/quarter combinations into a single column
DF <- data.frame(urlLocation, crossing(yNames, qNames)) %>%
  unite(urlLocations, c(urlLocation, "yNames", "qNames"), sep="")

# Download legacy data (in parallel)
csvFiles <- pbapply::pblapply(DF$urlLocations, urlDownload, cl = cl)

```

```{r preprocess, include=FALSE}

# Define the columns to load
statCols <- c('count', 'min', 'max', 'median', 'variance')
mainCols <- 'median'
dropCols <- setdiff(statCols, mainCols)

# Read the newest data file and skip the first 4 lines
newTable <- data.table::fread(csvFile, skip = 4, drop = dropCols, showProgress = TRUE)

# Leave EU data, rename main column to Value
newTable <- newTable %>%
  rename(Value=all_of(mainCols)) %>%
  filter(Country %in% euCodes)

# Read legacy data files (in parallel)
parallel::clusterExport(cl, c("dropCols"))
parallel::clusterEvalQ(cl, pacman::p_load(data.table))
fileNamesQ <- list.files('.', '^.*Q\\d.csv$')
oldTable <- do.call(rbind, pbapply::pblapply(fileNamesQ, function(x) fread(x, skip = 4, drop = dropCols), cl = cl))
oldTable <- oldTable %>%
  rename(Value=all_of(mainCols)) %>%
  filter(Country %in% euCodes)

# Append old (2018-2021) and new (2022-2023) data tables, sort, remove duplicates
dataTableEU <- oldTable %>% bind_rows(newTable) %>% arrange(Country, City, Date) %>%
  unique()

# Descriptive statistics for daily values of selected variables
selectedVars <- c('temperature', 'humidity')
dataTableEU <- dataTableEU %>% filter(Specie %in% selectedVars)

# 2021-10-03 Barcelona fix
dataTableEU <- dataTableEU %>% group_by(Date, Country, City, Specie) %>% summarize(Value=mean(Value), .groups="drop")

# Create pivot table, calculate THI for each row, drop rows with missing THI values
selectedVars <- c('temperature', 'humidity')
dataTableTHI <- dataTableEU %>%
  pivot_wider(names_from=Specie, values_from=Value) %>%
  group_by(Date, Country) %>%
  mutate(across(all_of(selectedVars), ~ ifelse(is.na(.), median(.,na.rm=TRUE), .))) %>%
  mutate(THI = 0.8*temperature + (humidity/100)*(temperature-14.4) + 46.4) %>%
  filter(!is.na(THI)) %>%
  data.table

# Prepare data for selected city
myCity <- 'Kaunas'
myTitle <- paste("THI @", myCity)

myTable <- dataTableTHI %>% filter(City==myCity, Date>='2019-01-01') %>%
  select(Date, humidity, temperature, THI) %>% mutate(across('THI', round, 2))

xts_data <- xts(myTable[,c('temperature','humidity','THI')], order.by = as.Date(myTable$Date))
saveRDS(xts_data, 'waqi-covid19-airqualitydata.rds')

```

Calplot (static)
=======================================================================

```{r static}

# Plot THI time series calendar heatmap
xts_heatmap(xts_data[,'THI']) + labs(title=myTitle) + coord_fixed(1) + theme_minimal(base_size=7)

```

Calplot (dynamic)
=======================================================================

```{r dynamic}

# Plot THI time series calendar heatmap
gg <- xts_heatmap(xts_data[,'THI']) + labs(title=myTitle) + coord_fixed(1) + theme_minimal(base_size=12)
ggplotly(gg)

```

Series
=======================================================================

### **Temperature (°C)**

```{r temp}

# Plot temperature time series dynamics
dygraph(xts_data[,'temperature'], group = "Air-Data") %>%
 dySeries(c("temperature"), label = 'Temperature') %>%
 dyCrosshair(direction = "vertical")
```

### **Humidity (%)**

```{r humid}

# Plot humidity time series dynamics
dygraph(xts_data[,'humidity'], group = "Air-Data") %>%
 dySeries(c("humidity"), label = 'Humidity') %>%
 dyCrosshair(direction = "vertical")

```

### **Temperature-Humidity Index (THI)**

```{r THI}

# Plot THI time series dynamics
dygraph(xts_data[,'THI'], group = "Air-Data") %>%
 dySeries("THI") %>%
 dyRangeSelector() %>% dyCrosshair(direction = "vertical")

```


About
=======================================================================

Air Quality Open Data Platform and World Air Quality Index project team provides a [new dedicated dataset](https://aqicn.org/data-platform/covid19/), updated 3 times a day, and covering about 380 major cities in the world, from January 2015 until now. The data for each major city is based on the average (median) of several air quality monitoring stations. The data set provides min, max, median and standard deviation for each of the air pollutant species (PM10, PM2.5, NO2, Ozone, SO2, CO) as well as meteorological data (Temperature, Humidity, Wind).

THI measures the level of discomfort the average person (or animal) is thought to experience as a result of the combined effects of the temperature and humidity of the air. The heat stress index THI is a simple combination of temperature and humidity (in decimal form) and can be calculated by the following formula:

$THI = 0.8 \times Temperature + Humidity \times (Temperature - 14.4) + 46.4$

Habeeb, A. A., Gad, A. E., & Atta, M. A. (2018). Temperature-humidity indices as indicators to heat stress of climatic conditions with relation to production and reproduction of farm animals. *International Journal of Biotechnology and Recent Advances*, 1(1), 35-50. [https://doi.org/10.18689/ijbr-1000107](https://doi.org/10.18689/ijbr-1000107)

### **Descriptive statistics for THI**

```{r table}

knitr::kable(do.call(rbind, tapply(myTable$THI, year(myTable$Date), summary)), digits=2)

```