# -*- coding: utf-8 -*-
"""
analysis tools to compare mast and ZephIR wind data

read data, interpolation and combining datasets, to prepare datasets for further analysis

@author: Steven Knoop (steven.knoop@knmi.nl)
"""

from datetime import datetime,timedelta
import numpy as np
import pandas as pd
import datetime as dt
from scipy.interpolate import interp1d
import os

################################################################################################################################################
#tools
################################################################################################################################################

def findheights(df):
    listcolumns=df.columns.tolist()
    hlist=[s[1:] for s in listcolumns if 'u' in s]
    hlist.sort(key=int)
    
    return hlist

def findmatchingheights(h1,h2):
    "find common heights and sort"
    hlist=list(set(h1).intersection(h2))
    hlist.sort(key=int)
    
    return hlist

def datespan(startDate, endDate, delta=timedelta(days=1)):
    currentDate = startDate
    while currentDate <= endDate:
        yield currentDate
        currentDate += delta

def createemptydataframe(df):
    "create a dataframe at 10-min frequency to fill in gaps of sensordata"
    
    startdate=df.index[0].date()
    enddate=df.index[-1].date()  
    
    cols=df.columns.tolist()
    
    rng = pd.date_range(startdate,enddate, freq='10T')
    ts=pd.Series(np.nan,index=rng)
    tsempty=pd.concat([ts]*len(cols),axis=1)
    tsempty.columns = cols
    
    return tsempty

def heightinterpolation(dgo):
    """
    produce additional 11 and 39 m data from mast, based on interpolation

    scalar averaging!!!
    """
    df=dgo.copy()
    
    hlist1=findheights(df)
    hlist2=[11,39]

    hlist1=np.array([int(i) for i in hlist1])
    hlist2=np.array([int(i) for i in hlist2])
    
    for h1 in hlist1:
        df['d'+str(h1)+'X']=np.sin(np.radians(df['d'+str(h1)]))         #decompose wind direction on unit circle
        df['d'+str(h1)+'Y']=np.cos(np.radians(df['d'+str(h1)]))

    dg = pd.DataFrame(columns=['u'+str(int(i)) for i in hlist2]+['d'+str(int(i)) for i in hlist2], index=df.index)
  
    for index, row in df.iterrows():                    #iterate per row
        print(index)
        yu=[]
        ydX=[]
        ydY=[]
        for h1 in hlist1:  
            yu.append(row['u'+str(h1)])
            ydX.append(row['d'+str(h1)+'X'])
            ydY.append(row['d'+str(h1)+'Y'])
            
        fu=interp1d(hlist1, yu, kind='cubic')          #linear interpolation function
        fdX=interp1d(hlist1, ydX, kind='cubic')
        fdY=interp1d(hlist1, ydY, kind='cubic')
       
        for h2 in hlist2:
            dg['u'+str(h2)].loc[index]=fu(h2)               #find values of wind lidar heights
            dummyd=np.degrees(np.arctan(fdX(h2)/fdY(h2)))   #transform back to wind direction
            dg['d'+str(h2)].loc[index]=np.where(fdY(h2)<0,dummyd+180,dummyd) % 360
             
    dg2=pd.concat([df,dg],axis=1)
    
    del df, dg

    return dg2  

################################################################################################################################################
#apply height interpolated or scaling
################################################################################################################################################

def writeTowerdata_interpolated():

    dgo=readTowerdata_csv()

    dg=heightinterpolation(dgo)
    
    startdate=dg.index.to_numpy()[0]
    enddate=dg.index.to_numpy()[-1]
            
    subfilename=str(startdate)[:10].replace('-','')+'_'+str(enddate)[:10].replace('-','')

    dg.to_csv('Towerdata_Cabauw_interpolated_'+subfilename+'.csv')  
    
 
def write_AWS_data_scaled_11m():
    """
    scale 10m AWS data to 11m, using the 11m/10m ratio from the interpolated Tower data for wind speed
    and adding the difference between 10 and 11m wind direction of the Tower data to the AWS 10m
    """

    dg0=readTowerdata_interpolated_csv()
    dg1=readAWSwinddata_csv()
    
    startdate=max(dg0.index.to_numpy()[0],dg1.index.to_numpy()[0])
    enddate=min(dg0.index.to_numpy()[-1],dg1.index.to_numpy()[-1])

    dg0=dg0[startdate:enddate]
    dg1=dg1[startdate:enddate]

    dg=pd.DataFrame()
 
    """
    wind speed
    """
    dg['u11']=dg1['u10']*dg0['u11']/dg0['u10']
    dg['u11'][dg0['u10']==0]=np.nan
    
    """
    wind direction
    """
    dg['d11']=dg1['d10']+(dg0['d11']-dg0['d10'])
    dg['d11'][dg['d11']>360]=dg['d11']-360
    dg['d11'][dg['d11']<0]=dg['d11']+360
   
    startdate=dg.index.to_numpy()[0]
    enddate=dg.index.to_numpy()[-1]
    
    subfilename=str(startdate)[:10].replace('-','')+'_'+str(enddate)[:10].replace('-','')
    
    dg.to_csv('AWSwind_scaled_11m_'+subfilename+'.csv')      
    
################################################################################################################################################
#read data from the  csv files
################################################################################################################################################

def readZephIRdata_csv(ext=False):

    if ext:
        filename='ZephIRdata_Cabauw_ext_2019'
    else:
        filename='ZephIRdata_Cabauw_2018'        
    
    for file in os.listdir():
        if '.csv' in file:
            if filename in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    return df

def readAWSwinddata_csv():   

    for file in os.listdir():
        if '.csv' in file:
            if 'AWSwind_20' in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    return df


def readTowerdata_csv():

    for file in os.listdir():
        if '.csv' in file:
            if 'Towerdata_Cabauw_20' in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    return df

def readTowerVisdata_csv():

    for file in os.listdir():
        if '.csv' in file:
            if 'TowerVisdata_Cabauw_20' in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    return df

def readSurfacedata_csv():

    for file in os.listdir():
        if '.csv' in file:
            if 'Surfacedata_Cabauw_20' in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    df = df[~df.index.duplicated()] 

    return df

def read_ceilometer_data_csv():

    for file in os.listdir():
        if '.csv' in file:
            if 'ceilometer_C1_Cabauw_20' in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    return df 

def readTowerdata_interpolated_csv():

    for file in os.listdir():
        if '.csv' in file:
            if 'Towerdata_Cabauw_interpolated_20' in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    return df

def readAWSdata_scaled_11m_csv():

    for file in os.listdir():
        if '.csv' in file:
            if 'AWSwind_scaled_11m_20' in file: 
                df=pd.read_csv(file,index_col=0, parse_dates=True)

    return df

def readTowerdata_interpolated_AWSdata_scaled_csv():
    """
    combine Tower data (interpolated to get 39m) with AWS data (scaled to get 11m, and original 10m for wind sector selection)
    drop 10m and 11m data from Tower data
    """

    df1=readTowerdata_interpolated_csv()
    df2=readAWSdata_scaled_11m_csv()
    df3=readAWSwinddata_csv()

    df1=df1.drop(['d11', 'u11', 'd10', 'u10'], axis=1)

    startdate=max(df1.index.to_numpy()[0],df2.index.to_numpy()[0])
    enddate=min(df1.index.to_numpy()[-1],df2.index.to_numpy()[-1])

    df=pd.concat([df1,df2,df3],axis=1)
    df=df[startdate:enddate]
    
    return df

################################################################################################################################################
#prepare data for DQ analysis with proper wind sector or other applications
################################################################################################################################################

start=dt.datetime.combine(dt.date(2018,2,15),dt.time(0,0,0))   
end=dt.datetime.combine(dt.date(2020,2,29),dt.time(23,59,59))

def ComparisonDataPreparationWSW():

    """
    take only the free wind sector
    select on wind directon lowest height (10m)
    """
    
    df=readZephIRdata_csv()
    dg=readTowerdata_interpolated_AWSdata_scaled_csv()
  
    df=df[start:end] 
    dg=dg[start:end] 
            
    hlist=findmatchingheights(findheights(df),findheights(dg))
    
    startdate=max(df.index.to_numpy()[0],dg.index.to_numpy()[0])
    enddate=min(df.index.to_numpy()[-1],dg.index.to_numpy()[-1])

    df=df[startdate:enddate]
    dg=dg[startdate:enddate]

    cols=['u','d']

    crit=(dg['d10']>200) & (dg['d10']<250)
    
    df1=pd.DataFrame()
    dg1=pd.DataFrame() 
    
    for col in cols:
        for h in hlist:
            df1[col+h]=df[col+h].copy()
            dg1[col+h]=dg[col+h].copy()

            df1[col+h][~crit]=np.nan  
            dg1[col+h][~crit]=np.nan
            
    subfilename=str(startdate)[:10].replace('-','')+'_'+str(enddate)[:10].replace('-','')+"_WSW_AWS_scaled"

    return df1,dg1,hlist,subfilename

if __name__=='__main__':

    ComparisonDataPreparationWSW()