# SELECT SENSOR-ID
id <- "id12776226"

# 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)))
temp[temp < -999] <- NA
rhum <- eval(parse(text=paste0("oklabaggrhum$",id)))
rhum[rhum < -999] <- NA
time <- as.POSIXlt(oklabaggpm10$time,tz="UTC")


# TIME PERIOD
t1 <- min(time)
t2 <- max(time)

# 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")
pm10ref <- edmaggpm10$idOPC004
pm10reftime <- as.POSIXlt(edmaggpm10$time,tz="UTC")
pm25ref <- edmaggpm25$idOPC004
pm25reftime <- as.POSIXlt(edmaggpm25$time,tz="UTC")

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

# PLOT PM10
message("plotting time series for pm10 ...")
png(paste0(id,"_pm10_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
pm10max <- min( max(pm10,na.rm=TRUE) , 60 )
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("PM10 oklab 5min-mean (",tmin," - ",tmax," UTC)" ))
ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , 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 time series for pm25 ...")
png(paste0(id,"_pm25_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
pm25max <- min( max(pm25,na.rm=TRUE) , 50 )
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("PM25 oklab 5min-mean (",tmin," - ",tmax," UTC)" ))
ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , 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()

# PLOT TEMP
message("plotting time series for temp ...")
png(paste0(id,"_temp_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
tempmin <- min(temp,na.rm=TRUE)
tempmax <- max(temp,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("Temp oklab 5min-mean (",tmin," - ",tmax," UTC)" ))
ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , 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 time series for rhum ...")
png(paste0(id,"_rhum_timeser.png"),width=1900, height=600, units="px", res=300, pointsize=3)
rhummin <- min(rhum,na.rm=TRUE)
rhummax <- max(rhum,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("Rhum oklab 5min-mean (",tmin," - ",tmax," UTC)" ))
ihours <- round( seq( as.POSIXct(tstart), as.POSIXct(tstop),by=3600) , 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
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 SCATTER PLOTS
message("plotting scatter plot for temp ...")
png(paste0(id,"_pm10_scatter.png"),width=800, height=800, units="px", res=300, pointsize=3)
plot(okm$pm10,okm$pm10ref,xlim=c(0,pm10max),ylim=c(0,pm10max),xlab="PM10 SDS011 (µg/m³)",ylab="PM10 EDM164 (µg/m³)",main=paste("OKlab sensor ",id," versus reference"))
dev.off()