consolidation of tables and vis

# coding: utf-8

# In[343]:

##the following will combine the datatables from DWD in Stuttgart, as
##it is the only station which collects Solar Radiation data for BW

##for the first iteration only the data from Bosch sensors in BW will


import os
import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt

from datetime import datetime as dt, timedelta

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity='all'


wkdir = 'C:/Users/MAR8RNG/BOSCH-PROJECTS/1-Humidity/python/condensed_CSVs/'
dwddir = wkdir+'dwddata/'
bscdir = wkdir+'basic/'

#directory where the dwd data is
os.chdir(dwddir)

#4928 refers to the fact that only solar radiation data from 
#station 4928 - Stuttgart was used for these calculations
txt_list = [i for i in os.listdir(dwddir) if '4928' in i]
csv_list = [i for i in os.listdir(bscdir) if 'basic' in i]

for txt in txt_list:
    
    #get pressure data - DWD provides the data in hPa, and should be converted to kPa for Penman-Monteith use
    #also note that the airpressure should be adjusted for the Height at where the farm stands
    if 'pressure' in txt:
        with open(txt) as f:
            
            "in Pressure"
            df = pd.read_table(f, delimiter=';')
            
            #convert the time into a pandas maneagable form
            df['DateTime'] = pd.to_datetime(df.MESS_DATUM, format='%Y%m%d%H')
            
            #drop irrelevant columns
            df.drop(['STATIONS_ID', 'QN_8', 'eor', 'MESS_DATUM'], axis=1, inplace=True)
            
            #take data only from 2017
            df = df[df['DateTime'] > '2017']
            
            #change units from hectoPascals into kiloPascals
            df['P0'] = df['P0']/10
            df['P'] = df['P']/10
            
            #set DateTimeIndex for merging DataFrames
            df.set_index('DateTime', inplace=True)
            
            #Main Dataframe which will hold all the DWD data
            main = df
    
    #get solar data - DWD provides the data in J/cm^2, and should be converted to MJ/m^2 * (day or hour)
    #Rs = given from the data (FG_LBERG)
    if 'solar' in txt:
        with open(txt) as f:
            
            "in solar"
            df = pd.read_table(f, delimiter=';')
            
            df['DateTime'] = pd.to_datetime(df.MESS_DATUM, format='%Y%m%d%H:%M')
            df.drop(['STATIONS_ID', 'MESS_DATUM', 'QN_592',
                     'eor', 'ATMO_LBERG', 'ZENIT', 
                     'MESS_DATUM_WOZ'], axis=1, inplace=True)
            
            df=df[df['DateTime'] > '2017']
            
            #conversion of SunLight Minutes into hours
            df['SD_LBERG'] = df['SD_LBERG'] / 60
            
            #conversion of J/cm^2 units MJ/m^2 -> (1MJ/10^6 J)*(100cm/1m)^2 = 
            # = conversion factor (1/100) MJ/m^2 * hour
            df['FD_LBERG'] = df['FD_LBERG'] / 100
            df['FG_LBERG'] = df['FG_LBERG'] / 100
            
            #original timestamp needed to be changed so that
            #time indexing and merging works properly 
            df['DateTime'] = df['DateTime'].dt.round('H')
            
            #DataFrames with DatetimeIndex are to be merged, not concatenated
            df.set_index('DateTime', inplace=True)         
            main = main.merge(df, how='outer', left_index=True, right_index=True)

    #get wind data - DWD provides the data as the mean wind speed in m/s
    if 'wind' in txt:
        with open(txt) as f:
            
            "in Wind"
            df = pd.read_table(f, delimiter=';')
            
            df['DateTime'] = pd.to_datetime(df.MESS_DATUM, format='%Y%m%d%H')
            df.drop(['STATIONS_ID', 'MESS_DATUM', 'QN_3', 'D', 'eor'], axis=1, inplace=True)            
            
            df = df[df['DateTime'] > '2017']
            
            #conversion of windspeed measured at 10m 
            #to the windspeed at two meters (m/s)
            df['F'] = df['F']*(4.87/math.log(67.8*10 - 5.42))

            df.set_index('DateTime', inplace=True)
            
            'Main df'
            main = main.merge(df, how='outer', left_index=True, right_index=True)


main.rename(columns={'FD_LBERG':'Rso', 
           'FG_LBERG':'Rs', 'SD_LBERG':'SunHour', 'F':'U2'}, inplace=True)

main[3000:3005]


#temp, humidity can be obtained from own data - OWN data should be resampled PER HOUR first


# In[471]:

os.chdir(bscdir)
import math

def Delta(T):
    """Delta slope calculation"""
    eT = lambda x: .6108*math.exp(17.27*x/(x+237.3))

    top=4098*(T.apply(eT))
    bot=(237.3+T).pow(2)

    return top/bot


def Rn_minus_G(T, HM1, time, Rs, Lm, Lat):
    
    """Net Solar radiation minus Soil heatflux Calculation"""
    
    h = 315 #the station-height of the solar measurement, this is the station-height for Stuttgart
    
    #Net Shortwave and Longwave radiation
    # Rns(1 - a)*Rs, .23 = reference crop albedo
    Rns = (1-.23)*Rs
    
    Ra = sol_Ra(Lm, Lat, time)
    
    #This changes the negative values given by the calculations
    #during nighttime hours, at night there is Zero (0) solar radiation
    Ra=Ra.where(Ra.values > 0, 0)
    
    Rso = sol_Rso(h, Ra)
    
    
    #concatenate Rs Rso
    
    
    
    Rnl = sol_Rnl(T,HM1,Rs,Rso)
    
    #Net Radiation
    Rn = Rns - Rnl
    
    
    
    #Soil heatflux hourly calculation
    G=Ghr(time, Rn)
   
    print('Rs, Rns, Ra, Rso, Rnl, Rn, G')
    RSC = pd.concat([Rs, Rns, Ra, Rso, Rnl, Rn, G], axis=1, verify_integrity=True)
    print(RSC[1816:1817])
    
    return Rn - G


def sol_Rso(h, Ra):
    """Rso: Rs on a cloudless day calculator"""
    return (.75 + h*2*math.pow(10, -5))*Ra

def sol_Rnl(T,HM1,Rs,Rso):
    """Net Longwave Calculation"""

    #sb, Stefan-Boltzmann hourly conversion (/24) MJ K^-4 * m^-2 * hr^-1
    sb=2.043*math.pow(10, -10)
    T4 = (273.16+T).pow(4) #Temp in Kelvin raised to the fourth power

    #vapor pressure conversions
    eT = lambda x: .6108*math.exp(17.27*x/(x+237.3))
    ea = T.apply(eT)*(HM1/100)
    
    Rsso=Rso
    #correction of Rso values less than 0,
    #since Rso is dependend on Ra and Ra is 0 at night
    Rso = Rso.where(Rs.values < Rso.values, Rs.values)
    Rso = Rso.where(Rso.values>0, 1)
    
    print('Rso < Rs values to Rs, 0 to 1: Rs, RsO, Rso, Rs/Rso')
    rso3 = pd.concat([Rs, Rsso, Rso, Rs/Rso], axis=1)
    print(rso3[1816:1840])
    
    return sb*T4*(.34 - .14*ea.apply(np.sqrt))*(1.35*(Rs/Rso) - .35)

def sol_Ra(Lm, Lat, time):
   
    J = time.dt.dayofyear
    
    #Lz--the Longitude of the Center of the TimeZone (i.e. CEST) 
    #where the measurements were taken in deegres West of Greenwich(DWG)
    Lz = pd.Series(data= 360 - 15.35, index=Lm.index)
   
    #Lm--the Longitude of the measurement site in DWG, Series
    Lm = 360-Lm
    #print('Lm')
    #print(Lm)
    #J = Transformation of the time Series into Julian day
    
    #Gsc .0820 MJ * m^-2 * min-1
    Gsc= .0820
    
    #dr = inverse relative distance Earth-Sun
    
    dr = 1+(.033*(J*2*math.pi/365).apply(math.cos))
    #dr = 1+.033*math.cos(2*math.pi*J/365)
    
    #solar time angles at end and beginning of period (Radian)
    w2, w1, w0, Sc = omegas(time, J, Lz, Lm)
    
    #print('Lm', 'Lat','Lz', 'w2', 'w1', 'w0', 'Sc')
    #oms = pd.concat([Lm, Lat, Lz, w2, w1, w0, Sc], axis=1, names=['w2', 'w1', 'w0', 'Sc'])
    #print(oms[:100])
    #lat and long in Radians
    #ld -- solar declination  dependend only on the day of the year
    #lp -- the latitude (Ln) of the measurements in Radians
    ld = lild(J)
    lp = lilp(Lat)

    
    return ((12*60*Gsc*dr)*((w2-w1)*lp.apply(math.sin)*ld.apply(math.sin)+
                           lp.apply(math.cos)*ld.apply(math.cos)*
                           (w2.apply(math.sin) - w1.apply(math.sin))))/math.pi


def omegas(time, J, Lz, Lm):
    
    #helper function
    def time_to_decimal(t):
        "helper function timestamp to decimal time conversion"
        h, m = [int(i) for i in t.split(':')]
        return h + m/60
        
    #t -- the conversion of the DateTime series to hour (%H) format

    #Since the time is the stamp of the end of data collection period, this calculation subtracts .5 hours
    t0=time - timedelta(hours=.5) #midpioint of the hour (hour - .5 i.e. 13-to-14, => 14:00 - .5 = 13:30)
    
    #variable holding the decimal converted timestamps
    t0=t0.dt.strftime('%H:%M').apply(time_to_decimal)
    
    t1=1 #one hour time step

    #Seasonal correction factor Series
    b= 2*math.pi*(J-81)/364

    #seasonal correction factor (hourly) Series
    Sc = .1645*((2*b).apply(math.sin)) - .1255*(b.apply(math.cos)) -.025*(b.apply(math.sin))

    w0 = math.pi*((t0 + .06667*(Lz-Lm) + Sc) -12)/12

    return w0+(math.pi*t1/24) , w0-(math.pi*t1/24), w0, Sc
    
    
def lild(J):
    #little delta, or solar declination
    
    #return .409*math.sin((J*2*math.pi/365)-1.39)
    return .409*((J*2*math.pi/365)-1.39).apply(math.sin)
    


def lilp(Lat):
    #little phi, conversion of latitude into Radians
    return Lat*math.pi/180

def Ghr(time, Rn):
    """Soil Heat flux calculator"""
    
    #sunrise -sunset hours from 1st of the month (sr,ss) to end of the month
    #Jan-1st SR 8:17, Jan-31st SR 07:47
    #Jan-1st SS 16:03, Jan-31st SS 16:52  
    daylight = {'Jan':('08:00', '16:25'),
            'Feb':('07:20', '17:20'),
            'Mar':('06:45', '18:40'),
            'Apr':('06:05', '20:05'),
            'May':('05:10', '20:55'),
            'Jun':('04:50', '21:25'),
            'Jul':('05:10', '21:15'),
            'Aug':('05:50', '20:30'),
            'Sep':('06:45', '19:25'),
            'Oct':('07:00', '17:30'),
            'Nov':('07:25', '16:15'),
            'Dec':('08:05', '16:00')}

    #this maps the month of 'time' as to the dictionary daylight
    #then unzips the tuple to daystart and dayend lists which 
    #are then converted to pandas.Series objects for use in boolean mask
    daystart, dayend = zip(*time.dt.strftime('%b').map(daylight))
    
    #conversion to pandas series
    daystart=pd.Series(pd.to_datetime(daystart).time)
    dayend=pd.Series(pd.to_datetime(dayend).time)
    
    t=pd.Series(time.dt.time)
    light_mask = t.between(daystart, dayend)
    

    #True values get multiplied by .1 (The inversion is because 
    #pd.Series.where() only changes the values which are negative)
    Ghr = Rn.where(light_mask, Rn*.5)

    #False values multiplied by *.5
    Ghr = Ghr.where(np.invert(light_mask), Ghr*.1)
    
    #sbs = light_mask.merge(Ghr, how='outer', left_index=True, right_index=True)
    #print('Ghr')
    #sbs =  pd.concat([time, Rn, Ghr, light_mask], axis=1)
    #print(sbs[:100])
    return Ghr


def cd_get(time):
    "Cd: short crop denominator constant dependend on night vs daylight"
    
    Cd = pd.Series(data=.24, index=time.index)
    
    daylight = {'Jan':('08:00', '16:25'),
            'Feb':('07:20', '17:20'),
            'Mar':('06:45', '18:40'),
            'Apr':('06:05', '20:05'),
            'May':('05:10', '20:55'),
            'Jun':('04:50', '21:25'),
            'Jul':('05:10', '21:15'),
            'Aug':('05:50', '20:30'),
            'Sep':('06:45', '19:25'),
            'Oct':('07:00', '17:30'),
            'Nov':('07:25', '16:15'),
            'Dec':('08:05', '16:00')}

    #this maps the month of 'time' as to the dictionary daylight
    #then unzips the tuple to daystart and dayend lists which 
    #are then converted to pandas.Series objects for use in boolean mask
    daystart, dayend = zip(*time.dt.strftime('%b').map(daylight))
    
    #conversion to pandas series
    daystart=pd.Series(pd.to_datetime(daystart).time)
    dayend=pd.Series(pd.to_datetime(dayend).time)
    
    t=pd.Series(time.dt.time)
    daylight_mask = t.between(daystart, dayend)
    
    return Cd.where(daylight_mask, .96)
    
def Gamma(P0):
    """Psychrometric Constant calculation"""
    
    return .000665*P0


def es_minus_ea(T, HM1):
    """Vapor Pressure Calculation"""
    
    #vapor pressure calculation
    eT = lambda x: .6108*math.exp(17.27*x/(x+237.3))

    es = T.apply(eT)
    ea = es*(HM1/100)

    return es - ea


ETo_Inst = []    
for csv in csv_list:
    imei = csv[-8:-4]
    imei
    with open(csv) as f:
        df=pd.read_csv(f, parse_dates=['DateTime'])
        
        #drop irrelevant data
        df = df[df['Volts'] > -254]
        df.drop(['IMEI', 'Volts'], inplace=True, axis=1)
        
        #resample own data at hourly intervals to merge onto DWD data 
        df.set_index('DateTime', inplace=True)
        df=df.resample('H').agg({'DigiTemp':'mean', 'HM1':'mean',
                                 'T1':'mean', 'VWC':'mean', 'Lat':'mean', 'Lng':'mean'})
        
        #create and fill change in VWC column at hourly interval
        df.reset_index(inplace=True)
        df['dVWC'] = None
        df.loc[(df['DateTime'].shift(-1)-df['DateTime'])==timedelta(hours=1),
              'dVWC'] = df['VWC'].shift(-1) - df['VWC']
        df['dVWC'] = df['dVWC'].shift(+1)
        df = df[1:] #drops first row as dVWC is NaN
        
        #WaterLoss column for coparision with ETo
        #df['WLoss'] = None
        
        df.set_index('DateTime', inplace=True)
        
        ##MERGE WITH Main DataFrame
        
        mdf = main.merge(df, how='outer', left_index=True, right_index=True)        
        mdf.reset_index(inplace=True)
        
        mdf.dropna(axis=0, subset=['T1', 'Rs', 'Rso'], inplace=True)
        mdf.reset_index(inplace=True)
        
        mdf['ETo'] = None
       
        #the vectorized individual parts of the Penman-Monteith Equation
        #321\        
        
        delta = Delta(mdf['T1'])

        rn_minus_g = Rn_minus_G(mdf['T1'], mdf['HM1'], mdf['DateTime'], mdf['Rs'], mdf['Lng'], mdf['Lat'])

        gamma = Gamma(mdf['P0'])

        es_less_ea = es_minus_ea(mdf['T1'], mdf['HM1'])

        Cd = cd_get(mdf['DateTime'])
        
        Tk = 37/(mdf['T1'] + 273.16)
        U2 = mdf['U2']
        
        mdf['ETo'] = ((.408 * delta * rn_minus_g) + (gamma *Tk*U2*es_less_ea)) / (delta + gamma*(1+(Cd*U2)))
        
        
        ETo_Inst.append((imei,mdf))
        


# In[462]:

mains = ETo_Inst
for imei, df in mains:
    df=df.drop('index', axis=1)
    
    df=df[df['DateTime'] > 'June-21-2017']
    df=df[df['DateTime' ]< 'June-22-2017']
    df


# In[ ]:




# In[472]:

import numpy as np
import seaborn as sns

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


from matplotlib import pyplot as plt, gridspec
from matplotlib.backends.backend_pdf import PdfPages

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

thisdir='C:/Users/MAR8RNG/BOSCH-PROJECTS/1-Humidity/python/condensed_CSVs/basic/'
thatdir='C:/Users/MAR8RNG/BOSCH-PROJECTS/1-Humidity/python/condensed_CSVs/'
os.chdir(thisdir)

   
dwd_ETo = []
txt_list = [i for i in os.listdir(thisdir) if 'dwd' in i]
for txt in txt_list:
    imei = txt[-8:-4]
    imei
    with open(txt) as f:
        df = pd.read_table(f, parse_dates=['Datum'],delimiter=';')
        df = df.loc[:, 'Datum':'VGSL']
        df.rename(columns={'VGSL':'DWD_ETo'}, inplace=True)
        df['Datum'] = pd.to_datetime(df['Datum'], format='%Y%m%d%H')
        df.set_index('Datum', inplace=True)
        
        dwd_ETo.append(df)

merge_count = 0            
for imei, df in ETo_Inst:
    print(imei)
    #convert the values that are less than Zero or NaN
    #df.loc[(df['ETo'] < 0), 'ETo'] = 0
    df=df.fillna(0)
    
    #convert positive+ change in water into 0 to only account for daily water loss
    df.loc[(df['dVWC'] > 0), 'dVWC'] = 0
    
    #Renaming for later use as Legend names
    df['VWCmax'] = df['VWC']
    df['VWCmin'] = df['VWC']
    
    #set index datetime for resampling purposes
    df.set_index('DateTime', inplace=True)
    df = df.resample('D').agg({'VWC':'mean', 'VWCmax':'max',
                               'VWCmin':'min', 'dVWC':'sum',
                               'ETo':'sum',})
    
    #convert negative dVWC volumes into absolute (mm) WaterLoss values
    df['mmWaterLoss'] = -df['dVWC'] * 1000
    df['VWC'] = df['VWC'] * 1000
    
    df.rename(columns={'ETo':'Own_ETo', 'VWC':'mmWater_mean'}, inplace=True)
    
    ###MergeWith DWD_ETo
    df=df.merge(dwd_ETo[merge_count], how='outer', left_index=True, right_index=True)
    merge_count+=1
    
    
    ##Drop Rows where there are NAs for OWN_ETo has NAs since 
    ##after the merge some days will have this calculation due to unavailable data from 
    ##the beginning of the year (see original data collection)
    df.dropna(axis=0, subset=['Own_ETo'], inplace=True)
    
    df.reset_index(inplace=True)
    
    'resample daily min, max'
    df['Own_ETo'].idxmin()
    df['Own_ETo'].idxmax()
    
    df.rename(columns={'index':'DateTime'}, inplace=True)
    #df['DateTime'] = pd.to_datetime(df['DateTime'])
    
    
    ########################################################################
    #df = df[df['mmWaterLoss'] < 20] ########################################
    ########################################################################
    
    ########################################################################
    df = df[df['mmWaterLoss'] < 150] ########################################
    ########################################################################
    
    ########################################################################
    #df = df[df['mmWaterLoss'] < 50] ########################################
    ########################################################################
    
    with PdfPages(thatdir+'figures/'+'OWN_ETo_'+imei+'.pdf') as pdf:

        fig = plt.figure(figsize=(16,36))
        grid = gridspec.GridSpec(6,1)

        
        """p3d = fig.add_subplot(111, projection='3d')
        p3d.scatter(xs=df['mmWaterLoss'], ys=df['Own_ETo'], zs=df['DateTime'], zdir='x')
        plt.show()
        #p3D = plt.subplot(grid[])
        break"""
        p1 = plt.subplot(grid[0,0]) #WaterLoss vs oETo
        p2 = plt.subplot(grid[1,0]) #WaterLoss vs dETo
        p3 = plt.subplot(grid[2,0]) #time vs. mean VWC, WaterLoss
        p4 = plt.subplot(grid[3,0]) #oETo vs. dETo
        p5 = plt.subplot(grid[4,0]) #time vs. oETo vs. dETo
        
        
        
        fz= 20
        #WaterLoss vs OWN_ETo
        #sns.jointplot(x=df['mmWaterLoss'], y=df['Own_ETo'])
        sns.regplot(x=df['mmWaterLoss'], y=df['Own_ETo'], ax=p1)
        #df.plot.scatter(x='mmWaterLoss', y='Own_ETo', ax=p1, color='Green', fontsize=fz)
        p1.set_title('WaterLoss vs Own_ETo', fontsize=fz)
        p1.set_ylabel('Own ETo (mm/day) ' + '{}'.format(imei), fontsize=fz)
        p1.set_xlabel('mm WaterLoss /day', fontsize=fz)
        
        #WaterLoss vs DWD_ETo
        sns.regplot(x=df['mmWaterLoss'], y=df['DWD_ETo'], ax=p2)
        #df.plot.scatter(x='mmWaterLoss', y='DWD_ETo', color='Blue', ax=p2, fontsize=fz)
        p2.set_title('WaterLoss vs. DWD_ETo', fontsize=fz)
        p2.set_ylabel('DWD ETo (mm/day) ' + '{}'.format(imei), fontsize=fz)
        p2.set_xlabel('mm WaterLoss /day', fontsize=fz)
        
        #WaterLoss and Mean VWC(mm) vs Time
        df.plot(x='DateTime', y=['mmWater_mean','mmWaterLoss'], ax=p3, fontsize=fz, color=['orange', 'blue'])
        p3.set_title('Water Loss & Content vs Time', fontsize=fz)
        p3.set_xlabel('DateTime', fontsize=fz)
        p3.set_ylabel('Loss & Content /day '+ '{}'.format(imei), fontsize=fz)
        p3.legend(fontsize=fz)
        
        #own ETo vs DWD ETo
        sns.regplot(x=df['Own_ETo'], y=df['DWD_ETo'], ax=p4)
        #df.plot.scatter(x='Own_ETo', y='DWD_ETo', color='red', ax=p4, fontsize=fz)
        p4.set_title('Calculated Own Data ETo vs Deutsche Wetterdienst ETo', fontsize=fz)
        p4.set_xlabel('Own ETo (mm) ', fontsize=fz)
        p4.set_ylabel('DWD ETo (mm) ' + '{}'.format(imei), fontsize=fz)
    
        #time vs dETo vs. DWD ETo
        df.plot(x='DateTime', y=['DWD_ETo', 'Own_ETo'], ax=p5, fontsize=fz, color=['blue', 'green'])
        #style can be set but color must be removed, very confusing to read
        p5.set_title('Time vs Own_ETo and DWD_ETo ')
        p5.set_xlabel('DateTime', fontsize=fz)
        p5.set_ylabel('ETo (mm) ' +'{}'.format(imei), fontsize=fz)
        p5.legend(fontsize=fz)
        
        fig.tight_layout()
        plt.show()

        pdf.savefig(fig)
        plt.close('all')

    ###OMIT THE VALUES THAT ARE GREATER THAN 50 and compare....
    ###You have to consider that on certain days after an 
    #irrigation event the water loss will be greater given
    #that percolation is by far a bigger quantity and happens at a 
    #faster rate than evaporation

    #df.plot(x=df['DateTime'], y=['ETo'])
    #plt.show()


# In[ ]: