library(randomcoloR) library(data.table) library(lubridate) Sys.setenv(TZ = "UTC") ################################################################### # SELECT TIME PERIOD # seconds for 2 days daystep <- 60*60*24 step <- daystep*7 args <- commandArgs(trailingOnly = TRUE) message("ARGUMENTS = ",length(args)) if (length(args)==0) { t2 <- as.POSIXlt(Sys.time(),tz="UTC") t2 <- t2 + 60*60 t2 <- round_date(as.POSIXct(t2),"5 minutes") #t2 <- Sys.time() t1 <- t2 - step #"args:" length(args) } else { t1 <- as.POSIXlt(args[1]) t2 <- as.POSIXlt(args[2]) } # SET TIME PERIOD HERE: #t1 <- as.POSIXct("2022-11-22 00:00:00 CET") t1 <- as.POSIXct("2023-05-23 14:00:00 UTC") t2 <- as.POSIXlt(Sys.time(),tz="UTC") t2 <- round_date(as.POSIXct(t2),"5 minutes") 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") aggregation <- "60 sec" aggregation <- "1 min" aggregation <- "5 min" ################################################################### # GET MATCHING FILE NAMES FOR SELECTED DAYS #slist <- c("14155848","547446","2398386","14722068","2519718","5682602","14573303", #"12776226","547238","550119","14595225","5682100","5630365","660436","14462360","14595225") slist <- c("5628651","662424","13145084","14620483","14594730","550119","14595225","12776226","2398386","14722068","14573303","547238","12777335","14155848") #slist <- c("13145084","550119","14595225","14722068","14725731","12777335","14155848","2398386") #slist <- c("14462360","547446","5682100", # "660436","2519718","5630365", # "14722068","14155848","12776226") #"14155848","547446","2398386","14722068","2519718", # "5682602","12776226","547238","550119","14595225", # "5682100","660436","14462360", # "14573303","5630365") slist <- sort(slist) #listo <- as.numeric(slist) #slist <- slist[order(listo)] 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 asopc 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$pm10[x$pm10 < -99 ] <- NA x$pm25[x$pm25 < -99 ] <- NA #message(str(asopc)) #message(typeof(asopc)) if( typeof(x) == "character"){ next } if( is.data.frame(x) == F ){ next } #message("nrow = ",nrow(asopc)) if (nrow(x)<1) { next } #asopc$opctype <- opctype x$id <- id #asopc$shttype <- shttype #asopc$shtid <- shtid #asopc[nrow(asopc)+1,] <- NA #asopc <- merge(asopc,asopc_tmp,all.x =TRUE) 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) } #asopc$time <- as.POSIXct(asopc$UTC,format="%Y-%m-%d %H:%M:%S") #asopc$time <- as.POSIXct(asopc$datetime,format="%Y-%m-%d %H:%M:%S") - 60*60 oklab$time <- as.POSIXlt(oklab$datetime,format="%Y-%m-%d %H:%M:%S",tz="UTC") oklab[as.character(oklab)=="nan"] <- NA #asopc$PM10 <- as.numeric(asopc$PM10) #asopc$PM25 <- as.numeric(asopc$PM25) #asopc$PM1 <- as.numeric(asopc$PM1) #asopc$Tout <- as.numeric(asopc$Tout) #asopc$rH <- as.numeric(asopc$rH) #asopc$longitude <- as.numeric(asopc$longitude) #asopc$latitude <- as.numeric(asopc$latitude) #asopc$height <- as.numeric(asopc$height) message("reading done!") asopc <- oklab ################################################################### # CREATE DATA FRAME OF aggregation INTERVALL TIME STEPS #secs <- seq.POSIXt(from=as.POSIXct(t1),to=as.POSIXct(t2),by="sec") #secs <- data.frame( time=seq.POSIXt(from=as.POSIXct(t1),to=as.POSIXct(t2),by=aggregation)) #secs <- data.frame( time=seq(from=as.POSIXct(t1),to=as.POSIXct(t2),by=aggregation)) 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 <- c(0.0,asopc$pm10[asopc$id==ids[i]]) # idserpm25 <- c(0.0,asopc$pm25[asopc$id==ids[i]]) # idsertemp <- c(0.0,asopc$temp[asopc$id==ids[i]]) # idserrhum <- c(0.0,asopc$rhum[asopc$id==ids[i]]) idserpm10 <- asopc$pm10[asopc$id==ids[i]] idserpm25 <- asopc$pm25[asopc$id==ids[i]] idsertemp <- asopc$temp[asopc$id==ids[i]] idserrhum <- asopc$rhum[asopc$id==ids[i]] #idser[1] <- NA #idtime <- c(as.POSIXct(t1),asopc$time[asopc$opcid==ids[i]]) # idtime <- c(as.POSIXct(runningdate[1]),asopc$time[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=idsertemp) tmprhum <- data.frame(time=idtime,var=idserrhum) # FILTER tmppm10$var[tmppm10$var > 60 ] <- NA #tmp$cuts <- cut( idtime,breaks=aggregation) #cuts <- cut( idtime,breaks=timesteps$time) secstime1 <- as.numeric(timesteps$time) idtime1 <- as.numeric(idtime) cuts <- cut( idtime1,breaks=secstime1,labels=format(timesteps$time[1:length(secstime1)-1],"%Y-%m-%d %H:%M:%S UTC")) 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( sum( is.na(tmppm10$var) | is.na(tmppm10$cuts) ) == length(tmppm10$cuts) ){ #if( all(is.na(tmppm10$var)) ){ 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[1:length(reg$var)-1] 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[1:length(reg$var)-1] aggpm25$var <- reg$var names(aggpm25)[length(names(aggpm25))] <- paste0("id",ids[i]) } message(paste("tmptemp: ",str(tmptemp))) if( sum( is.na(tmptemp$var) | is.na(tmptemp$cuts) ) == length(tmptemp$cuts) ){ #if( all(is.na(tmptemp$var)) ){ 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[1:length(reg$var)-1] aggtemp$var <- reg$var names(aggtemp)[length(names(aggtemp))] <- paste0("id",ids[i]) } message(paste("tmprhum: ",length(tmprhum))) if( sum( is.na(tmprhum$var) | is.na(tmprhum$cuts) ) == length(tmprhum$cuts) ){ #if( all(is.na(tmprhum$var)) ){ aggrhum$var <- NA names(aggrhum)[length(names(aggrhum))] <- 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[1:length(reg$var)-1] aggrhum$var <- reg$var names(aggrhum)[length(names(aggrhum))] <- paste0("id",ids[i]) } } message("aggregation done!") ################################################################### ## GET LUEB REFERENCE DATA (in CET = Winterzeit = UTC + 1) ## aktuelle Daten offensichtlich in Lokalzeit! #tshift <- - 60*60*2 - 30 ## 1/2 hour since LUEB hourly means end at the given time stamp but represent the hour before # #path <- format(t2,"https://saqn.geo.uni-augsburg.de/lueb/%Y_%m/%Y_%m_%d_Bourges-Platz_PM10.csv") #bou <- read.csv(path,sep=";",header=F,stringsAsFactors=F,na.strings="---") #colnames(bou) <- c("tdate","ttime","pm10") #bou$time <- as.POSIXct(paste(bou$tdate," ",bou$ttime,sep=""),format="%d.%m.%Y %H:%M") #bou$time <- bou$time + tshift #bou$pm10[bou$pm10==":"] <- NA #bou$pm10 <- as.numeric(bou$pm10) # #path <- format(t2,"https://saqn.geo.uni-augsburg.de/lueb/%Y_%m/%Y_%m_%d_Karlstrasse_PM10.csv") #kar <- read.csv(path,sep=";",header=F,stringsAsFactors=F,na.strings="---") #colnames(kar) <- c("tdate","ttime","pm10") #kar$time <- as.POSIXct(paste(kar$tdate," ",kar$ttime,sep=""),format="%d.%m.%Y %H:%M") #kar$time <- kar$time + tshift #kar$pm10[kar$pm10==":"] <- NA #kar$pm10 <- as.numeric(kar$pm10) # #path <- format(t2,"https://saqn.geo.uni-augsburg.de/lueb/%Y_%m/%Y_%m_%d_Koenigsplatz_PM10.csv") #koe <- read.csv(path,sep=";",header=F,stringsAsFactors=F,na.strings="---") #colnames(koe) <- c("tdate","ttime","pm10") #koe$time <- as.POSIXct(paste(koe$tdate," ",koe$ttime,sep=""),format="%d.%m.%Y %H:%M") #koe$time <- koe$time + tshift #koe$pm10[koe$pm10==":"] <- NA #koe$pm10 <- as.numeric(koe$pm10) # #path <- format(t2,"https://saqn.geo.uni-augsburg.de/lueb/%Y_%m/%Y_%m_%d_LfU_PM10.csv") #lfu <- read.csv(path,sep=";",header=F,stringsAsFactors=F,na.strings="---") #colnames(lfu) <- c("tdate","ttime","pm10") #lfu$time <- as.POSIXct(paste(lfu$tdate," ",lfu$ttime,sep=""),format="%d.%m.%Y %H:%M") #lfu$time <- lfu$time + tshift #lfu$pm10[lfu$pm10==":"] <- NA #lfu$pm10 <- as.numeric(lfu$pm10) ################################################################### # LOAD LUEB, GRIMM AND ASOPC DATA load("/zugspitze/lappwald/saqn/recent/lueb.Rdata") luebmax <- max(c(lfu$pm10,koe$pm10,kar$pm10,bou$pm10),na.rm=T) luebmax <- ceiling(1+luebmax/10.0)*10 message("luebmax =",luebmax) load("/zugspitze/lappwald/saqn/recent/edm.Rdata") load("/zugspitze/lappwald/saqn/recent/asopc.Rdata") # delete measurement error t1error <- as.POSIXct("2022-11-27T03:29:44 CET") t2error <- as.POSIXct("2022-11-30T23:59:59 CET") asopcaggrhum$id177410813[asopcaggrhum$time >= t1error & asopcaggrhum$time <= t2error] <- NA asopcaggtemp$id177410813[asopcaggtemp$time >= t1error & asopcaggtemp$time <= t2error] <- NA asopcaggtemp$time <- as.POSIXlt(asopcaggtemp$time + 2*60*60,tz="UTC") asopcaggrhum$time <- as.POSIXlt(asopcaggrhum$time + 2*60*60,tz="UTC") #asopcaggtemp$time <- with_tz(asopcaggtemp$time,"UTC") #asopcaggrhum$time <- with_tz(asopcaggrhum$time,"UTC") ################################################################### # PLOT TIME SERIES PM10 png("oklab_pm10_exclim.png",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")) #message("tstart =",tstart) #message("tstop =",tstop) aggrmax <- max(aggpm10[,3:ncol(aggpm10)], na.rm=TRUE) allmax <- min(max(aggrmax,luebmax),60) message("ylim =",0," ",allmax) plot(aggpm10$time,aggpm10[,3],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") #lines(bou$time,bou$pm10,col="magenta",t="l",lwd=1.5) #lines(kar$time,kar$pm10,col="red",t="l",lwd=1.5) #lines(koe$time,koe$pm10,col="blue",t="l",lwd=1.5) #lines(lfu$time,lfu$pm10,col="green",t="l",lwd=1.5) #lines(edmaggpm10$time,edmaggpm10$idOPC002,col="orange",t="l",lwd=1.5) lines(edmaggpm10$time,edmaggpm10$idOPC004,col="black",t="l",lwd=1.5) #lines(edmaggpm10$time,edmaggpm10$idOPC006,col="yellow",t="l",lwd=1.5) #lines(edmaggpm10$time,edmaggpm10$idSN19010,col="grey",t="l",lwd=1.5) #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') for (i in seq_along(ids)) { LWD <- 0.3 if (ids[i] == "SN19014" || ids[i] == "SN19015" || ids[i] == "SN19017" || ids[i] == "SN19019" || ids[i] == "OPC004" || ids[i] == "OPC006" || ids[i] == "SN19010") { LWD <- 1.0 } if (ids[i] == "OPC004") pallette[i] <- "#000000" # black if (ids[i] == "SN19010") pallette[i] <- "#808080" # grey LTY <- 1 if (ids[i] == "SN19014" ) { LTY <- 3 LWD <- 0.7 } lines(aggpm10$time,aggpm10[,i+1],xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col=pallette[i],t="l",lwd=LWD,lty=LTY) } #legend("topright",c("LUEB-bou","LUEB-kar","LUEB-koe","LUEB-lfu","OPC002","OPC004","OPC006","SN19010"), # col=c("magenta","red","blue","green","orange","black","yellow","grey"),lty=1,lwd=0.5) legend("topright",c("OPC004"),col=c("black"),lty=1,lwd=0.5) legend("topleft",ids,col=pallette,lty=1,lwd=0.5) dev.off() message("pm10 done") ################################################################### # PLOT TIME SERIES PM25 png("oklab_pm25_exclim.png",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")) #message("tstart =",tstart) #message("tstop =",tstop) aggrmax <- max(aggpm25[,2:ncol(aggpm25)], na.rm=TRUE) allmax <- min(max(aggrmax,luebmax),50) plot(aggpm10$time,aggpm25[,2],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") #lines(bou$time,bou$pm25,col="magenta",t="l",lwd=0.6) #lines(kar$time,kar$pm25,col=2,t="l",lwd=0.6) #lines(koe$time,koe$pm25,col=3,t="l",lwd=0.6) #lines(lfu$time,lfu$pm25,col=4,t="l",lwd=0.6) for (i in seq_along(ids)) { LWD <- 0.3 if (ids[i] == "SN19014" || ids[i] == "SN19015" || ids[i] == "SN19017" || ids[i] == "SN19019" || ids[i] == "OPC-004") { LWD <- 1.0 } lines(aggpm25$time,aggpm25[,i+1],xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col=rainbow(length(ids))[i],t="l",lwd=LWD) } #lines(pm25lueb$time,pm25lueb$"Augsburg/Bourges-Platz",col="magenta",t="l",lwd=1.5) #lines(pm25lueb$time,pm25lueb$"Augsburg/LfU",col="green",t="l",lwd=1.5) #lines(edmaggpm25$time,edmaggpm25$idOPC002,col="orange",t="l",lwd=1.5) lines(edmaggpm25$time,edmaggpm25$idOPC004,col="black",t="l",lwd=1.5) #lines(edmaggpm25$time,edmaggpm25$idOPC006,col="yellow",t="l",lwd=1.5) #lines(edmaggpm25$time,edmaggpm25$idSN19010,col="grey",t="l",lwd=1.5) #legend("topright",c("LUEB-bou","LUEB-lfu","OPC002","OPC004","OPC006","SN19010"),col=c("magenta","green","orange","black","yellow","grey"),lty=1,lwd=0.5) legend("topright",c("OPC004"),col=c("black"),lty=1,lwd=0.5) legend("topleft",ids,col=rainbow(length(ids)),lty=1,lwd=0.5) dev.off() message("pm25 done") ################################################################### # PLOT TIME SERIES RHUM png("oklab_rhum_exclim.png",width=1900, height=600, units="px", res=300, pointsize=3) tmin <- min(aggrhum$time) tmax <- max(aggrhum$time) #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 plot(aggrhum$time,aggrhum[,2],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") #lines(bou$time,bou$pm10,col=1,t="l",lwd=0.6) #lines(kar$time,kar$pm10,col=2,t="l",lwd=0.6) #lines(koe$time,koe$pm10,col=3,t="l",lwd=0.6) #lines(lfu$time,lfu$pm10,col=4,t="l",lwd=0.6) message(paste("plotting","ref")) LWD <- 1.0 lines(asopcaggrhum$time,asopcaggrhum$id177410813,col="black",lwd=LWD) message(paste("plotting","ok")) #lines(asopcaggrhum$time,asopcaggrhum$id176430103,col="darkgrey",lwd=LWD) message(paste("plotting","ok")) for (i in seq_along(ids)) { LWD <- 0.3 if (ids[i] == "SN19014" || ids[i] == "SN19015" || ids[i] == "SN19017" || ids[i] == "SN19019" || ids[i] == "OPC-004") { LWD <- 1.0 } message(paste("plotting",ids[i])) lines(aggrhum$time,aggrhum[,i+1],xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col=rainbow(length(ids))[i],t="l",lwd=LWD) } #legend("topright",c("bou","kar","koe","lfu"),col=seq(1,4),lty=1,lwd=0.5) #legend("topright",c("id177410813","id176430103"),col=c("black","darkgrey"),lty=1,lwd=1.0) legend("topright",c("id177410813"),col=c("black"),lty=1,lwd=1.0) legend("topleft",ids,col=rainbow(length(ids)),lty=1,lwd=0.5) dev.off() message("rhum done") ################################################################### # PLOT TIME SERIES TEMP png("oklab_temp_exclim.png",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")) #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) refmax <- max(asopcaggtemp$id177410813, na.rm=TRUE) refmin <- min(asopcaggtemp$id177410813, na.rm=TRUE) allmax <- ceiling(max(aggrmax,refmax,na.rm=TRUE)) #min(max(aggrmax,luebmax),70) allmin <- floor(min(aggrmin,refmin,na.rm=TRUE)) #min(max(aggrmax,luebmax),70) plot(aggtemp$time,aggtemp[,2],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") abline(v= as.POSIXct(ihours) ,lwd=0.2,col="gray") abline(h=seq(allmin,allmax,by=1),lwd=0.2,col="gray") #lines(bou$time,bou$pm10,col=1,t="l",lwd=0.6) #lines(kar$time,kar$pm10,col=2,t="l",lwd=0.6) #lines(koe$time,koe$pm10,col=3,t="l",lwd=0.6) #lines(lfu$time,lfu$pm10,col=4,t="l",lwd=0.6) LWD <- 1.0 lines(asopcaggtemp$time,asopcaggtemp$id177410813,col="black",lwd=LWD) #lines(asopcaggtemp$time,asopcaggtemp$id176430103,col="darkgrey",lwd=LWD) for (i in seq_along(ids)) { LWD <- 0.3 if (ids[i] == "SN19014" || ids[i] == "SN19015" || ids[i] == "SN19017" || ids[i] == "SN19019" || ids[i] == "OPC-004") { LWD <- 1.0 } lines(aggtemp$time,aggtemp[,i+1],xlim=c(as.POSIXct(tstart),as.POSIXct(tstop)),col=rainbow(length(ids))[i],t="l",lwd=LWD) } #legend("topright",c("id177410813","176430103"),col=c("black","darkgrey"),lty=1,lwd=1.0) legend("topright",c("id177410813"),col=c("black"),lty=1,lwd=1.0) legend("topleft",ids,col=rainbow(length(ids)),lty=1,lwd=0.5) dev.off() message("temp done") oklabaggpm10 <- aggpm10 oklabaggpm25 <- aggpm25 oklabaggtemp <- aggtemp oklabaggrhum <- aggrhum save(oklab,oklabaggpm10,oklabaggpm25,oklabaggrhum,oklabaggtemp,file="oklab.Rdata") message("writing oklab.Rdata done!")