#library(randomcoloR)
library(data.table)
Sys.setenv(TZ = "UTC")
options(digits=2)

###################################################################
# SELECT TIME PERIOD
# seconds for  2 days
daystep <- 60*60*24
step <- daystep*2
args <- commandArgs(trailingOnly = TRUE)
message("ARGUMENTS = ",length(args))

# SELECT ID
slist <- c(as.character(args[1]))
#slist <- c("14462360")

#if (length(args)==0) {
  t2 <- as.POSIXlt(Sys.time(),tz="UTC")
  t2 <- t2 + 60*60
  #t2 <- Sys.time()
  t1 <- t2 - step
  #"args:" length(args)
#} else {
#  t1 <- as.POSIXlt(args[1])
#  t2 <- as.POSIXlt(args[2])
#}

#*********************************************************
# SET TIME PERIOD HERE (if omitted: last two days)
t1 <- as.POSIXct("2022-05-17 00:00:00 CET")
t2 <- as.POSIXct("2022-07-25 23:59:59 CET")

# SET AGGREGATION INTERVAL:
# aggregation <- "60 sec"
aggregation <- "5 min"
# SET LIST OF OKLAB SENSOR IDs:
#slist <- c("14155848","547446","2398386","14722068","2519718",
#           "5682602","12776226","547238","550119","14595225",
#           "5682100","660436","14462360",
#           "14573303","5630365")
#slist <- c("14155848")
#slist <- c("547446")
#slist <- c("2398386")
#slist <- c("14722068")
#slist <- c("2519718")
#slist <- c("5682602")
#slist <- c("12776226")
#slist <- c("547238")
#slist <- c("550119")
#slist <- c("14595225")
#slist <- c("5682100")
#slist <- c("660436")
#slist <- c("14462360")
#slist <- c("14573303")
#slist <- c("5630365")

#*********************************************************



###################################################################
###################################################################
# GET MATCHING OKLAB FILE NAMES FOR SELECTED DAYS
t1 <- as.POSIXlt(format(t1,"%Y-%m-%d %H:%M:00"))
t2 <- as.POSIXlt(format(t2,"%Y-%m-%d %H:%M:00"))
#t2 <- t2 + 60*60
message("PERIOD = ",t1," ",t2)
runningdate <- seq.Date( from=as.Date(t1), to=as.Date(t2), by="day")

slist <- sort(slist)

flist <- NULL
for (j in seq_along(slist)) {
for (i in seq_along(runningdate)) {
  name1 <- runningdate[i]
  #name1 <- "2021-11-03"
  #name1 <- format(t1,"%Y-%m-%d")
  name2 <- slist[j]
  #name3 <- runningdate[3]
  name3 <- "txt"
  url <- paste("https://saqn.geo.uni-augsburg.de/oklab.cgi?name1=",name1,"&name2=",name2,"&name3=",name3,sep="")
  #message("URL = ",url)
  flist_tmp <- try(fread(url,stringsAsFactors=F,header=F))
  if(inherits(flist_tmp, "try-error")){
    message("no files")
  }else{
    #message("FILE LIST = ",flist_tmp)
    flist <- rbind(flist,flist_tmp)
  }
}}

###################################################################
# READ MATCHING FILES DATA IN OBJECT oklab
header <- c("datetime","pm10","pm25","temp","rhum","v1","v2","v3","v4")

oklab <- data.frame(
  id=character(),datetime=character(),pm10=double(),pm25=double(),temp=double(),rhum=double(),
  v1=numeric(),v2=numeric(),v3=numeric(),v4=numeric(),
  stringsAsFactors=FALSE)

id <- ""
for (i in seq_along(flist$V1)) {
  id <- strsplit(flist$V1[i],"_",fixed=T)[[1]][2]
  id <- strsplit(id,".",fixed=T)[[1]][1]
  fname <- paste("https://saqn.geo.uni-augsburg.de/oklab/",flist$V1[i],sep="")
  message(id," FILE ADRESS = ",fname)
    x <- try(fread(fname,header=FALSE,stringsAsFactors=FALSE,sep=",",dec=".",skip=0,col.names=header,fill=T))
    x$temp[x$temp < -99 | x$temp > 99 ] <- NA
    x$rhum[x$rhum < -99 ] <- NA
    x$rhum[x$rhum > 99 ] <- NA
    x$pm10[x$pm10 < -99 ] <- NA
    x$pm25[x$pm25 < -99 ] <- NA
    if( typeof(x) == "character"){
      next
    }
    if( is.data.frame(x) == F ){
      next
    }
    if (nrow(x)<1) {
       next
    }
    x$id <- id   
    oklab_tmp <- data.frame(id=x$id, datetime=x$datetime, pm10=x$pm10, pm25=x$pm25, temp=x$temp, rhum=x$rhum )
    oklab <- rbind(oklab,oklab_tmp)
}
oklab$time <- as.POSIXlt(oklab$datetime,format="%Y-%m-%d %H:%M:%S",tz="UTC")
oklab[as.character(oklab)=="nan"] <- NA



# --------------------------------------------------------------------------
# DELETE WRONG DATA

tx1 <- as.POSIXlt("2022-07-05 15:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
tx2 <- as.POSIXlt("2022-07-05 20:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
oklab$pm10[oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$pm25[oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$temp[oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$rhum[oklab$time>=tx1 & oklab$time<=tx2] <- NA
#tx1 <- as.POSIXlt("2022-07-12 15:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
#tx2 <- as.POSIXlt("2022-07-13 20:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
#oklab$pm10[oklab$time>=tx1 & oklab$time<=tx2] <- NA
#oklab$pm25[oklab$time>=tx1 & oklab$time<=tx2] <- NA
#oklab$temp[oklab$time>=tx1 & oklab$time<=tx2] <- NA
#oklab$rhum[oklab$time>=tx1 & oklab$time<=tx2] <- NA

for (i in seq_along(slist)) {

if (slist[i] == "12776226"){
tx1 <- as.POSIXlt("2022-05-25 04:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
tx2 <- as.POSIXlt("2022-05-31 16:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
oklab$pm10[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$pm25[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$temp[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$rhum[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
tx1 <- as.POSIXlt("2022-06-04 20:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
tx2 <- as.POSIXlt("2022-06-06 11:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
oklab$pm10[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$pm25[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$temp[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$rhum[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
}
if (slist[i] == "14722068"){
tx1 <- as.POSIXlt("2022-05-17 00:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
tx2 <- as.POSIXlt("2022-06-14 12:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
oklab$pm10[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$pm25[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$temp[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$rhum[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
}
if (slist[i] == "2398386"){
tx1 <- as.POSIXlt("2022-05-17 00:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
tx2 <- as.POSIXlt("2022-06-02 00:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
oklab$pm10[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$pm25[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$temp[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
oklab$rhum[oklab$id==id & oklab$time>=tx1 & oklab$time<=tx2] <- NA
}
}



###################################################################
###################################################################
# PM REFERENCE DATA: GET MATCHING FILE NAMES FOR SELECTED DAYS
flist <- NULL
for (i in seq_along(runningdate)) {
  name1 <- runningdate[i]
  name2 <- "OPC-004"
  name3 <- "txt"
  url <- paste("https://saqn.geo.uni-augsburg.de/grimm.cgi?name1=",name1,"&name2=",name2,"&name3=",name3,sep="")
  #message("URL = ",url)
  flist_tmp <- try(read.table(url,stringsAsFactors=F))
  if(inherits(flist_tmp, "try-error")){
    message("no files")
  }else{
    flist <- rbind(flist,flist_tmp)
  }
}
edmheader <- c("datetime","UTC","longitude","latitude","height","PM10","PM4","PM25","PM1","um0_25","um0_28","um0_30",
            "um0_35","um0_40","um0_45","um0_50","um0_58","um0_65","um0_70","um0_80","um1_00","um1_30","um1_60","um2_00",
            "um2_50","um3_00","um3_50","um4_00","um5_00","um6_50","um7_50","um8_50","um10_00","um12_50","um15_00","um17_50",
            "um20_00","um25_00","um30_00","um32_00","Tout","rH","pressure","serialno","softwareversion")
opcid <- ""
for (i in seq_along(flist$V1)) {
  opcid <- strsplit(flist$V1[i],"/",fixed=T)[[1]][2]
  opcid <- gsub("-", "", opcid)
  fname <- paste("https://saqn.geo.uni-augsburg.de/grimm/",flist$V1[i],sep="")
  message(opcid," FILE ADRESS = ",fname)
  x <- try(read.table(fname,header=F,stringsAsFactors=F,comment.char="*",sep=";",dec=",",skip=0,col.names=edmheader,fill=T))
    if( typeof(x) == "character"){
      next
    }
    if( is.data.frame(x) == F ){
      next
    }
    if (nrow(x)<1) {
       next
    }
  x$opcid <- opcid
    oklab_tmp <- data.frame(id=x$opcid, datetime=as.character(x$UTC), pm10=x$PM10, pm25=x$PM25, temp=x$Tout, rhum=x$rH )
    # FILE ERRORS:
    oklab_tmp$datetime[grepl(",",x$UTC)] <- NA
    oklab_tmp$time <- as.POSIXlt(oklab_tmp$datetime,format="%Y-%m-%d %H:%M:%S",tz="UTC")
    oklab <- rbind(oklab,oklab_tmp)
}
terr <- as.POSIXct("2022-05-17 08:15:00 CET")
oklab$pm10[ oklab$id=="OPC004" & oklab$datetime < terr] <- NA
oklab$pm25[ oklab$id=="OPC004" & oklab$datetime < terr] <- NA



###################################################################
###################################################################
# TEMP REFERENCE DATA: GET MATCHING FILE NAMES FOR SELECTED DAYS
flist <- NULL
for (i in seq_along(runningdate)) {
  name1 <- runningdate[i]
  name2 <- "177410813"
  name3 <- ""
  url <- paste("https://saqn.geo.uni-augsburg.de/asopc.cgi?name1=",name1,"&name2=",name2,"&name3=",name3,sep="")
  #message("URL = ",url)
  flist_tmp <- try(read.table(url,stringsAsFactors=F))
  if(inherits(flist_tmp, "try-error")){
    message("no files")
  }else{
    flist <- rbind(flist,flist_tmp)
  }
}
for (i in seq_along(flist$V1)) {
  opcid <- "177410813"
  fname <- paste("https://saqn.geo.uni-augsburg.de/asopc/",flist$V1[i],sep="")
  message(opcid," FILE ADRESS = ",fname)
  x <- try(fread(fname,stringsAsFactors=F,header=T))
    if( typeof(x) == "character"){
      next
    }
    if( is.data.frame(x) == F ){
      next
    }
    if (nrow(x)<1) {
       next
    }
  x$opcid <- opcid
    oklab_tmp <- data.frame(id=x$opcid, datetime=as.character(x$time_UTC), pm10=as.numeric(x$PM_10), pm25=as.numeric(x$PM_2.5), 
                            temp=as.numeric(x$SHT_Temp), rhum=as.numeric(x$SHT_Rhum))
    # FILE ERRORS:
    oklab_tmp$time <- as.POSIXlt(oklab_tmp$datetime,format="%Y-%m-%d %H:%M:%S",tz="UTC")
    oklab <- rbind(oklab,oklab_tmp)
}


fname <- paste0("oklab_",slist[1],"_oklab.Rdata")
#save(oklab,file=fname,compress=T,compression_level=9)
#load(fname)

message("reading done!")
asopc <- oklab


###################################################################
###################################################################
# CREATE DATA FRAME OF aggregation INTERVALL TIME STEPS
timesteps <- data.frame( time=seq(from=as.POSIXct(t1),to=as.POSIXct(t2),by=aggregation))
aggpm10 <- timesteps
aggtemp <- timesteps
aggrhum <- timesteps
aggpm25 <- timesteps

ids <- sort(unique(asopc$id))
for (i in seq_along(ids)) {
   #message(" ")
   message("id =",ids[i])
   idserpm10 <- asopc$pm10[asopc$id==ids[i]]
   idserpm25 <- asopc$pm25[asopc$id==ids[i]]
   idsertemp <- as.numeric(asopc$temp[asopc$id==ids[i]])
   idserrhum <- as.numeric(asopc$rhum[asopc$id==ids[i]])
   idtime <- asopc$time[asopc$id==ids[i]]

   tmppm10 <- data.frame(time=idtime,var=idserpm10)
   tmppm25 <- data.frame(time=idtime,var=idserpm25)
   tmptemp <- data.frame(time=idtime,var=as.numeric(idsertemp))
   tmprhum <- data.frame(time=idtime,var=as.numeric(idserrhum))

   # FILTER
   tmppm10$var[tmppm10$var > 60 ] <- NA
   tmptemp$var[tmptemp$var < -30 ] <- NA

   cuts <- cut( idtime,breaks=timesteps$time)
   if ( all(is.na(cuts)) ) {
     aggpm10$var <- NA
     names(aggpm10)[length(names(aggpm10))] <- paste0("id",ids[i])
     aggpm25$var <- NA
     names(aggpm25)[length(names(aggpm25))] <- paste0("id",ids[i])
     aggtemp$var <- NA
     names(aggtemp)[length(names(aggtemp))] <- paste0("id",ids[i])
     aggrhum$var <- NA
     names(aggrhum)[length(names(aggrhum))] <- paste0("id",ids[i])
     next
   }
   tmppm10$cuts <- cuts
   tmppm25$cuts <- cuts
   tmptemp$cuts <- cuts
   tmprhum$cuts <- cuts

   #message(paste("tmppm10: ",length(tmppm10)))
   #if( all(is.na(tmppm10$var)) ){
   if( sum( is.na(tmppm10$var) | is.na(tmppm10$cuts) ) == length(tmppm10$cuts) ){
   aggpm10$var <- NA   
   names(aggpm10)[length(names(aggpm10))] <- paste0("id",ids[i])
   } else {
   agg <- aggregate(. ~ cuts, tmppm10, mean)
   agg$time <- as.POSIXct(agg$cuts)
   reg <- merge(timesteps,agg,by="time",all.x=T)
   aggpm10$var <- reg$var
   names(aggpm10)[length(names(aggpm10))] <- paste0("id",ids[i])
   }

   #message(paste("tmppm25: ",length(tmppm25)))
   if( sum( is.na(tmppm25$var) | is.na(tmppm25$cuts) ) == length(tmppm25$cuts) ){
   #if( all(is.na(tmppm25$var)) ){
   aggpm25$var <- NA
   names(aggpm25)[length(names(aggpm25))] <- paste0("id",ids[i])
   } else {
   #message("pm25:",str(tmppm25))
   agg <- aggregate(. ~ cuts, tmppm25, mean)
   agg$time <- as.POSIXct(agg$cuts)
   reg <- merge(timesteps,agg,by="time",all.x=T)
   aggpm25$var <- reg$var
   names(aggpm25)[length(names(aggpm25))] <- paste0("id",ids[i])
   }

   #message(paste("tmptemp: ",str(tmptemp)))
   #if( all(is.na(tmptemp$var)) ){
   if( sum( is.na(tmptemp$var) | is.na(tmptemp$cuts) ) == length(tmptemp$cuts) ){
   aggtemp$var <- NA
   names(aggtemp)[length(names(aggtemp))] <- paste0("id",ids[i])
   } else {
   #message("temp:",str(tmptemp))
   agg <- aggregate(. ~ cuts, tmptemp, mean)
   agg$time <- as.POSIXct(agg$cuts)
   reg <- merge(timesteps,agg,by="time",all.x=T)
   aggtemp$var <- reg$var
   names(aggtemp)[length(names(aggtemp))] <- paste0("id",ids[i])
   }

   #message(paste("tmprhum: ",length(tmprhum)))
   #if( all(is.na(tmprhum$var)) ){
   if( sum( is.na(tmprhum$var) | is.na(tmprhum$cuts) ) == length(tmprhum$cuts) ){
   aggrhum$var <- NA
   names(aggrhum)[length(names(aggrgum))] <- paste0("id",ids[i])
   } else {
   #message("rhum:",str(tmprhum))
   agg <- aggregate(. ~ cuts, tmprhum, mean)
   agg$time <- as.POSIXct(agg$cuts)
   reg <- merge(timesteps,agg,by="time",all.x=T)
   aggrhum$var <- reg$var
   names(aggrhum)[length(names(aggrhum))] <- paste0("id",ids[i])
   }

}
message("aggregation done!")


fname <- paste0("oklab_",slist[1],"_oklab.Rdata")
save(oklab,aggpm10,aggpm25,aggtemp,aggrhum,file=fname,compress=T,compression_level=9)
#load(fname)

next

###################################################################
# FIT AND COMPARE SENSOR VALUES TO REFERENCE
message("correction and fitting to reference data ...")
#library(qmap)
library(randomForest)
#pdf(file="scatter.pdf",width=1900, height=600, pointsize=3)
pdfname <- paste0("oklab_",slist[1],"_scatter.pdf")
pdf(file=pdfname, pointsize=10)
for (i in seq_along(ids)) {
for (j in seq_along(slist)) {
if ( ids[i] == slist[j] ) {

  ############################################################
  # PM10 --------------------------------
  message()
  message("pm10 ...")

  pm10org <- aggpm10[,i+1]
  pm10ref <- aggpm10$'idOPC004'
  size <- c(0,max(pm10ref[!is.na(pm10ref)]))
  plot(pm10ref,pm10org,xlab="PM10 reference EDM164 OPC004",ylab=paste0("PM10 SDS011 ",slist[j]),xlim=size,ylim=size)
  abline(0,1)
  
#  # humidity correction 
#  rh <- aggrhum[,i+1] / 100.0
#  kap <- 0.5
#  c <- 1 + ( kap/1.65 ) / (-1 + 1/rh)
#  pm10cor <- aggpm10[,i+1] / c
#  points(pm10ref,pm10cor,col="blue")
  
#  # fitting by linear model
#  m <- lm(aggpm10$idOPC004 ~ pm10cor)
#  pm10reg <- m$coefficients[1] + pm10cor *  m$coefficients[2]
#  points(pm10ref,pm10reg,col="green")
  
#  # fitting quantile mapping
#  ref <- pm10ref[ !is.na(pm10ref) & !is.na(pm10cor)]
#  cor <- pm10cor[ !is.na(pm10ref) & !is.na(pm10cor)]
#  pm10qmap <- fitQmap( ref, cor, method="QUANT",wet.day=F)
#  pm10fit <- pm10cor
#  pm10fit[!is.na(pm10cor)] <- doQmap(pm10cor[!is.na(pm10cor)],pm10qmap)
#  points(pm10ref,pm10fit,col="red")
  
#  # fitting by random forests using humidity corrected sds011 values
#  set.seed(2)
#  ref <- pm10ref[ !is.na(pm10ref) & !is.na(pm10cor) ]
#  obs <- pm10cor[ !is.na(pm10ref) & !is.na(pm10cor) ]
#  # random subsets for 2/3 calibration and 1/3 validation
#  mask <- runif(length(ref))<0.66
#  ref_cal <- ref[mask]
#  ref_val <- ref[!mask]
#  obs_cal <- obs[mask]
#  obs_val <- obs[!mask]
#  #rfdata <- data.frame(ref=ref,obs=obs)
#  rfdata <- data.frame(ref=ref_cal,obs=obs_cal)
#  rf <-randomForest(ref~.,data=rfdata, ntree=1000) 
#  print(rf)
#  rfdata <- data.frame(ref=ref_val,obs=obs_val)
#  pm10rf <- predict(rf, rfdata)
#  points(rfdata$ref,pm10rf,col="darkgray")


  # fitting by random forests using sds011 and humidity as predictor
  set.seed(2)
  rhum <- aggrhum[,i+1]
  hum <- rhum   [ !is.na(pm10ref) & !is.na(pm10org) & !is.na(rhum)]
  ref <- pm10ref[ !is.na(pm10ref) & !is.na(pm10org) & !is.na(rhum)]
  obs <- pm10org[ !is.na(pm10ref) & !is.na(pm10org) & !is.na(rhum)]
  pm10rftime <- as.POSIXlt(aggpm10$time[ !is.na(pm10ref) & !is.na(pm10org) & !is.na(rhum)])
  
  # random subsets for 2/3 calibration and 1/3 validation
  #mask <- runif(length(ref))<0.66
  #mask <- seq(1:length(ref))<length(ref)*0.66
  mask <- seq(1:length(ref))<=length(ref) # all true
  tc1 <- as.POSIXlt("2022-07-12 15:45:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  tc2 <- as.POSIXlt("2022-07-13 18:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  mask[as.numeric(pm10rftime) >= as.numeric(tc1) & as.numeric(pm10rftime) <= as.numeric(tc2)] <- FALSE
  
  ref_cal <- ref[mask]
  ref_val <- ref[!mask]
  obs_cal <- obs[mask]
  obs_val <- obs[!mask]
  hum_cal <- hum[mask]
  hum_val <- hum[!mask]
  pm10time_cal <- pm10rftime[mask]
  pm10time_val <- pm10rftime[!mask]
  
  #rfhdata <- data.frame(ref=ref,obs=obs,hum=hum)
  rfdata <- data.frame(ref=ref_cal,obs=obs_cal,hum=hum_cal)
  rf <-randomForest(ref~.,data=rfdata, ntree=1000) 
  print(rf)
  rfdata <- data.frame(obs=obs_val,hum=hum_val)
  pm10rf <- predict(rf, rfdata)
  points(ref_val,pm10rf,col="red")
  
  # error estimation
  bias <- mean(pm10rf) - mean(ref_val)
  biaspercent <- ( bias / mean(ref_val) ) *100
  mae <- mean(abs(pm10rf - ref_val),na.rm=T)
  cor2 <- cor(pm10rf,ref_val,use="complete.obs")**2
  prec <- sqrt( sum((pm10rf - ref_val)**2)/length(ref_val) )
  precpercent <- ( prec / mean(ref_val) ) *100
  accuracy <- sqrt(prec**2 - bias**2) # Foken

  message("bias  = ",bias)
  message("bias% = ",biaspercent)
  message("mae   = ",mae)
  message("cor²  = ",cor2)
  message("prec  = ",prec)
  message("prec% = ",precpercent)
  
  bias <- as.character( format(round(bias, 2), nsmall = 2) )
  biaspercent <- as.character( format(round(biaspercent, 2), nsmall = 2) )
  mae <- as.character( format(round(mae, 2), nsmall = 2) )
  cor2 <- as.character( format(round(cor2, 2), nsmall = 2) )
  prec <- as.character( format(round(prec, 2), nsmall = 2) )
  precpercent <- as.character( format(round(precpercent, 2), nsmall = 2) )
  
#  mtext(paste("bias =",bias," MAE =",mae," r² =",cor2,side=3))
  mtext(paste0("Bias=",bias," (",biaspercent,"%)"," MAE=",mae," r²=",cor2," Precision=",prec," (",precpercent,"%)"),side=3)

  
  rfpm10 <- rf
  #rfpm10data <- rfdata
  rfpm10data <- data.frame(time=pm10time_val,ref=ref_val,obs=obs_val,cal=pm10rf)
  fname <- paste0("oklab_",slist[1],"_calib_pm10.Rdata")
  save(rfpm10,rfpm10data,file=fname,compress=T,compression_level=9)

  ############################################################
  # PM25 --------------------------------
  message()
  message("pm25 ...")

  pm25org <- aggpm25[,i+1]
  pm25ref <- aggpm25$'idOPC004'
  size <- c(0,max(pm25ref[!is.na(pm25ref)]))
  plot(pm25ref,pm25org,xlab="PM2.5 reference EDM164 OPC004",ylab=paste0("PM2.5 SDS011 ",slist[j]),xlim=size,ylim=size)
  abline(0,1)
  
  # humidity correction 
  rh <- aggrhum[,i+1] / 100.0
  kap <- 0.5
  c <- 1 + ( kap/1.65 ) / (-1 + 1/rh)
  pm25cor <- aggpm25[,i+1] / c
  points(pm25ref,pm25cor,col="blue")
  
#  # fitting by linear model
#  m <- lm(pm25ref ~ pm25cor)
#  pm25reg <- m$coefficients[1] + pm25cor *  m$coefficients[2]
#  points(pm25ref,pm25reg,col="green")
  
#  # quantile mapping
#  ref <- pm25ref[ !is.na(pm25ref) & !is.na(pm25cor)]
#  cor <- pm25cor[ !is.na(pm25ref) & !is.na(pm25cor)]
#  pm25qmap <- fitQmap( ref, cor, method="QUANT",wet.day=F)
#  pm25fit <- pm25cor
#  pm25fit[!is.na(pm25cor)] <- doQmap(pm25cor[!is.na(pm25cor)],pm25qmap)  
#  points(pm25ref,pm25fit,col="red")
#  mae <- as.character(mean(abs(pm25fit - pm25ref),na.rm=T))
#  cor2 <- as.character(cor(pm25fit,pm25ref,use="complete.obs")**2)
 
  # fitting by random forests using sds011 and humidity as predictor
  set.seed(2)
  rhum <- aggrhum[,i+1]
  hum <- rhum   [ !is.na(pm25ref) & !is.na(pm25org) & !is.na(rhum)]
  ref <- pm25ref[ !is.na(pm25ref) & !is.na(pm25org) & !is.na(rhum)]
  obs <- pm25org[ !is.na(pm25ref) & !is.na(pm25org) & !is.na(rhum)]
  pm25rftime <- aggpm25$time[ !is.na(pm25ref) & !is.na(pm25org) & !is.na(rhum)]

  # random subsets for 2/3 calibration and 1/3 validation
  #mask <- runif(length(ref))<0.66
  #mask <- seq(1:length(ref))<length(ref)*0.66
  mask <- seq(1:length(ref))<=length(ref) # all true
  tc1 <- as.POSIXlt("2022-07-12 15:45:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  tc2 <- as.POSIXlt("2022-07-13 18:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  mask[pm25rftime>=tc1 & pm25rftime<=tc2] <- FALSE



  ref_cal <- ref[mask]
  ref_val <- ref[!mask]
  obs_cal <- obs[mask]
  obs_val <- obs[!mask]
  hum_cal <- hum[mask]
  hum_val <- hum[!mask]
  pm25time_cal <- pm25rftime[mask]
  pm25time_val <- pm25rftime[!mask]

  rfdata <- data.frame(ref=ref_cal,obs=obs_cal,hum=hum_cal)
  rf <-randomForest(ref~.,data=rfdata, ntree=1000) 
  print(rf)
  rfdata <- data.frame(obs=obs_val,hum=hum_val)
  pm25rf <- predict(rf, rfdata)
  points(ref_val,pm25rf,col="red")
   
  # error estimation
  bias <- mean(pm25rf) - mean(ref_val)
  biaspercent <- ( bias / mean(ref_val) ) *100
  mae <- mean(abs(pm25rf - ref_val),na.rm=T)
  cor2 <- cor(pm25rf,ref_val,use="complete.obs")**2
  prec <- sqrt( sum((pm25rf - ref_val)**2)/length(ref_val) )
  precpercent <- ( prec / mean(ref_val) ) *100
  accuracy <- sqrt(prec**2 - bias**2) # Foken

  message("bias  = ",bias)
  message("bias% = ",biaspercent)
  message("mae   = ",mae)
  message("cor²  = ",cor2)
  message("prec  = ",prec)
  message("prec% = ",precpercent)
  
  bias <- as.character( format(round(bias, 2), nsmall = 2) )
  biaspercent <- as.character( format(round(biaspercent, 2), nsmall = 2) )
  mae <- as.character( format(round(mae, 2), nsmall = 2) )
  cor2 <- as.character( format(round(cor2, 2), nsmall = 2) )
  prec <- as.character( format(round(prec, 2), nsmall = 2) )
  precpercent <- as.character( format(round(precpercent, 2), nsmall = 2) )

  
#  mtext(paste("bias =",bias," MAE =",mae," r² =",cor2,side=3))
  mtext(paste0("Bias=",bias," (",biaspercent,"%)"," MAE=",mae," r²=",cor2," Precision=",prec," (",precpercent,"%)"),side=3)

  
  rfpm25 <- rf
  #rfpm25data <- rfdata
  rfpm25data <- data.frame(time=pm25time_val,ref=ref_val,obs=obs_val,cal=pm25rf)
  fname <- paste0("oklab_",slist[1],"_calib_pm25.Rdata")
  save(rfpm25,rfpm25data,file=fname,compress=T,compression_level=9)


  ############################################################
  # TEMP --------------------------------
  message()
  message("temp ...")
  temporg <- aggtemp[,i+1]
  tempref <- aggtemp$id177410813

  size <- c(min(c(tempref[!is.na(tempref)],temporg[!is.na(temporg)]))   ,  max(c(tempref[!is.na(tempref)],temporg[!is.na(temporg)])))
  plot(aggtemp$id177410813,aggtemp[,i+1],xlab="Temp reference SHT 177410813",ylab=paste0("Temp BME280 ",slist[j]),xlim=size,ylim=size)
  abline(0,1)
  
  m <- lm(aggtemp$id177410813 ~ aggtemp[,i+1])
  mcor <- m$coefficients[1] + aggtemp[,i+1] *  m$coefficients[2]
  points(aggtemp$id177410813,mcor,col="red")
  mae <- as.character(mean(abs(mcor - aggtemp$id177410813),na.rm=T))
  cor2 <- as.character(cor(mcor,aggtemp$id177410813,use="complete.obs")**2)
  
  # fitting by random forests using sds011 and humidity as predictor
  set.seed(2)
  obs <- temporg[ !is.na(temporg) & !is.na(tempref)]
  ref <- tempref[ !is.na(temporg) & !is.na(tempref)]
  temprftime <- aggtemp$time[ !is.na(temporg) & !is.na(tempref)]
  
  # random subsets for 2/3 calibration and 1/3 validation
  #mask <- runif(length(ref))<0.66
  mask <- seq(1:length(ref))<length(ref)*0.66
  mask <- seq(1:length(ref))<=length(ref) # all true
  tc1 <- as.POSIXlt("2022-07-12 15:45:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  tc2 <- as.POSIXlt("2022-07-13 18:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  mask[temprftime>=tc1 & temprftime<=tc2] <- FALSE
  
  
  ref_cal <- ref[mask]
  ref_val <- ref[!mask]
  obs_cal <- obs[mask]
  obs_val <- obs[!mask]
  temptime_cal <- temprftime[mask]
  temptime_val <- temprftime[!mask]  
  
  
  rfdata <- data.frame(ref=ref_cal,obs=obs_cal)
  rf <-randomForest(ref~.,data=rfdata, ntree=1000) 
  print(rf)
  rfdata <- data.frame(ref=ref_val,obs=obs_val)
  temprf <- predict(rf, rfdata)
  points(ref_val,temprf,col="red")    
  
  bias <- mean(temprf) - mean(ref_val)
  biaspercent <- ( bias / mean(ref_val) ) *100
  mae <- mean(abs(temprf - ref_val),na.rm=T)
  cor2 <- cor(temprf,ref_val,use="complete.obs")**2
  prec <- sqrt( sum((temprf - ref_val)**2)/length(ref_val) )
  precpercent <- ( prec / mean(ref_val) ) *100
  accuracy <- sqrt(prec**2 - bias**2) # Foken

  message("bias  = ",bias)
  message("bias% = ",biaspercent)
  message("mae   = ",mae)
  message("cor²  = ",cor2)
  message("prec  = ",prec)
  message("prec% = ",precpercent)


  bias <- as.character( format(round(bias, 2), nsmall = 2) )
  biaspercent <- as.character( format(round(biaspercent, 2), nsmall = 2) )
  mae <- as.character( format(round(mae, 2), nsmall = 2) )
  cor2 <- as.character( format(round(cor2, 2), nsmall = 2) )
  prec <- as.character( format(round(prec, 2), nsmall = 2) )
  precpercent <- as.character( format(round(precpercent, 2), nsmall = 2) )


  #mtext(paste("bias =",bias,"MAE =",mae," r² =",cor2,side=3))
  mtext(paste0("Bias=",bias," (",biaspercent,"%)"," MAE=",mae," r²=",cor2," Precision=",prec," (",precpercent,"%)"),side=3)

  
  rftemp <- rf
  #rftempdata <- rfdata
  rftempdata <- data.frame(time=temptime_val,ref=ref_val,obs=obs_val,cal=temprf)
  fname <- paste0("oklab_",slist[1],"_calib_temp.Rdata")
  save(rftemp,rftempdata,file=fname,compress=T,compression_level=9)

  ############################################################
  # RHUM --------------------------------
  message()
  message("rhum ...")
  
  rhumorg <- aggrhum[,i+1]
  rhumref <- aggrhum$id177410813

  size <- c(min(c(rhumorg[!is.na(rhumorg)],rhumref[!is.na(rhumref)]))   ,max(c(rhumorg[!is.na(rhumorg)],rhumref[!is.na(rhumref)])))
  plot(aggrhum$id177410813,aggrhum[,i+1],xlab="Rhum reference SHT 177410813",ylab=paste0("Rhum BME280 ",slist[j]),xlim=size,ylim=size)
  abline(0,1)
 
  #m <- lm(aggrhum$id177410813 ~ aggrhum[,i+1])
  #mcor <- m$coefficients[1] + aggrhum[,i+1] *  m$coefficients[2]
  #points(aggrhum$id177410813,mcor,col="red")
  #mae <- as.character(mean(abs(mcor - aggrhum$id177410813),na.rm=T))
  #cor2 <- as.character(cor(mcor,aggrhum$id177410813,use="complete.obs")**2)
  #mtext(paste("MAE =",mae," r² =",cor2,side=3))
  
  # fitting by random forests using sds011 and humidity as predictor
  set.seed(2)
  obs <- rhumorg[ !is.na(rhumorg) & !is.na(rhumref)]
  ref <- rhumref[ !is.na(rhumorg) & !is.na(rhumref)]
  rhumrftime <- aggrhum$time[ !is.na(rhumorg) & !is.na(rhumref)]
  
  # random subsets for 2/3 calibration and 1/3 validation
  #mask <- runif(length(ref))<0.66
  #mask <- seq(1:length(ref))<length(ref)*0.66
  mask <- seq(1:length(ref))<=length(ref) # all true
  tc1 <- as.POSIXlt("2022-07-12 15:45:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  tc2 <- as.POSIXlt("2022-07-13 18:00:00",format="%Y-%m-%d %H:%M:%S",tz="UTC")
  mask[rhumrftime>=tc1 & rhumrftime<=tc2] <- FALSE


  ref_cal <- ref[mask]
  ref_val <- ref[!mask]
  obs_cal <- obs[mask]
  obs_val <- obs[!mask]
  rhumtime_cal <- rhumrftime[mask]
  rhumtime_val <- rhumrftime[!mask]  
  
  rfdata <- data.frame(ref=ref_cal,obs=obs_cal)
  rf <-randomForest(ref~.,data=rfdata, ntree=1000) 
  print(rf)
  rfdata <- data.frame(ref=ref_val,obs=obs_val)
  rhumrf <- predict(rf, rfdata)
  
  points(ref_val,rhumrf,col="red")    
  bias <- mean(rhumrf) - mean(ref_val)
  biaspercent <- ( bias / mean(ref_val) ) *100
  mae <- mean(abs(rhumrf - ref_val),na.rm=T)
  cor2 <- cor(rhumrf,ref_val,use="complete.obs")**2
  prec <- sqrt( sum((rhumrf - ref_val)**2)/length(ref_val) )
  precpercent <- ( prec / mean(ref_val) ) *100
  accuracy <- sqrt(prec**2 - bias**2) # Foken
  
  message("bias  = ",bias)
  message("bias% = ",biaspercent)
  message("mae   = ",mae)
  message("cor²  = ",cor2)
  message("prec  = ",prec)
  message("prec% = ",precpercent)
  message("acc   = ",accuracy)


  bias <- as.character( format(round(bias, 2), nsmall = 2) )
  biaspercent <- as.character( format(round(biaspercent, 2), nsmall = 2) )
  mae <- as.character( format(round(mae, 2), nsmall = 2) )
  cor2 <- as.character( format(round(cor2, 2), nsmall = 2) )
  prec <- as.character( format(round(prec, 2), nsmall = 2) )
  precpercent <- as.character( format(round(precpercent, 2), nsmall = 2) )


  #mtext(paste("bias =",bias,"MAE =",mae," r² =",cor2,side=3))
  mtext(paste0("Bias=",bias," (",biaspercent,"%)"," MAE=",mae," r²=",cor2," Precision=",prec," (",precpercent,"%)"),side=3)
  
  rfrhum <- rf
  #rfrhumdata <- rfdata
  rfrhumdata <- data.frame(time=rhumtime_val,ref=ref_val,obs=obs_val,cal=rhumrf)
  fname <- paste0("oklab_",slist[1],"_calib_rhum.Rdata")
  save(rfrhum,rfrhumdata,file=fname,compress=T,compression_level=9)
  
}}}
dev.off()


###################################################################
# PLOT TIME SERIES PM10
message("plotting time series ...")

pngname <- paste0("oklab_",slist[1],"_pm10.png")
png(pngname,width=1900, height=600, units="px", res=300, pointsize=3)
tmin <- min(aggpm10$time)
tmax <- max(aggpm10$time)
message("tmin =",tmin)
message("tmax =",tmax)

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

#message("tstart =",tstart)
#message("tstop  =",tstop)
aggrmax <- max(aggpm10[,3:ncol(aggpm10)], na.rm=TRUE)
allmax <- min(max(aggrmax),60)
#message("ylim =",0," ",allmax)

ts <- aggpm10[,3]
plot(aggpm10$time[aggpm10$time>=tc1 & aggpm10$time<=tc2],ts[aggpm10$time>=tc1 & aggpm10$time<=tc2],xlab=paste("time UTC"),ylab="pm10 (µg/m³)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(0,allmax),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("PM10 oklab",aggregation,"mean (",tmin," - ",tmax," UTC)" ))

ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , units="hours")
axis.POSIXct(1, x=aggpm10$time, at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(0,allmax,by=10),lwd=0.2,col="gray")

#pallette <- rainbow(length(ids))
#pallette <- distinctColorPalette(length(ids))
pallette <- c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', '#ffffff', '#000000')

leglist <- list()
lines(aggpm10$time,aggpm10$idOPC004,col="grey30",t="l",lwd=1.5)
leglist <- c(leglist,"EDM164 OPC004")
for (i in seq_along(ids)) {
for (j in seq_along(slist)) {
  if( ids[i] == slist[j] ){
  leglist <- c(leglist,paste("SDS011",ids[i]))
  lines(aggpm10$time,aggpm10[,i+1],col="blue",t="l",lwd=0.3)
  }
}}
lines(pm10time_val,pm10rf,col="red",lwd=0.6)
leglist <- c(leglist,"Recalibration")
legend("topleft",legend=c(leglist), col=c("grey30","blue","red"), lty=1,lwd=0.5)

dev.off()
message("pm10 done")


###################################################################
# PLOT TIME SERIES PM25
pngname <- paste0("oklab_",slist[1],"_pm25.png")
png(pngname,width=1900, height=600, units="px", res=300, pointsize=3)
tmin <- min(aggpm25$time)
tmax <- max(aggpm25$time)
#message("tmin =",tmin)
#message("tmax =",tmax)
#tstart <- min(seq(tmin, tmax, by = "hour"))
#tstop <- max(seq(tmin, tmax, by = "hour"))
tstart <- min(seq(t1, t2, by = "hour"))
tstop <- max(seq(t1, t2, by = "hour"))
message("tstart = ",tstart)
message("tstop = ",tstop)


#message("tstart =",tstart)
#message("tstop  =",tstop)
aggrmax <- max(aggpm25[,2:ncol(aggpm25)], na.rm=TRUE)
allmax <- min(max(aggrmax),50)
#allmax <- aggrmax
ts <- aggpm25[,2]
plot(aggpm25$time[aggpm25$time>=tc1 & aggpm25$time<=tc2],ts[aggpm25$time>=tc1 & aggpm25$time<=tc2],xlab=paste("time UTC"),ylab="pm2.5 (µg/m³)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(0,allmax),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("PM2.5 oklab",aggregation,"mean (",tmin," - ",tmax," UTC)" ))
ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , units="hours")
axis.POSIXct(1, x=aggpm25$time, at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(0,allmax,by=10),lwd=0.2,col="gray")

pallette <- c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', '#ffffff', '#000000')

leglist <- list()
lines(aggpm25$time,aggpm25$idOPC004,col="grey20",t="l",lwd=1.5)
leglist <- c(leglist,"EDM164 OPC004")
for (i in seq_along(ids)) {
for (j in seq_along(slist)) {
  if( ids[i] == slist[j] ){
  leglist <- c(leglist,paste("SDS011",ids[i]))
  lines(aggpm25$time,aggpm10[,i+1],col="blue",t="l",lwd=0.3)
  }
}}
lines(pm25time_val,pm25rf,col="red",lwd=0.6)
leglist <- c(leglist,"SDS011 Recalibration")
legend("topleft",legend=c(leglist), col=c("grey20","blue","red"), lty=1,lwd=0.5)

dev.off()


###################################################################
# PLOT TIME SERIES RHUM
pngname <- paste0("oklab_",slist[1],"_rhum.png")
png(pngname,width=1900, height=600, units="px", res=300, pointsize=3)
#tmin <- min(aggrhum$time)
#tmax <- max(aggrhum$time)
tstart <- min(seq(t1, t2, by = "hour"))
tstop <- max(seq(t1, t2, by = "hour"))
message("tstart = ",tstart)
message("tstop = ",tstop)

#message("tmin =",tmin)
#message("tmax =",tmax)
tstart <- min(seq(tmin, tmax, by = "hour"))
tstop <- max(seq(tmin, tmax, by = "hour"))
#message("tstart =",tstart)
#message("tstop  =",tstop)
allmax <- 100
ts <- aggrhum[,2]
plot(aggrhum$time[aggrhum$time>=tc1 & aggrhum$time<=tc2],ts[aggrhum$time>=tc1 & aggrhum$time<=tc2],xlab=paste("time UTC"),ylab="rhum (%)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(0,allmax),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("RHUM oklab",aggregation,"mean (",tmin," - ",tmax," UTC)" ))
ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , units="hours")
axis.POSIXct(1, x=aggrhum$time, at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(0,allmax,by=10),lwd=0.2,col="gray")

#LWD <- 1.0
#lines(aggrhum$time,aggrhum$id177410813,col="black",lwd=LWD)
#lines(aggrhum$time,aggrhum$id176430103,col="darkgrey",lwd=LWD)

pallette <- c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', '#ffffff', '#000000')

leglist <- list()
lines(aggrhum$time,aggrhum$id177410813,col="grey20",t="l",lwd=1.5)
leglist <- c(leglist,"SHT85 177410813")

for (i in seq_along(ids)) {
for (j in seq_along(slist)) {
  if( ids[i] == slist[j] ){
  leglist <- c(leglist,paste("BME280",ids[i]))
  lines(aggrhum$time,aggrhum[,i+1],col="blue",t="l",lwd=0.3)
  }
}}
lines(rhumtime_val,rhumrf,col="red",lwd=0.6)
leglist <- c(leglist,"BME280 Recalibration")
legend("topleft",legend=c(leglist), col=c("grey20","blue","red"), lty=1,lwd=0.5)
dev.off()

###################################################################
# PLOT TIME SERIES TEMP
pngname <- paste0("oklab_",slist[1],"_temp.png")
png(pngname,width=1900, height=600, units="px", res=300, pointsize=3)
tmin <- min(aggtemp$time)
tmax <- max(aggtemp$time)
#message("tmin =",tmin)
#message("tmax =",tmax)

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

#message("tstart =",tstart)
#message("tstop  =",tstop)
#allmax <- 100

aggrmax <- max(aggtemp[,2:ncol(aggtemp)], na.rm=TRUE)
aggrmin <- min(aggtemp[,2:ncol(aggtemp)], na.rm=TRUE)
allmax <- ceiling(aggrmax) #min(max(aggrmax,luebmax),70)
allmin <- floor(aggrmin) #min(max(aggrmax,luebmax),70)

ts <- aggtemp[,2]
plot(aggtemp$time[aggtemp$time>=tc1 & aggtemp$time<=tc2],ts[aggtemp$time>=tc1 & aggtemp$time<=tc2],xlab=paste("time UTC"),ylab="temp (°C)",
xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),ylim=c(allmin,allmax),type="n",xaxs="i",yaxs="i",xaxt="n",
main=paste("TEMP oklab",aggregation,"mean (",tmin," - ",tmax," UTC)" ))
ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , units="hours")
#axis.POSIXct(1, x=aggtemp$time, at = ihours, format = "%a %m-%d %H:%M")
axis.POSIXct(1, x=seq(tstart,tstop,by=3600), at = ihours, format = "%a %m-%d %H:%M")
abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray")
abline(h=seq(allmin,allmax,by=1),lwd=0.2,col="gray")

leglist <- list()
lines(aggtemp$time,aggtemp$id177410813,col="grey20",lwd=1.5)
leglist <- c(leglist,"SHT85 177410813")
for (i in seq_along(ids)) {
for (j in seq_along(slist)) {
  if( ids[i] == slist[j] ){
  leglist <- c(leglist,paste("BME280",ids[i]))
  lines(aggtemp$time,aggtemp[,i+1],col="blue",t="l",lwd=0.3)
  }
}}
lines(temptime_val,temprf,col="red",lwd=0.6)
leglist <- c(leglist,"BME280 Recalibration")
legend("topleft",legend=c(leglist), col=c("grey20","blue","red"), lty=1,lwd=0.5)
dev.off()

oklabaggpm10 <- aggpm10
oklabaggpm25 <- aggpm25
oklabaggtemp <- aggtemp
oklabaggrhum <- aggrhum

save(oklab,oklabaggpm10,oklabaggpm25,oklabaggrhum,oklabaggtemp,file=paste0("oklab_",slist[1],".Rdata"))
message("writing oklab.Rdata done!")