---
title: "Heat-health indicators for Bhubaneswar and Jodhpur"
output: 
  flexdashboard::flex_dashboard:
    storyboard: true
    source: embed
---

```{r setup, include=FALSE}
library(flexdashboard)
```

<style type="text/css">
body{
  font-family: Verdana;
  font-size: 16px;
}
/* Headers */
h1,h2,h3,h4,h5,h6{
  font-size: 16pt;
}
</style>

<style type="text/css">
  li {font-size: 16px;}
</style>



### <font size="3"> **Introduction** </font>

This html page is an interactive graphical tool accompanying the paper "Deciphering heat-health indicators for the non-specialist: Data driven evidence and a tool from two Indian cities", authored by Vikram Dayal, Purnamita Dasgupta, Rupa Kumar Kolli, Girika Sharma, Rajib Chattopadhyay, Rajib Dasgupta, Roxy Mathew Koll. The paper has the details about the data and discusses the different indicators. 

You can view the source code for the tool (see top right). 

This page has interactive versions of the static graphs in the paper. You can follow the story by clicking on the tabs above or the arrow to the right.  

If you move the cursor over a graph, you will find some icons in the top right corner as you look at the screen. You can hover and see the underlying data, zoom and select a portion of the data, and reset axes to go back to the original view.  



```{r}
library(tidyverse)
library(patchwork)
library(plotly)
library(xts)
library(zoo)
```


```{r}
Bhub_data <- read_csv("Bhub_data.csv")
Jodh_data <- read_csv("Jodh_data.csv")
```



### <font size="3">  **Daily maximum and minimum temperature in Bhubaneswar in 2020.** </font>

```{r}
ggB1 <- Bhub_data |>
  filter(Date > "2019-12-31") |>
  ggplot() +
  geom_line(aes(x = Date, y = TemperatureCMax)) +
  geom_line(aes(x = Date, y = TemperatureCMin),
            col = "red") +
  labs(y = "Temperature (Degree Celsius)",
       title = "Bhubaneswar") +
  ylim(0, 50) +
  theme_bw()
ggplotly(ggB1)
```

> **This is part of Figure 2 in the paper. Data accessed with the R climate package.**

***


<font size="3"> If we use the cursor to hover over the lines, we can see the date and temperatures for individual days. 

Looking at the graph, we see that: </font>

- The maximum (black line) was above 30 degrees Celsius for much of the year. 

- The minimum (red line) never went below 10 degrees Celsius. 



### <font size="3">  **Daily maximum and minimum temperature in Jodhpur in 2020.** </font>

```{r}
ggJ1 <- Jodh_data |>
  filter(Date > "2019-12-31") |>
  ggplot() +
  geom_line(aes(x = Date, y = TemperatureCMax)) +
  geom_line(aes(x = Date, y = TemperatureCMin),
            col = "red") +
  labs(y = "Temperature (Degree Celsius)" , 
       title = "Jodhpur") +
  ylim(0,50) +
  theme_bw()
ggplotly(ggJ1)
```

> **This is part of Figure 2 in the paper. Data accessed with the R climate package.**

***


<font size="3"> If we use the cursor to hover over the lines, we can see the dates and temperatures for individual days. 

Looking at the graph we see that: 

- The maximum (black line) was above 40 degrees Celsius several days in the year. 

- The minimum (red line) went above 30 degrees Celsius several days in a year. 

Jodhpur is warmer than Bhubaneswar in the summer and colder in the winter. </font>


### <font size="3">  **Daily average relative humidity in Bhubaneswar in 2020.** </font>

```{r}
ggB2 <- Bhub_data |>
  filter(Date > "2019-12-31") |>
  ggplot() +
  #geom_point(aes(x = TemperatureCAvg, y = HrAvg)) +
  geom_line(aes(x = Date, y = HrAvg)) +
  labs(y = "Relative humidity (Percent)", title = "Bhubaneswar") +
  ylim(20, 100) +
  theme_bw()
ggplotly(ggB2)
```

> **This is part of Figure 3 in the paper. Data accessed with the R climate package.**

***

<font size="3"> In Bhubaneswar the average relative humidity stayed largely over 60 percent. It was over 80 percent for a substantial proportion of days. </font>



### <font size="3">  **Daily average relative humidity in Jodhpur in 2020.** </font>

```{r}
ggJ2 <- Jodh_data |>
  filter(Date > "2019-12-31") |>
  ggplot() +
  #geom_point(aes(x = TemperatureCAvg, y = HrAvg)) +
  geom_line(aes(x = Date, y = HrAvg)) +
  labs(y = "Relative humidity (Percent)", title = "Jodhpur") +
  ylim(20, 100) +
  theme_bw()
ggplotly(ggJ2)
```

> **This is part of Figure 3 in the paper. Data accessed with the R climate package.**

***

<font size="3"> In Jodhpur the average relative humidity was below 60 percent for several months. 

Jodhpur had far lower average relative humidity over the year than Bhubaneswar.

</font>

### <font size="3">  **Heat Index in Bhubaneswar in 2020.** </font>

```{r}
# compute (both Bhub and Jodh)
library(weathermetrics)
Bhub_data$HI <- heat.index(t = 
        Bhub_data$TemperatureCAvg,
        rh = Bhub_data$HrAvg,
        temperature.metric = "celsius",
          output.metric = "celsius")
Jodh_data$HI <- heat.index(t = 
        Jodh_data$TemperatureCAvg,
        rh = Jodh_data$HrAvg,
       temperature.metric = "celsius",
        output.metric = "celsius")
```

```{r}
ggHIB <- Bhub_data |>
  filter(Date > "2019-12-31") |>
  ggplot(aes(x = Date, y = HI)) +
  #geom_point() +
  geom_line() +
  #scale_y_continuous() +
  geom_hline(yintercept = 41, col = "red",
             linetype = "dashed") +
  geom_hline(yintercept = 32, col = "blue",
             linetype = "dashed") +
  labs(y = "Heat Index (degree Celsius)", title = "Bhubaneswar") +
  ylim(10,50) +
  theme_bw()
ggplotly(ggHIB)
```

> **This is part of Figure 4 in the paper. Data accessed with the R climate package.**

***

<font size="3"> In Bhubaneswar, for a lot of the year the heat index is between 32 degrees celsius and 41 degrees celsius (Danger category). 

</font>

### <font size="3">  **Heat Index in Jodhpur in 2020.** </font>

```{r}
ggHIJ <- Jodh_data |>
  filter(Date > "2019-12-31") |>
  ggplot(aes(x = Date, y = HI)) +
  #geom_point() +
  geom_line() +
  #scale_y_continuous() +
  geom_hline(yintercept = 41, col = "red",
             linetype = "dashed") +
  geom_hline(yintercept = 32, col = "blue",
             linetype = "dashed") +
  labs(y = "Heat Index (degree Celsius)", title = "Jodhpur") +
  ylim(10, 50) +
  theme_bw()
ggplotly(ggHIJ)
```

> **This is part of Figure 4 in the paper. Data accessed with the R climate package.**

***

<font size="3"> In Jodhpur, fewer days were in the Danger category but the highest value of the heat index was greater than that in Bhubaneswar. 

</font>


### <font size="3">  **Relationship between relative humidity and temperature in Bhubaneswar in 2020.** </font>

```{r}
ggB2B <- Bhub_data |>
  filter(Date > "2019-12-31") |>
  ggplot() +
  geom_point(aes(x = TemperatureCAvg, y = HrAvg),
             shape = 1) +
  #geom_line(aes(x = Date, y = HrAvg)) +
  labs(y = "Relative humidity (Percent)", 
       x = "Temperature (Degree Celsius)", 
       title = "Bhubaneswar") +
  ylim(20, 100) +
  annotate("text", x = c(35.56, 32.22, 30), 
           y = c(65, 85, 100), 
           label = "60", col = "red") +
  theme_bw()
ggplotly(ggB2B)
```

> **This is part of Figure 5 in the paper. Data accessed with the R climate package.**

***

<font size="3"> The outer edge of the scatter of points in 2020 was not that far from the alarming level of 60 degrees apparent temperature using the revised heat index of Lu and Romps (2022). 

### <font size="3">  **Relationship between relative humidity and temperature in Jodhpur in 2020.** </font>

```{r}
ggJ2B <- Jodh_data |>
  filter(Date > "2019-12-31") |>
  ggplot() +
  geom_point(aes(x = TemperatureCAvg, y = HrAvg),
             shape = 1) +
  #geom_line(aes(x = Date, y = HrAvg)) +
  labs(y = "Relative humidity (Percent)", 
       x = "Temperature (Degree Celsius)",
       title = "Jodhpur") +
  ylim(20, 100) +
  annotate("text", x = c(35.56, 32.22, 30), 
           y = c(65, 85, 100), 
           label = "60", col = "red") +
  theme_bw()
ggplotly(ggJ2B)
```

> **This is part of Figure 5 in the paper. Data accessed with the R climate package.**

***

<font size="3"> The outer edge of the scatter of points in 2020 was not that far from the alarming level of 60 degrees apparent temperature using the revised heat index of Lu and Romps (2022). 

</font>

### <font size="3">  **Excess heat factor in Bhubaneswar from 1951 to 2023** </font>

```{r}
# Pack
library(tidyverse)
library(xts)
library(patchwork)
library(gets)
```


```{r}
# Read
tempB_xts <- readRDS("tempB_xts.rds")
tempB_xts_df <- read_csv("tempB_xts.csv")
tempJ_xts <- readRDS("tempJ_xts.rds")
tempJ_xts_df <- read_csv("tempJ_xts.csv")
```

```{r}
#YearMonth
library(lubridate)
tempB_xts_df$Year <- year(ymd(tempB_xts_df$Index))
tempB_xts_df$Month <- month(ymd(tempB_xts_df$Index))
tempB_xts_df$ay <- mday(ymd(tempB_xts_df$Index))
library(lubridate)
tempJ_xts_df$Year <- year(ymd(tempJ_xts_df$Index))
tempJ_xts_df$Month <- month(ymd(tempJ_xts_df$Index))
tempJ_xts_df$ay <- mday(ymd(tempJ_xts_df$Index))
```

```{r}
## ----Decade
tempB_xts_df$Decades <- case_when(
  tempB_xts_df$Year > 1950 & tempB_xts_df$Year < 1960 ~ "1950s",
  tempB_xts_df$Year > 1959 & tempB_xts_df$Year < 1970 ~ "1960s",
  tempB_xts_df$Year > 1969 & tempB_xts_df$Year < 1980 ~ "1970s",
  tempB_xts_df$Year > 1979 & tempB_xts_df$Year < 1990 ~ "1980s",
  tempB_xts_df$Year > 1989 & tempB_xts_df$Year < 2000 ~ "1990s",
  tempB_xts_df$Year > 1999 & tempB_xts_df$Year < 2010 ~ "2000s",
  tempB_xts_df$Year > 2009 & tempB_xts_df$Year < 2020 ~ "2010s",
)
tempJ_xts_df$Decades <- case_when(
  tempJ_xts_df$Year > 1950 & tempJ_xts_df$Year < 1960 ~ "1950s",
  tempJ_xts_df$Year > 1959 & tempJ_xts_df$Year < 1970 ~ "1960s",
  tempJ_xts_df$Year > 1969 & tempJ_xts_df$Year < 1980 ~ "1970s",
  tempJ_xts_df$Year > 1979 & tempJ_xts_df$Year < 1990 ~ "1980s",
  tempJ_xts_df$Year > 1989 & tempJ_xts_df$Year < 2000 ~ "1990s",
  tempJ_xts_df$Year > 1999 & tempJ_xts_df$Year < 2010 ~ "2000s",
  tempJ_xts_df$Year > 2009 & tempJ_xts_df$Year < 2020 ~ "2010s",
)
```

```{r}
## ----t95B
tempB_xts <- tempB_xts["1951/2023"]
taB <- as.xts(coredata(tempB_xts)[,3],
             order.by = index(tempB_xts))
t95B <- quantile(taB, 0.95, na.rm = T)
```

```{r}
# rollB
lagtaB <- as.xts(c(NA, taB[1:(length(taB) - 1)]),
                order.by = index(tempB_xts))
t3rollB <- rollapply(taB,
                    width = 3,
                    align = "right",
                    FUN = mean)
ehi_sig_rollB <- t3rollB - t95B
t30rollB <- rollapply(lagtaB,
                     width = 30,
                     align = "right",
                     FUN = mean)
ehi_accl_rollB <- t3rollB - t30rollB
EHFrollB <- ehi_sig_rollB * pmax(1, ehi_accl_rollB)
EHFrollB <- ifelse(EHFrollB < 0, 0, EHFrollB)
```

```{r}
## ----EHFdatB
EHFdatB <- data.frame(EHFB = EHFrollB, 
            Date = index(EHFrollB))
EHFdatB <- EHFdatB |>
  mutate(Decade =
           case_when(
Date > "1950-12-31" & Date < "1960-01-01" ~ "1950s",
Date > "1959-12-31" & Date < "1970-01-01" ~ "1960s",
Date > "1969-12-31" & Date < "1980-01-01" ~ "1970s",
Date > "1979-12-31" & Date < "1990-01-01" ~ "1980s",
Date > "1989-12-31" & Date < "2000-01-01" ~ "1990s",
Date > "1999-12-31" & Date < "2010-01-01" ~ "2000s",
Date > "2009-12-31" & Date < "2020-01-01" ~ "2010s",
Date > "2019-12-31"  ~ "2020s",
           ))
```

```{r}
# EHFggB
ggEHF1B <- EHFdatB |>
  #filter(Decade != "1950s") |>
  #filter(Decade != "2020s") |>
  ggplot(aes(x = Date, y = EHFB)) +
  geom_line() + 
  labs(title = "Bhubaneswar",
       y = "Excess heat factor (Degree Celsius squared)") +
  #facet_wrap(~ Decade, scales = "free_x", ncol = 1) +
  #coord_flip() +
  theme_bw() 
ggplotly(ggEHF1B)
```

> **This is Figure 6 in the paper. Data from the India Meteorological Department.**

***

<font size="3"> In most years, the excess heat factor was low. 

By using rectangular selection, selecting a window across a portion of the x-axis, we can zoom in successively to the large spikes in the EHF, for example, that on the 9th of June, 1967. We can see the values of the EHF for the days and weeks around such spikes. 

</font>

### <font size="3">  **Excess heat factor in Jodhpur from 1951 to 2023** </font>

```{r}
# t95J
tempJ_xts <- tempJ_xts["1951/2023"]
taJ <- as.xts(coredata(tempJ_xts)[,3],
              order.by = index(tempJ_xts))
t95J <- quantile(taJ, 0.95, na.rm = T)
```

```{r}
## ----rollJ
lagtaJ <- as.xts(c(NA, taJ[1:(length(taJ) - 1)]),
                 order.by = index(tempJ_xts))
t3rollJ <- rollapply(taJ,
                     width = 3,
                     align = "right",
                     FUN = mean)
ehi_sig_rollJ <- t3rollJ - t95J
t30rollJ <- rollapply(lagtaJ,
                      width = 30,
                      align = "right",
                      FUN = mean)
ehi_accl_rollJ <- t3rollJ - t30rollJ
EHFrollJ <- ehi_sig_rollJ * pmax(1, ehi_accl_rollJ)
EHFrollJ <- ifelse(EHFrollJ < 0, 0, EHFrollJ)
```

```{r}
# EHFdatJ
EHFdatJ <- data.frame(EHFJ = EHFrollJ, 
                Date = index(EHFrollJ))
EHFdatJ <- EHFdatJ |>
  mutate(Decade =
           case_when(
             Date > "1950-12-31" & Date < "1960-01-01" ~ "1950s",
             Date > "1959-12-31" & Date < "1970-01-01" ~ "1960s",
             Date > "1969-12-31" & Date < "1980-01-01" ~ "1970s",
             Date > "1979-12-31" & Date < "1990-01-01" ~ "1980s",
             Date > "1989-12-31" & Date < "2000-01-01" ~ "1990s",
             Date > "1999-12-31" & Date < "2010-01-01" ~ "2000s",
             Date > "2009-12-31" & Date < "2020-01-01" ~ "2010s",
             Date > "2019-12-31"  ~ "2020s",
           ))
```

```{r}
# EHFggJ
ggEHF1J <- EHFdatJ |>
  #filter(Decade != "1950s") |>
  #filter(Decade != "2020s") |>
  ggplot(aes(x = Date, y = EHFJ)) +
  geom_line() + 
  labs(title = "Jodhpur",
          y = "Excess heat factor (Degree Celsius squared)") +
  #facet_wrap(~ Decade, scales = "free_x", ncol = 1) +
  #coord_flip() +
  theme_bw() 
ggplotly(ggEHF1J)
```

> **This is Figure 7 in the paper. Data from the India Meteorological Department.**

***

<font size="3"> In most years, the excess heat factor was low. 

By using rectangular selection, selecting a window across a portion of the x-axis, we can zoom in successively to the large spikes in the EHF, for example, that on the 22nd of May, 2016. We can see the values of the EHF for the days and weeks around such spikes. 

</font>

### <font size="3">  **Cumulative distribution function of excess heat factor in Bhubaneswar from 1951 to 2023** </font>

```{r}
EHFdatB0 <- EHFdatB |>
  filter(EHFB > 0) |>
  select(EHFB)
Q80B <- quantile(EHFdatB0$EHFB, probs = 0.8)
MaxB <- max(EHFdatB$EHFB)
ggecdfB <- EHFdatB |>
  filter(EHFB > 0) |>
  ggplot(aes(x = EHFB)) +
  stat_ecdf(geom = "step") +
  labs(x = "Excess Heat Factor (Degree Celsius squared)",
       y = "Empirical cumulative distribution function",
       title = "Bhubaneswar") +
  geom_hline(yintercept = 0.8, 
             linetype = "dashed") +
  geom_vline(xintercept = Q80B, 
             linetype = "dashed") +
  #xlim(0, max(EHFdatJ$EHFJ)) +
  theme_bw()
ggplotly(ggecdfB)
```

> **This is part of Figure 8 in the paper. Data from the India Meteorological Department.**

***

<font size="3"> The empirical cumulative distribution function for a location resembles a Generalised Pareto Distribution function. Each location will have a different empirical cumulative distribution function (Nairn et al. 2018).

</font>

### <font size="3">  **Cumulative distribution function of excess heat factor in Jodhpur from 1951 to 2023** </font>

```{r}
EHFdatJ0 <- EHFdatJ |>
  filter(EHFJ > 0) |>
  select(EHFJ)
Q80J <- quantile(EHFdatJ0$EHFJ, probs = 0.8)
MaxJ <- max(EHFdatJ$EHFJ)
ggecdfJ <- EHFdatJ |>
  filter(EHFJ > 0) |>
  ggplot(aes(x = EHFJ)) +
  stat_ecdf(geom = "step") +
  labs(x = "Excess Heat Factor (Degree Celsius squared)",
       y = "Empirical cumulative distribution function",
       title = "Jodhpur") +
  geom_hline(yintercept = 0.8, 
             linetype = "dashed") +
  geom_vline(xintercept = Q80J, 
             linetype = "dashed") +
  #xlim(0, max(EHFdatJ$EHFJ)) +
  theme_bw()
ggplotly(ggecdfJ)
```

> **This is part of Figure 8 in the paper. Data from the India Meteorological Department.**

***

<font size="3"> The empirical cumulative distribution function for a location resembles a Generalised Pareto Distribution function. Each location will have a different empirical cumulative distribution function (Nairn et al. 2018); we can see that for Jodhpur differs from that for Bhubaneswar. 

</font>



### <font size="3">  **Twenty highest excess heat factor days in Bhubaneswar from 1951 to 2023** </font>

```{r}
EHFhighB <- EHFdatB |>
  select("EHFB", "Date") |>
  slice_max(EHFB, n = 20)
  #filter(EHFB > 12)
EHFhighB <- data.frame(EHFB = EHFhighB$EHFB, 
                       Date = EHFhighB$Date)
EHFhighJ <- EHFdatJ |>
  select("EHFJ", "Date") |>
  slice_max(EHFJ, n = 20)
#filter(EHFB > 12)
EHFhighJ <- data.frame(EHFJ = EHFhighJ$EHFJ, 
                       Date = EHFhighJ$Date)
```

```{r}
ggEHFhighB <- EHFhighB |>
  #arrange(desc(EHFB)) |>
  ggplot(aes(x = EHFB, y = factor(Date))) +
  geom_point(stat = "identity", 
             size=2, color="red") +
  labs(title = "Bhubaneswar",
       y = 'Date', 
       x = "Excess heat factor \n (Degree Celsius squared)") +
  theme_bw()
ggplotly(ggEHFhighB)
```

> **This is part of Figure 9 in the paper. Data from the India Meteorological Department.**

***

<font size="3"> In Bhubaneswar the excess heat factor was highest on 9 June 1967 (17.2 degrees celsius squared) and next highest the day before (14.5 degrees celsius squared) (Figure 9). After this, the days with the next highest excess heat factors were 9 June 2012 (13.9 degrees celsius) and 8 June 2012 (11.5 degrees celsius). In Bhubaneswar the consecutive days from 7 June 1967 to 11 June 1967 figure in the twenty highest excess factor days. 

</font>

### <font size="3">  **Twenty highest excess heat factor days in Jodhpur from 1951 to 2023** </font>

```{r}
ggEHFhighJ <- EHFhighJ |>
  #arrange(desc(EHFJ)) |>
  ggplot(aes(x = EHFJ, y = factor(Date))) +
  geom_point(stat = "identity", 
             size=2, color="red") +
  labs(title = "Jodhpur",
       y = 'Date', 
       x = "Excess heat factor \n (Degree Celsius squared)") +
  theme_bw()
ggplotly(ggEHFhighJ)
```

> **This is part of Figure 9 in the paper. Data from the India Meteorological Department.**

***

<font size="3"> The date with the highest excess heat factor in Jodhpur was 28 April 1958, with a value of 25.7 degrees celsius squared (Figure 9). The next day, 29 April 1958, also had a high value of 20 degrees celsius squared. In Jodhpur the consecutive days from 21 May 2016 to 23 May 2016 figure in the twenty highest excess factor days. 

</font>