library(psychrolib)
SetUnitSystem("SI")


# Step 1: Create a dataframe with hourly temperature data for a full year

# Let's assume the dataframe is named "temp_data" and has two columns: "datetime" (with hourly timestamps) and "temperature" (with hourly temperature readings).
df <- read.csv("data2020.csv",col.names=c("DateTime", "DB", "RH","DP","P","WIND","WIND-DIR","SUN","RAIN","START"))

df=df[complete.cases(df), ]

date_format <- "%d/%m/%Y %H:%M"
df$DateTime=as.POSIXct(df$DateTime, format = date_format)
df$DB=(df$DB/10)
df$RH=(df$RH/100)
df$P=(df$P*100)
df$WB=round(GetTWetBulbFromRelHum(df$DB, df$RH, df$P),1)

df=df[format(df$DateTime, "%m") %in% c("06", "07", "08", "09"),]
total_hours=nrow(df)

dfdb=aggregate(WB ~ DB, data = df, FUN = mean)
dfdb$count <- sapply(dfdb$DB, function(x) sum(df$DB > x))
dfdb$pct_exceeded=(dfdb$count/total_hours)*100
dfdb <- subset(dfdb, pct_exceeded <= 0.02)
nrow(dfdb)
head(dfdb[order(dfdb$pct_exceeded), ],n=20)

#indexdb=which.min(abs(dfdb$pct_exceeded - 1))
#print(dfdb[indexdb,'DB'])


dfwb=aggregate(DB ~ WB, data = df, FUN = mean)
dfwb$count <- sapply(dfwb$WB, function(x) sum(df$WB > x))
dfwb$pct_exceeded=(dfwb$count/total_hours)*100
dfwb <- subset(dfwb, pct_exceeded <= 0.02)
nrow(dfwb)
head(dfwb[order(dfwb$pct_exceeded), ],n=50)
#indexwb=which.min(abs(dfwb$pct_exceeded - 1))
#print(dfwb[indexwb,'WB'])

library(ggplot2)
library(ggpsychro)


p = ggpsychro(tdb_lim = c(-15, 40)) +
    geom_grid_relhum(alpha = 1.0, label.alpha = 1.0, label.size = 6, label.fontface = 2) +
    scale_relhum(minor_breaks = NULL) +
    geom_grid_wetbulb(size = 1.0, color = "black", alpha = 1.0, label_loc = NA) +
    scale_wetbulb(breaks = seq(25, 30, by = 5), minor_breaks = NULL) +
    geom_grid_vappres(label.size = 5) +
    scale_vappres(breaks = seq(6000, 7000, by = 500), limits = c(6000, 7000)) +
    geom_grid_specvol() +
    scale_specvol(labels = NULL) +
    geom_point(aes(dfdb$DB, wetbulb=dfdb$WB),stat = "wetbulb",color="red")+
    geom_vline(xintercept = min(dfdb$DB), color = "red" )+
    geom_line(aes(x = 25:min(dfdb$DB), wetbulb = min(dfwb$WB)), stat = "wetbulb",color="blue")+
    geom_point(aes(dfwb$DB, wetbulb=dfwb$WB),stat = "wetbulb",color="blue")+
    geom_grid_wetbulb()+
    scale_wetbulb(breaks = seq(0, 30, by = 2), minor_breaks = NULL) 

ggsave("chartcombop.pdf", plot = p, width = 6, height = 4)


# Step 8: Create a list of unique combinations of "RH_pct_exceeded" and "pct_exceeded"
#combo_list <- unique(temp_pct[, c("RH_pct_exceeded", "pct_exceeded")])