---
title: "Consolidated Code"
author: "EPatton"
date: "2024-04-09"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
#install.packages("/Users/erikpatton/Downloads/heatmetrics_1.0.tar",repos=NULL,type="source")
#install.packages("ecmwfr")
#library(heatmetrics)
library(wbgt)
library(ncdf4)
library(keyring)
#library(foreach)
library(ecmwfr)
library(lubridate)
library(humidity)
library(tidyverse)
```
STEP #00: Make a series for files to deposit historical data using the Terminal in macOS. 

This command:
mkdir -p "/Users/erikpatton/Library/Mobile Documents/com~apple~CloudDocs/R Files/Consolidated code"/{1994..2023}


Step #0: Set up ERA5 data download
##Note the user will have to replace the user name and key with their own API credentials
```{r, eval=FALSE}
wf_set_key(user ="XXX", key = "XXX", service = "cds")
user <- "XXX"
```
Download historical ERA 5 data for study location (Duke University)
```{r, eval=FALSE}
tzone = "America/New_York"

yr <- c(2019:2021)
#ph <- paste0("/Users/erikpatton/Library/Mobile Documents/com~apple~CloudDocs/R Files/Consolidated code/Durham/", yr)
area <- "36/-79/35.99/-78.99"

vr <- c("2m_temperature", "surface_solar_radiation_downwards", 
               "10m_u_component_of_wind", "10m_v_component_of_wind", 
               "surface_pressure","2m_dewpoint_temperature")
mn <- c(sprintf("0%d",1:9),10:12)
dy <- c(sprintf("0%d",1:9), 10:31)
tm <- c(paste0("0",0:9,":00"),paste0(10:23,":00"))

  for (i in 1:length(yr)) {
  for (j in 1:length(mn)) {

    ph <- paste0("/Users/erikpatton/Library/Mobile Documents/com~apple~CloudDocs/R Files/Consolidated code/Durham/", yr[i])
    request <- list(
      "dataset_short_name" = "reanalysis-era5-land",
      "product_type" = "reanalysis",
      "variable" = vr,
      "year" = yr[i],
      "month" = mn[j],
      "day" = dy,
      "time" = tm,
      "area" = area,
      "format" = "netcdf",
      "target" = paste0("era5_",yr[i],"_",mn[j],".nc")
    )
    
    file <- wf_request(
      user     = "174423",   # user ID (for authentication)
      request  = request,  # the request
      transfer = TRUE,     # download the file
      path     = ph       # store data in current working directory
    )
    
  }
}

```

Step #1: Create a data frame of the required weather variables (temperature, relative humidity, solar radiation, surface pressure, and windspeed)
#Step #1
```{r}
#create blank data frames to read data onto
Durham.temp <- NULL
Durham.UWind <- NULL
Durham.VWind <- NULL
Durham.Solar <- NULL
Durham.Pressure <- NULL
Durham.dewpoint <- NULL
Durham.time <- NULL

#select the years and time zone for the study location
yr <- c(1994:1999,2018)
mn <- c(sprintf("0%d",1:9),10:12)
tzone = "America/New_York"
#Latitude and Longitude for the study location
Lat=36.0
Lon=-79.0


for (k in 1:length(yr)){
for (i in mn){
        read_in <- nc_open(paste0("/Users/erikpatton/Library/Mobile Documents/com~apple~CloudDocs/R Files/Consolidated code/Durham/",yr[k],"/era5_",yr[k],"_",i,".nc"))
        
        assign(("new_temp_row"),ncvar_get(read_in,"t2m"))
        Durham.temp <- c(Durham.temp,as.numeric(new_temp_row),sep="")

        assign(("new_UWind_row"),ncvar_get(read_in,"u10"))
        Durham.UWind <- c(Durham.VWind,as.numeric(new_UWind_row),sep="")
        
        assign(("new_VWind_row"),ncvar_get(read_in,"v10"))
        Durham.VWind <- c(Durham.VWind,as.numeric(new_VWind_row),sep="")

        assign(("new_Solar_row"),ncvar_get(read_in,"ssrd"))
        Durham.Solar <- c(Durham.Solar,as.numeric(new_Solar_row),sep="")

        assign(("new_Pressure_row"),ncvar_get(read_in,"sp"))
        Durham.Pressure <- c(Durham.Pressure,as.numeric(new_Pressure_row),sep="")
        
        assign(("new_dewpoint_row"),ncvar_get(read_in,"d2m"))
        Durham.dewpoint <- c(Durham.dewpoint,as.numeric(new_dewpoint_row),sep="")
        
        assign(("new_year_row"),ncvar_get(read_in,"time"))
        Durham.time <- c(Durham.time,as.numeric(new_year_row),sep="")
}}

Durham.temp <- as.data.frame(Durham.temp)

#create a datetime column in local time
Durham.time <- as.numeric(Durham.time)
seconds_since_origin <- Durham.time * 3600
origin_date <- ymd_hms("1900-01-01 00:00:00.0", tz = "UTC")
result_datetime <- origin_date + seconds(seconds_since_origin)
local_datetime <- with_tz(result_datetime, tzone = tzone)

Durham.UWind <- as.data.frame(Durham.UWind)

Durham.VWind <- as.data.frame(Durham.VWind)

Durham.Solar <- as.data.frame(Durham.Solar)

Durham.Pressure <- as.data.frame(Durham.Pressure)

Durham.dewpoint <- as.data.frame(Durham.dewpoint)

Durham.ERA5 <- cbind(Durham.temp,Durham.UWind,Durham.VWind,Durham.Solar,Durham.Pressure,Durham.dewpoint)


##covert to numeric 
Durham.ERA5[] <- lapply(Durham.ERA5, function(x) as.numeric(as.character(x)))  

#add in the datetime column
Durham.ERA5$local_datetime <- local_datetime

##Remove empty rows between months
Durham.ERA5 <- na.omit(Durham.ERA5)

##Covert Joules to Watts
Durham.ERA5$Durham.Solar <- Durham.ERA5$Durham.Solar / (60 * 24 * 24)
Durham.ERA5$Durham.WindSpeed <- sqrt(Durham.ERA5$Durham.UWind^2 + Durham.ERA5$Durham.VWind^2)

##convert dewpoint to relative humidity
Durham.ERA5$relative_humidity <- humidity::RH(Durham.ERA5$Durham.temp,Durham.ERA5$Durham.dewpoint,isK=TRUE)


##Create hour and date column
Durham.ERA5$hour_local <- hour(Durham.ERA5$local_datetime)
Durham.ERA5$year <- year(Durham.ERA5$local_datetime)
Durham.ERA5$date <- date(Durham.ERA5$local_datetime)
Durham.ERA5$julian <- yday(Durham.ERA5$local_datetime)

#remove the U and V components column
Durham.ERA5 <- Durham.ERA5%>%
  select(-c(Durham.UWind,Durham.VWind,Durham.dewpoint))%>%
  mutate(Durham.temp=Durham.temp-273.15)

Durham.ERA5$month <- month(Durham.ERA5$local_datetime)
Durham.ERA5$day <- day(Durham.ERA5$local_datetime)

Durham.ERA5$Durham.Pressure <- Durham.ERA5$Durham.Pressure / 100
#View(Durham.ERA5)
```

Step 2: Create a WBGT column
#Step #2
```{r}
#clean up minimum windspeed to reflect 0.62m/s at 10m (this will get adjusted to ~0.5m/s at 2m elevation). This is required because the Liljegren method will give bad values at very low windspeeds, and natural processes create at least 0.5m/s of air movement over skin surface
Durham.ERA5$Durham.WindSpeed <- ifelse(Durham.ERA5$Durham.WindSpeed < 0.62, 0.62, Durham.ERA5$Durham.WindSpeed)

#Create a column for the difference between local time and GMT. This is needed as an input to the "gmt" variable in the wbgt.
Durham.ERA5$local_datetime_posix <- with(Durham.ERA5, as.POSIXct(local_datetime, tz = tzone))
# Create a new POSIXct object with GMT timezone
Durham.ERA5$gmt_datetime_posix <- with(Durham.ERA5, as.POSIXct(local_datetime, tz = "GMT"))

# Calculate the difference in hours and create a new column
Durham.ERA5$time_difference <- with(Durham.ERA5, hour(local_datetime_posix) - hour(gmt_datetime_posix))
#replace the first and last 5 values with corrected values (-5)
Durham.ERA5$time_difference[Durham.ERA5$time_difference > 0] <- -5

library(wbgt)

##set variables required
gmt <- Durham.ERA5$time_difference #four hour difference between Ft Moore time and GMT
avg <- rep(60,nrow(Durham.ERA5))
Lat <- rep(Lat,nrow(Durham.ERA5))
Lon <- rep(Lon,nrow(Durham.ERA5))
zspeed <- rep(10,nrow(Durham.ERA5)) #this is the height of the wind speed measurement in meters
#dT is the difference in temperature between measured and reference height. set to 0.
dT <- rep(0,nrow(Durham.ERA5)) #numeric vector of differential between zspeed and reference height
urban <- rep(0,nrow(Durham.ERA5)) #1 for urban, 0 for rural settings
minute <- rep(1,nrow(Durham.ERA5))

##Calculate the WBGT

WBGT_Liljegren <- wbgt(Durham.ERA5$year,
                                        Durham.ERA5$month,
                                        Durham.ERA5$day,
                                        Durham.ERA5$hour_local,
                                        minute,
                                        gmt,                          
                                        avg,
                                        Lat,
                                        Lon,
                                        Durham.ERA5$Durham.Solar,
                                        Durham.ERA5$Durham.Pressure,
                                        Durham.ERA5$Durham.temp,
                                        Durham.ERA5$relative_humidity,
                                        Durham.ERA5$Durham.WindSpeed,
                                        zspeed,
                                        dT,
                                        urban)

Durham.ERA5$WBGT_Liljegren <- WBGT_Liljegren$Twbg
```

Step #3: Prepare the 30-year average variables to use for mapping GCM data to sub-daily WBGT
#Step #3
```{r}
#Create a scaling factor to go from mean daily solar radiation to maximum daily solar radiation and daily mean WindSpeed to max and min WindSpeed

#find the mean values
Durham.ERA5.MeantoMax.Scaling.mean <- Durham.ERA5%>%
  group_by(julian)%>%
  summarize(Durham.Solar.Mean=mean(Durham.Solar),Durham.WindSpeed.Mean=mean(Durham.WindSpeed))
 
#find the mean values at maximum WBGT
Durham.ERA5.MeantoMax.Scaling.max <- Durham.ERA5%>%
  group_by(date)%>%
  slice_max(WBGT_Liljegren)%>%
  group_by(julian)%>%
  summarize(Durham.Solar.Max=mean(Durham.Solar),Durham.WindSpeed.Max=mean(Durham.WindSpeed))%>%
  select(-c(julian))

#find the mean values at min WBGT
Durham.ERA5.MeantoMax.Scaling.min <- Durham.ERA5%>%
  group_by(date)%>%
  slice_min(WBGT_Liljegren)%>%
  group_by(julian)%>%
  summarize(Durham.Solar.Min=mean(Durham.Solar),Durham.WindSpeed.Min=mean(Durham.WindSpeed))%>%
  select(-c(julian))

#Create a final data frame with 365 days
Durham.ERA5.MeantoMax.Scaling <- cbind(Durham.ERA5.MeantoMax.Scaling.mean,Durham.ERA5.MeantoMax.Scaling.max,Durham.ERA5.MeantoMax.Scaling.min)

Durham.ERA5.MeantoMax.Scaling <- Durham.ERA5.MeantoMax.Scaling%>%
  mutate(Solar.MeanToMax.Scaling = Durham.Solar.Max - Durham.Solar.Mean)%>%
  mutate(WindSpeed.MeanToMax.Scaling = Durham.WindSpeed.Max - Durham.WindSpeed.Mean)%>%
  mutate(WindSpeed.MeanToMin.Scaling = Durham.WindSpeed.Min - Durham.WindSpeed.Mean)%>%
  select(julian,Solar.MeanToMax.Scaling,WindSpeed.MeanToMax.Scaling,WindSpeed.MeanToMin.Scaling)

##Clean up the environment
remove(Durham.ERA5.MeantoMax.Scaling.max,Durham.ERA5.MeantoMax.Scaling.mean,Durham.ERA5.MeantoMax.Scaling.min,Durham.dewpoint,Durham.Pressure,Durham.Solar,Durham.temp,Durham.UWind,Durham.VWind,read_in)

#Create linear models to predict surface pressure at max, mean, and min WBGT
library(EnvStats)

##Create a year of data at mean daily WBGT
LM.dataframe.meanWBGT <- Durham.ERA5%>%
  select(Durham.temp:Durham.Pressure,Durham.WindSpeed,relative_humidity,WBGT_Liljegren,julian)%>%
  group_by(julian)%>%
  mutate(mean_WBGT=mean(WBGT_Liljegren))%>%
  mutate(difference=abs(WBGT_Liljegren-mean_WBGT))%>%
    group_by(julian)%>%
  slice(which.min(difference))

##create the linear model to predict surface pressure
lm.press_at_mean.WBGT <- lm(data=LM.dataframe.meanWBGT,Durham.Pressure~Durham.temp+Durham.Solar+Durham.WindSpeed+relative_humidity)
summary(lm.press_at_mean.WBGT)  

##refine the linear model dataframe to test predicted values vs observed. The variables that go into the linear model need to be the same name as the variables used in to create the linear model
LM.dataframe.meanWBGT.test <- LM.dataframe.meanWBGT%>%
  select(Durham.temp,Durham.Solar,relative_humidity,Durham.WindSpeed)
##create the predicted pressure values
LM.dataframe.meanWBGT$test <- as.numeric(predict(lm.press_at_mean.WBGT,LM.dataframe.meanWBGT))

#cor.test(LM.dataframe.meanWBGT$Durham.Pressure,LM.dataframe.meanWBGT$test)

##Create a year of data at max daily WBGT
LM.dataframe.maxWBGT <- Durham.ERA5%>%
  select(Durham.temp:Durham.Pressure,Durham.WindSpeed,relative_humidity,WBGT_Liljegren,julian)%>%
    group_by(julian)%>%
  slice_max(WBGT_Liljegren)

##create the linear model to predict surface pressure
lm.press_at_max.WBGT <- lm(data=LM.dataframe.maxWBGT,Durham.Pressure~Durham.temp+Durham.Solar+Durham.WindSpeed+relative_humidity)
summary(lm.press_at_max.WBGT)  

##refine the linear model dataframe to test predicted values vs observed. The variables that go into the linear model need to be the same name as the variables used in to create the linear model
LM.dataframe.maxWBGT.test <- LM.dataframe.maxWBGT%>%
  select(Durham.temp,Durham.Solar,relative_humidity,Durham.WindSpeed)
##create the predicted pressure values
LM.dataframe.maxWBGT$test <- as.numeric(predict(lm.press_at_max.WBGT,LM.dataframe.maxWBGT))

#cor.test(LM.dataframe.maxWBGT$Durham.Pressure,LM.dataframe.maxWBGT$test)


##Create a year of data at min daily WBGT
LM.dataframe.minWBGT <- Durham.ERA5%>%
  select(Durham.temp:Durham.Pressure,Durham.WindSpeed,relative_humidity,WBGT_Liljegren,julian)%>%
    group_by(julian)%>%
  slice_min(WBGT_Liljegren)

##create the linear model to predict surface pressure
lm.press_at_min.WBGT <- lm(data=LM.dataframe.minWBGT,Durham.Pressure~Durham.temp+Durham.Solar+Durham.WindSpeed+relative_humidity)
summary(lm.press_at_min.WBGT)  

##refine the linear model dataframe to test predicted values vs observed. The variables that go into the linear model need to be the same name as the variables used in to create the linear model
LM.dataframe.minWBGT.test <- LM.dataframe.minWBGT%>%
  select(Durham.temp,Durham.Solar,relative_humidity,Durham.WindSpeed)
##create the predicted pressure values
LM.dataframe.minWBGT$test <- as.numeric(predict(lm.press_at_min.WBGT,LM.dataframe.minWBGT))

#cor.test(LM.dataframe.minWBGT$Durham.Pressure,LM.dataframe.minWBGT$test)
```

Step #4: Create the historical 30-year hourly WBGT average
#Step #4
```{r}
Durham.ERA5.Hourly.Avg <- Durham.ERA5%>%
  group_by(month,day,hour_local)%>%
  summarise_at(c("WBGT_Liljegren"), mean,na.rm=TRUE)%>%
  mutate(datetime=make_datetime(year=2024,month=month,day=day,hour=hour_local))

#display the average WBGT curve over the year
ggplot(Durham.ERA5.Hourly.Avg)+geom_line(aes(x=datetime,y=WBGT_Liljegren))
```