rdmtinez/dwd_ETo vs dwd_data_ETo

## 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')

	# 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.810 - 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: .6108math.exp(17.27x/(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 + h2math.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: .6108math.exp(17.27x/(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 sbT4(.34 - .14ea.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(J2*math.pi/365).apply(math.cos))
	#dr = 1+.033math.cos(2math.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 ((1260Gscdr)((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= 2math.pi(J-81)/364

	#seasonal correction factor (hourly) Series
	Sc = .1645((2b).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.pit1/24) , w0-(math.pit1/24), w0, Sc


	def lild(J):
	#little delta, or solar declination

	#return .409math.sin((J2*math.pi/365)-1.39)
	return .409((J2*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: .6108math.exp(17.27x/(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 TkU2es_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')