dwd_ETo vs dwd_data_ETo

# coding: utf-8

# This uses DWD data (Temp, Humidity, Wind, Solar) to Calculate ETo
# Then compares it against the ETo given by DWD. The hope is to check
# whether the OwnData_ETo script is sufficiently good.

# In[122]:

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/Desktop/dwdata_ETo/'
os.chdir(wkdir)

txt_list = [i for i in os.listdir(wkdir) if 'derived' not in i]
dwd_eto = [i for i in os.listdir(wkdir) if 'derived' in i]

for txt in txt_list:
    if 'pressure' in txt:
        with open(txt) as f:
            'in pressure'
            df = pd.read_table(f, delimiter=';')
            
            df['DateTime'] = pd.to_datetime(df['MESS_DATUM'], format='%Y%m%d%H')
            
            df.drop(['MESS_DATUM', 'STATIONS_ID', 'QN_8', 'eor'], axis=1, inplace=True)
            
            df = df[df['DateTime'] > '2017']
  
            #hectopascals to kilopascals conversion
            df['P'] = df['P']/10
            df['P0'] = df['P0']/10
            
            df.set_index('DateTime', inplace=True)
            
            'press df'
            df1=df
            
    if 'wind' in txt:
        'in wind'
        with open(txt) as f:
            df = pd.read_table(f, delimiter=';')

            df['DateTime'] = pd.to_datetime(df['MESS_DATUM'], format='%Y%m%d%H')
            
            df.drop(['MESS_DATUM', 'STATIONS_ID', 'QN_3', 'D', 'eor'], axis=1, inplace=True)
            
            df = df[df['DateTime'] > '2017']
            
            #conversion of windspeed from measured height to 2m above ground
            df['F'] = df['F']*(4.87/math.log(67.8*10 - 5.43))
            
            df.rename(columns={'F':'U2'}, inplace=True)                               
            
            df.set_index('DateTime', inplace=True)
                
            df2=df
            
    if 'temperature' in txt:
        'in temperature'
        with open(txt) as f:
            df = pd.read_table(f, delimiter=';')
            
            df['DateTime'] = pd.to_datetime(df['MESS_DATUM'], format='%Y%m%d%H')
            
            df.drop(['MESS_DATUM', 'STATIONS_ID', 'QN_9', 'eor'], axis=1, inplace=True)
            
            df = df[df['DateTime'] > '2017']
            
            df.rename(columns={'TT_TU':'Tc', 'RF_TU':'H%'}, inplace=True)
            
            df.set_index('DateTime', inplace=True)
                        
            df3 = df
            
    if 'solar' in txt:
        'in solar'
        with open(txt) as f:
            df = pd.read_table(f, delimiter=';')
            
            df['DateTime'] = pd.to_datetime(df['MESS_DATUM_WOZ'], format='%Y%m%d%H:%M')
            
            df.drop(['MESS_DATUM', 'STATIONS_ID', 'QN_592',
                     'MESS_DATUM_WOZ', 'SD_LBERG', 'eor'], axis=1, inplace=True)
            
            df=df[df['DateTime'] > '2017']
            
            df['FD_LBERG'] = df['FD_LBERG']/100
            df['FG_LBERG'] = df['FG_LBERG']/100
            
            df.rename(columns={'ATMO_LBERG':'LW_in', 'FD_LBERG':'Rso',
                               'FG_LBERG':'Rs', 'ZENIT':'Zenith'}, inplace=True)
            
            df.set_index('DateTime', inplace=True)
            
            df4=df

mdf = df4.merge(df3, how='outer', left_index=True, right_index=True)
mdf = mdf.merge(df2, how='outer', left_index=True, right_index=True)
main = mdf.merge(df1, how='outer', left_index=True, right_index=True)

'mdf final'
main

for txt in dwd_eto:
    with open(txt) as f:
        df = pd.read_table(f, delimiter=';')
        
        df=df.loc[:, 'Datum':'VGSL']
        
        df.rename(columns={'VGSL':'DWD_ETo', 'Datum':'DateTime'}, inplace=True)
        
        df['DateTime'] = pd.to_datetime(df['DateTime'], format='%Y%m%d')
        
        df.set_index('DateTime', inplace=True)

        dfETo = df

'dfETo'
dfETo[:10]
        

###DWD ETp Data is given per daily
            


# In[125]:


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


mdf = main
mdf.dropna(axis=0, subset=['Rs', 'Rso'], inplace=True)
mdf.reset_index(inplace=True)


###create lng and lat series
mdf['Lng'] = 9.2
mdf['Lat'] = 48.83
delta = Delta(mdf['Tc'])
rn_minus_g = Rn_minus_G(mdf['Tc'], mdf['H%'], mdf['DateTime'], mdf['Rs'], mdf['Lng'], mdf['Lat'])
gamma = Gamma(mdf['P0'])
es_less_ea = es_minus_ea(mdf['Tc'], mdf['H%'])
Tk = 37/(mdf['Tc'] + 273.16)
U2 = mdf['U2']

mdf['ETo'] = ((.408 * delta * rn_minus_g) + (gamma *Tk*U2*es_less_ea)) / (delta + gamma*(1+(.34*U2)))

mdf.set_index('DateTime', inplace=True)
mdf=mdf.resample('D').agg({'Rs':'sum', 'Rso':'sum', 'ETo':'sum'})


# In[ ]:




# In[129]:

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

ETo_C = mdf.merge(dfETo, how='outer', left_index=True, right_index=True)
ETo_C.reset_index(inplace=True)

ETo_C[:10]
ETo_C = ETo_C[ETo_C['ETo'] < 6]

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

p1 = plt.subplot(grid[0,0])
p2 = plt.subplot(grid[0,0])


sns.regplot(x=ETo_C['DWD_ETo'], y=ETo_C['ETo'], ax=p1)
p1.set_title('Given ETo vs. Calcualted ETo Using DWD Data for Stuttgart')
p1.set_ylabel('Cal. ETo (mm H2O')
p1.set_xlabel('Given ETo (mm H20)')

#sns.pairplot(ETo_C,  x_vars='DateTime', y_vars=['ETo', 'DWD_ETo'])

#ETo_C.plot(x='DateTime', y=['ETo', 'DWD_ETo'], ax=p2)

plt.show()
plt.close('all')