message("")
message("#########################################")
message("#########################################")
message("#########################################")

# IDS="id5628651 id662424 id13145084 id14620483 id14594730 id550119 id14595225 id12776226 id2398386 id14722068
#     id14573303 id547238 id12777335 id14155848"

# SELECT SENSOR-ID
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
#id <- "id5628651"
id <- "id662424"
#id <- "id12776226"
#id <- "id2398386"
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

args <- commandArgs(trailingOnly = TRUE)
message("ARGUMENTS = ",length(args))
if (length(args)==1) {
  id <- args[1]
}
message("ID = ",id)

# READ DATA
load("oklab.Rdata")
load("edm.Rdata")
load("asopc.Rdata")
message("reading done!")

# COPY DATA TO STANDARD VARIABLE NAMES
pm10 <- eval(parse(text=paste0("oklabaggpm10$",id)))
pm25 <- eval(parse(text=paste0("oklabaggpm25$",id)))
temp <- eval(parse(text=paste0("oklabaggtemp$",id)))
rhum <- eval(parse(text=paste0("oklabaggrhum$",id)))
time <- as.POSIXlt(oklabaggpm10$time,tz="UTC")

# FILTER
pm10[pm10<1|pm10>85] <- NA
pm25[pm25<1|pm25>60] <- NA
temp[temp < -99 | temp > 60] <- NA
rhum[rhum < 1 | rhum > 99] <- NA


t1 <- min(time) # 2023-05-23 16:00:00 UTC
t2 <- max(time) # present
t2 <- as.POSIXct("2023-08-22 06:00:00 CET")

#########################################################

message("#########################################")
message("reading periods.xlsx ...")
library("readxl")
url <- "https://uniaugsburg-my.sharepoint.com/:x:/g/personal/andreas_philipp_uni-a_de/ETgvm1h53WlHowKks387oqgBRptiupQzXp4Pbu2e1rY-4A?e=7GxXEg&download=1"
x <- read_excel("periods.xlsx",na="NA")
#str(x)
x$SensorID <- as.character(x$SensorID)
x$t1excl1 <- as.character(x$t1excl1)
x$t2excl1 <- as.character(x$t2excl1)
x$t1excl2 <- as.character(x$t1excl2)
x$t2excl2 <- as.character(x$t2excl2)
x$t1excl3 <- as.character(x$t1excl3)
x$t2excl3 <- as.character(x$t2excl3)
#str(x)

message("reading period dates ...")
#seq_along(x$SensorID)
for( i in seq_along(x$SensorID) ) {
 # message("i = ",i)
 # message("id = ",x$SensorID[i]) 
 sid <- paste0("id",x$SensorID[i])
 if( sid == id ){
  message("checking ", sid )
  t1train_th <- as.POSIXct(x$t1train_th[i])
  t2train_th <- as.POSIXct(x$t2train_th[i])
  t1train_pm <- as.POSIXct(x$t1train_pm[i])
  t2train_pm <- as.POSIXct(x$t2train_pm[i])
  t1exp <- as.POSIXct(x$t1exp[i])
  t2exp <- as.POSIXct(x$t2exp[i])

  if( ! is.na(x$t1excl1[i]) & ! is.na(x$t2excl1[i]) ){
    t1exclude <- as.POSIXct(x$t1excl1[i])
    t2exclude <- as.POSIXct(x$t2excl1[i])
    temp [ time >= t1exclude & time <= t2exclude ] <- NA
    rhum [ time >= t1exclude & time <= t2exclude ] <- NA
    pm10 [ time >= t1exclude & time <= t2exclude ] <- NA
    pm25 [ time >= t1exclude & time <= t2exclude ] <- NA
  }

  if( ! is.na(x$t1excl2[i]) & ! is.na(x$t2excl2[i]) ){
    t1exclude <- as.POSIXct(x$t1excl2[i])
    t2exclude <- as.POSIXct(x$t2excl2[i])
    temp [ time >= t1exclude & time <= t2exclude ] <- NA
    rhum [ time >= t1exclude & time <= t2exclude ] <- NA
    pm10 [ time >= t1exclude & time <= t2exclude ] <- NA
    pm25 [ time >= t1exclude & time <= t2exclude ] <- NA
  }

  if( ! is.na(x$t1excl3[i]) & ! is.na(x$t2excl3[i]) ){
    t1exclude <- as.POSIXct(x$t1excl3[i])
    t2exclude <- as.POSIXct(x$t2excl3[i])
    temp [ time >= t1exclude & time <= t2exclude ] <- NA
    rhum [ time >= t1exclude & time <= t2exclude ] <- NA
    pm10 [ time >= t1exclude & time <= t2exclude ] <- NA
    pm25 [ time >= t1exclude & time <= t2exclude ] <- NA
  }
} # if sid
} # for i

message("t1train_th = ",t1train_th)
message("t2train_th = ",t2train_th)
message("t1train_pm = ",t1train_pm)
message("t2train_pm = ",t2train_pm)
message("t1exp      = ",t1exp)
message("t2exp      = ",t2exp)

## <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
## DEFINE TRAINING PERIOD FOR TEMP & RHUM:
#t1train_th <- as.POSIXct("2023-05-23 16:00:00 UTC")
#t2train_th <- as.POSIXct("2023-06-26 23:59:59 UTC")
## DEFINE TRAINING PERIOD FOR PM10 & PM25:
#t1train_pm <- as.POSIXct("2023-05-23 16:00:00 UTC")
#t2train_pm <- as.POSIXct("2023-06-26 23:59:59 UTC")
## DEFINE WHICH PERIOD SHOULD BE RECALIBRATED (OBSERVATION PERIOD OF FIELD EXPERIMENT)
#t1exp <- as.POSIXct("2023-06-27 00:00:00 UTC")
#t2exp <- as.POSIXct("2023-07-03 06:00:00 UTC")
## <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#### OPTIONAL: exclude periods for calibration (repeat for more than one exclusion period):
#t1exclude <- as.POSIXct("2023-05-23 16:00:00 UTC")
#t2exclude <- as.POSIXct("2023-05-23 17:00:00 UTC")
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
#temp [ time >= t1exclude & time <= t2exclude ] <- NA
#rhum [ time >= t1exclude & time <= t2exclude ] <- NA
#pm10 [ time >= t1exclude & time <= t2exclude ] <- NA
#pm25 [ time >= t1exclude & time <= t2exclude ] <- NA
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<



# SELECT REFERENCE DATA
#tempref <-asopcaggtemp$id177410813
#tempreftime <-  as.POSIXlt(asopcaggtemp$time + 2*60*60,tz="UTC")
#rhumref <- asopcaggrhum$id177410813
#rhumreftime <-  as.POSIXlt(asopcaggrhum$time + 2*60*60,tz="UTC")

tempref <- edmaggtemp$idOPC004
tempreftime <- as.POSIXlt(edmaggtemp$time,tz="UTC")
rhumref <- edmaggrhum$idOPC004
rhumreftime <- as.POSIXlt(edmaggrhum$time,tz="UTC")


pm10ref <- edmaggpm10$idOPC004
pm10reftime <- as.POSIXlt(edmaggpm10$time,tz="UTC")
pm25ref <- edmaggpm25$idOPC004
pm25reftime <- as.POSIXlt(edmaggpm25$time,tz="UTC")

# FILTER
pm10ref[pm10ref<1|pm10ref>85] <- NA
pm25ref[pm25ref<1|pm25ref>60] <- NA
tempref[tempref < -99 | tempref > 60] <- NA
rhumref[rhumref < 1 | rhumref > 99] <- NA


# PREPARE PLOTS
#tmin <- min(time)
#tmax <- max(time)
tmin <- t1train_pm
tmax <- max(time)

message("tmin =",tmin)
message("tmax =",tmax)
tstart <- min(seq(tmin, tmax, by = "hour"))
tstop <- max(seq(tmin, tmax, by = "hour"))

col1 <- "lightblue"
col2 <- "lavender"


message("tstart =",tstart)
message("tstop =",tstop)


pm10max <- max(c(pm10,pm10ref),na.rm=TRUE)
pm10min <- 0
pm25max <- max(c(pm25,pm25ref),na.rm=TRUE) 
pm25min <- 0


# PLOT PM10
message("plotting raw time series for pm10 ...")
png(paste0(id,"_pm10_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
plot(time,pm10,xlab=paste("time UTC"),ylab="pm10 (µg/m³)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(0,pm10max),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("raw data PM10 oklab 5min-mean (",tmin," - ",tmax," UTC)" ))

rect(xleft=t1train_pm,ybottom=0,xright=t2train_pm,ytop=pm10max,col=col1)
rect(xleft=t1exp,ybottom=0,xright=t2exp,ytop=pm10max,col=col2)

ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600*24) , units="hours")
axis.POSIXct(1, x=time, at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(0,pm10max,by=10),lwd=0.2,col="gray")
lines(pm10reftime,pm10ref,col="black",t="l",lwd=1.5)
lines(time,pm10,xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col="red",t="l",lwd=0.3,lty=1)
legend("topright",c("OPC004"),col=c("black"),lty=1,lwd=0.5)
legend("topleft",id,col="red",lty=1,lwd=0.3)
dev.off()

# PLOT PM25
message("plotting raw time series for pm25 ...")
png(paste0(id,"_pm25_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
plot(time,pm25,xlab=paste("time UTC"),ylab="pm25 (µg/m³)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(0,pm25max),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("raw data PM25 oklab 5min-mean (",tmin," - ",tmax," UTC)" ))

rect(xleft=t1train_pm,ybottom=0,xright=t2train_pm,ytop=pm25max,col=col1)
rect(xleft=t1exp,ybottom=0,xright=t2exp,ytop=pm25max,col=col2)

ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600*24) , units="hours")
axis.POSIXct(1, x=time, at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(0,pm25max,by=10),lwd=0.2,col="gray")
lines(pm25reftime,pm25ref,col="black",t="l",lwd=1.5)
lines(time,pm25,xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col="red",t="l",lwd=0.3,lty=1)
legend("topright",c("OPC004"),col=c("black"),lty=1,lwd=0.5)
legend("topleft",id,col="red",lty=1,lwd=0.3)
dev.off()


tmin <- t1train_th
tmax <- max(time)

# PLOT TEMP
message("plotting raw time series for temp ...")
png(paste0(id,"_temp_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
tempmin <- min(c(temp,tempref),na.rm=TRUE)
tempmax <- max(c(temp,tempref),na.rm=TRUE)
plot(time,temp,xlab=paste("time UTC"),ylab="Temperature (°C)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(tempmin,tempmax),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("raw data Temp oklab 5min-mean (",tmin," - ",tmax," UTC)" ))

message("tstart =",tstart)
message("tstop =",tstop)


rect(xleft=t1train_th,ybottom=tempmin,xright=t2train_th,ytop=tempmax,col=col1)
rect(xleft=t1exp,ybottom=tempmin,xright=t2exp,ytop=tempmax,col=col2)

ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600*24) , units="hours")
axis.POSIXct(1, x=time, at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(tempmin,tempmax,by=5),lwd=0.2,col="gray")
lines(tempreftime,tempref,col="black",t="l",lwd=1.5)
lines(time,temp,xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col="red",t="l",lwd=0.3,lty=1)
legend("topright",c("SHT85"),col=c("black"),lty=1,lwd=0.5)
legend("topleft",id,col="red",lty=1,lwd=0.3)
dev.off()

# PLOT RHUM
message("plotting raw time series for rhum ...")
png(paste0(id,"_rhum_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
rhummin <- min(c(rhum,rhumref),na.rm=TRUE)
rhummax <- max(c(rhum,rhumref),na.rm=TRUE)
plot(time,rhum,xlab=paste("time UTC"),ylab="Relative Humidity (%)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(rhummin,rhummax),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("raw data Rhum oklab 5min-mean (",tmin," - ",tmax," UTC)" ))

rect(xleft=t1train_th,ybottom=rhummin,xright=t2train_th,ytop=rhummax,col=col1)
rect(xleft=t1exp,ybottom=rhummin,xright=t2exp,ytop=rhummax,col=col2)

ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600*24) , units="hours")
axis.POSIXct(1, x=time, at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(rhummin,rhummax,by=5),lwd=0.2,col="gray")
lines(rhumreftime,rhumref,col="black",t="l",lwd=1.5)
lines(time,rhum,xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col="red",t="l",lwd=0.3,lty=1)
legend("topright",c("SHT85"),col=c("black"),lty=1,lwd=0.5)
legend("topleft",id,col="red",lty=1,lwd=0.3)
dev.off()

# CREATE MERGED DATA FRAME "okm" OF SAME LENGTH (including NAs)
ok <- data.frame(time=time,pm10=pm10,pm25=pm25,temp=temp,rhum=rhum)
as <- data.frame(time=tempreftime,tempref=tempref,rhumref=rhumref)
edm <- data.frame(time=pm10reftime,pm10ref=pm10ref,pm25ref=pm25ref)
okm <- merge(ok,edm,by="time",all.x=TRUE)
okm <- merge(okm,as,by="time",all.x=TRUE)


# CREATE RAW DATA SCATTER PLOTS
 # ACCURACY METRICS & SCATTER PLOT PM25
 mae <- mean(abs(okm$pm10-okm$pm10ref),na.rm=T)
 bias <- mean(okm$pm10-okm$pm10ref,na.rm=T)
 rmse <- sqrt(sum( (okm$pm10-okm$pm10ref)**2,na.rm=T)/length(okm$pm10[!is.na(okm$pm10)]))
 r2 <- cor(okm$pm10,okm$pm10ref,method="spearman")
message("plotting raw scatter plot for pm10 ...")
png(paste0(id,"_pm10_scatter_raw.png"),width=800, height=800, units="px", res=300, pointsize=3)
plot(okm$pm10ref,okm$pm10,xlim=c(0,pm10max),ylim=c(0,pm10max),ylab="PM10 SDS011 (µg/m³)",xlab="PM10 EDM164 (µg/m³)",main=paste("PM10 raw data ",id," versus reference"))
abline(0,1)
dev.off()

message("plotting raw scatter plot for pm25 ...")
png(paste0(id,"_pm25_scatter_raw.png"),width=800, height=800, units="px", res=300, pointsize=3)
plot(okm$pm25ref,okm$pm25,xlim=c(0,pm25max),ylim=c(0,pm25max),ylab="PM25 SDS011 (µg/m³)",xlab="PM25 EDM164 (µg/m³)",main=paste("PM25 raw data ",id," versus reference"))
abline(0,1)
dev.off()


message("plotting raw scatter plot for temp ...")
png(paste0(id,"_temp_scatter_raw.png"),width=800, height=800, units="px", res=300, pointsize=3)
plot(okm$tempref,okm$temp,xlim=c(tempmin,tempmax),ylim=c(tempmin,tempmax),ylab="TEMP BME128 (°C)",xlab="TEMP SHT85 (°C)",main=paste("Temp raw data ",id," versus reference"))
abline(0,1)
dev.off()

message("plotting raw scatter plot for rhum ...")
png(paste0(id,"_rhum_scatter_raw.png"),width=800, height=800, units="px", res=300, pointsize=3)
plot(okm$rhumref,okm$rhum,xlim=c(rhummin,rhummax),ylim=c(rhummin,rhummax),ylab="RHUM BME128 (°C)",xlab="RHUM SHT85 (°C)",main=paste("Rhum raw data ",id," versus reference"))
abline(0,1)
dev.off()