library(psychrolib)
SetUnitSystem("SI")
#library(lubridate)


# 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("data.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=subset(df, year(DateTime) == 2019)
#df=subset(df, format(as.Date(DateTime), "%Y") == "2019")
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)

percent=((100-0.1)/100)
#print(quantile(df$DB, percent))
#print(quantile(df$WB, percent))

quantile(df$DB,((100-0.4)/100),names = FALSE)
quantile(df$WB,((100-0.4)/100),names = FALSE)


results = data.frame("year" = integer(),"Aexdb" = numeric(),"Aexwb"=numeric(),"Bexdb" = numeric(),"Bexwb"=numeric())

for (year in 1999:2022) {
  # create a subset of the data for the current year
  year_data <- subset(df, format(DateTime, "%Y") == as.character(year)) 
  Aexdb=quantile(year_data$DB,((100-0.4)/100),names = FALSE)
  Aexwb=quantile(year_data$WB,((100-0.4)/100),names = FALSE)
  Bexdb=quantile(year_data$DB,((100-1.0)/100),names = FALSE)
  Bexwb=quantile(year_data$WB,((100-1.0)/100),names = FALSE)

#  print(paste("In",year,"the WB", percent,"excedance is",quantile(year_data$WB,percent,type=1)))
  results = rbind(results,data.frame("year"=year,"Aexdb"=Aexdb,"Aexwb"=Aexwb,"Bexdb"=Bexdb,"Bexwb"=Bexwb)) 
}
print(results)
print(max(results$Aexdb))
print(max(results$Aexwb))

#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")])