########################Wind lidar########################################
#C.Burckhardt
#2021-03
###########################################################################

####packages####
library(readr)
library(plotrix)
library(openair)

####direction####
wd<-"Z:/burc/Data/WindLidar"
setwd(wd)

####data####
#vertical wind speed
WP <- read_delim("PAY/CNR202006_19_20.csv", 
                 ";", escape_double = FALSE, col_names = FALSE, 
                 na = "NA", trim_ws = TRUE)
#reduce to 1000m and exclude every 5 measurements (vertical)
WP<-WP[,c(1:9)]
WP<-as.data.frame(WP)
WP<-WP[seq(5, nrow(WP),5),]
#time
WPtime <- read_csv("PAY/Time202006_19_20.csv", 
                  col_names = FALSE, col_types = cols(X1 = col_datetime(format = "%d-%b-%Y %H:%M:%S")))
WPtime<-WPtime[seq(5, nrow(WPtime),5),]
WPtime <- as.POSIXct(WPtime$X1,format = "%d-%b-%Y %H:%M:%S",tz="UTC")

#height
WPheight <- read_csv("PAY/Z202006_19_20.csv", 
                     col_names = FALSE)
WPheight <- WPheight$X1[1:9]

#prepare matrix
colnames(WP)<-WPheight
rownames(WP)<-WPtime
WP<- WP[!is.na(WP[,1]),]
#new time
WPtime <- as.POSIXct(rownames(WP),tz="UTC")
# get PBL height
WP_ana <- t(WP)
WP_dif<-matrix(NA,nrow = dim(WP_ana)[1]-1, ncol = dim(WP_ana)[2])
rownames(WP_dif)<- as.numeric(colnames(WP))[1:8]+50
WP_PBL <- data.frame(date=character(), PBLH=numeric())
for (i in 1:dim(WP_ana)[2]) {
  WP_dif[,i]<-diff(WP_ana[,i])
  n<-which(grepl(max(abs(WP_dif[,i])),WP_dif[,i]))[1]
  WP_PBL[i,1]<-WPtime[i]
  WP_PBL[i,2]<-as.numeric(rownames(WP_dif))[n]
}

WP_PBL<-data.frame(date=as.POSIXct(WPtime),PBLH=as.numeric(colnames(WP)[max.col(WP,ties.method = "first")]))
WP_PBL <- timeAverage(WP_PBL,avg.time = "15 min",statistic = 
                        "mean",start.date = "2020-06-19 00:00:00")

write.table(WP_PBL,"Z:/burc/WindLidar/WL_PBL.csv",sep=";",row.names = F)

#prepare plot variables
xtcks <- pretty(WPtime)
ytcks <- pretty(WPheight)
ycloc <-c(ytcks[1],max(WPheight)-c(c(35,25,15,5)*max(WPheight)/100))
ybrks<-seq( -30, 30, 5)                      
if(max(WP,na.rm=T)>ybrks[length(ybrks)]) ybrks<-c(ybrks,max(WP,na.rm=T))  
if(min(WP,na.rm=T)<ybrks[1])            ybrks<-c(min(WP,na.rm=T),ybrks)    
yclab<-ybrks; yclab[!yclab%in%seq(-30, 30, 10)]<-NA


#par(mfrow=c(1,1),oma=c(2,0,2,0),mar=c(1,5,0,5),mgp=c(2.8,0.3,0),cex=1.2,las="1",tck=-0.02)
image(x=WPtime,y=WPheight,as.matrix(WP),axes=F,col = hcl.colors(12, "Oslo", rev = F),
      ylab="",xlab="",xlim = range(xtcks),ylim=range(ytcks))
axis.POSIXct(1,at=xtcks,format = "%d %b %H:%M")
axis(2,at=ytcks,labels = T)
text(max(WPtime),ycloc[4:3],expression(paste("CNR"),paste("[ dB ]")),xpd=NA,cex=c(1),pos=4)
color.legend(max(WPtime)+3600,ycloc[1],max(WPtime)+10800,ycloc[2],yclab,rect.col = hcl.colors(12, "Oslo", rev = F),
             align="rb",gradient="y",xpd=NA,cex=1.3)
points(WP_PBL$PBLH~WP_PBL$date,pch=25,col="red",bg="white", cex=0.8)