%matplotlib inline
import sys
import os
#sys.path.append('..') #append the parent directory to the system path
import pvlib
import pandas as pd
import matplotlib.pyplot as plt
try:
import seaborn
seaborn.set_context('paper')
except ImportError:
pass
datapath = os.path.join(pvlib.__path__[0], 'test', 'data')
fname = '703165TY.csv' #Use absolute path if the file is not in the local directory
TMY, meta = pvlib.tmy.readtmy3(filename=os.path.join(datapath, fname))
TMY.index.name = 'Time'
meta['SurfTilt'] = 30
meta['SurfAz'] = 0
meta['Albedo'] = 0.2
# create pvlib Location object
sand_point = pvlib.location.Location(meta['latitude'], meta['longitude'], tz='US/Alaska',
altitude=meta['altitude'], name=meta['Name'].replace('"',''))
print(sand_point)
# TMY data seems to be given as hourly data with time stamp at the end
# shift the index 30 Minutes back for calculation of sun postions
TMY = TMY.shift(freq='-30Min')
SAND POINT: latitude=55.317, longitude=-160.517, tz=US/Alaska, altitude=7.0
Calculate the solar position for all times in the TMY file. Requires either PyEphem or a compiled spa_c library.
solpos = pvlib.solarposition.get_solarposition(TMY.index, sand_point, method='pyephem')
solpos.head()
apparent_elevation | apparent_azimuth | elevation | azimuth | apparent_zenith | zenith | |
---|---|---|---|---|---|---|
Time | ||||||
1997-01-01 00:30:00-09:00 | -54.750103 | 328.861096 | -54.750103 | 328.861096 | 144.750103 | 144.750103 |
1997-01-01 01:30:00-09:00 | -57.530205 | 353.265280 | -57.530205 | 353.265280 | 147.530205 | 147.530205 |
1997-01-01 02:30:00-09:00 | -56.627712 | 18.751764 | -56.627712 | 18.751764 | 146.627712 | 146.627712 |
1997-01-01 03:30:00-09:00 | -52.329549 | 41.456666 | -52.329549 | 41.456666 | 142.329549 | 142.329549 |
1997-01-01 04:30:00-09:00 | -45.717933 | 60.026503 | -45.717933 | 60.026503 | 135.717933 | 135.717933 |
# append solar position data to TMY DataFrame
TMY['SunAz'] = solpos['apparent_elevation']
TMY['SunZen'] = solpos['apparent_zenith']
TMY['HExtra'] = pvlib.irradiance.extraradiation(doy=TMY.index.dayofyear)
TMY['AM'] = pvlib.atmosphere.relativeairmass(z=TMY.SunZen)
DFOut = pvlib.clearsky.disc(Time=TMY.index,GHI=TMY.GHI, SunZen=TMY.SunZen)
TMY['DNI_gen_DISC'] = DFOut['DNI_gen_DISC']
TMY['Kt_gen_DISC'] = DFOut['Kt_gen_DISC']
TMY['AM'] = DFOut['AM']
TMY['Ztemp'] = DFOut['Ztemp']
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-4-0523b848c9cd> in <module>() 5 TMY['HExtra'] = pvlib.irradiance.extraradiation(doy=TMY.index.dayofyear) 6 ----> 7 TMY['AM'] = pvlib.atmosphere.relativeairmass(z=TMY.SunZen) 8 9 DFOut = pvlib.clearsky.disc(Time=TMY.index,GHI=TMY.GHI, SunZen=TMY.SunZen) /home/will/git_repos/pythonic-PVLIB/pvlib/atmosphere.pyc in relativeairmass(z, model) 227 Expect={'z': ('array','num','x<=90','x>=0'), 228 'model': ('default','default=kastenyoung1989')} --> 229 var=pvt.Parse(Vars,Expect) 230 231 if ('kastenyoung1989') == var.model.lower(): /home/will/git_repos/pythonic-PVLIB/pvlib/pvl_tools.pyc in __init__(self, dct, Expect) 78 ''' 79 def __init__(self, dct, Expect): ---> 80 self.__dict__.update(self.parse_fcn(dct,Expect)) 81 82 def parse_fcn(self, kwargs, Expect): /home/will/git_repos/pythonic-PVLIB/pvlib/pvl_tools.pyc in parse_fcn(self, kwargs, Expect) 196 raise Exception('Error: Numeric input "'+arg+' " fails on logical test " '+ re.findall(reg,string)[0]+'"') 197 except: --> 198 if not(eval(lambdastring)(kwargs[arg])): 199 raise Exception('Error: Numeric input "'+arg+' " fails on logical test " '+ re.findall(reg,string)[0]+'"') 200 ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
models = ['Perez', 'Hay-Davies', 'Isotropic', 'King', 'Klucher', 'Reindl']
TMY['Perez'] = pvlib.pvl_perez(SurfTilt=meta['SurfTilt'],
SurfAz=meta['SurfAz'],
DHI=TMY.DHI,
DNI=TMY.DNI,
HExtra=TMY.HExtra,
SunZen=TMY.SunZen,
SunAz=TMY.SunAz,
AM=TMY.AM)
TMY['Hay-Davies'] = pvlib.pvl_haydavies1980(SurfTilt=meta['SurfTilt'],
SurfAz=meta['SurfAz'],
DHI=TMY.DHI,
DNI=TMY.DNI,
HExtra=TMY.HExtra,
SunZen=TMY.SunZen,
SunAz=TMY.SunAz)
TMY['Isotropic'] = pvlib.pvl_isotropicsky(SurfTilt=meta['SurfTilt'],
DHI=TMY.DHI)
TMY['King'] = pvlib.pvl_kingdiffuse(SurfTilt=meta['SurfTilt'],
DHI=TMY.DHI,
GHI=TMY.GHI,
SunZen=TMY.SunZen)
TMY['Klucher'] = pvlib.pvl_klucher1979(SurfTilt=meta['SurfTilt'],
SurfAz=meta['SurfAz'],
DHI=TMY.DHI,
GHI=TMY.GHI,
SunZen=TMY.SunZen,
SunAz=TMY.SunAz)
TMY['Reindl'] = pvlib.pvl_reindl1990(SurfTilt=meta['SurfTilt'],
SurfAz=meta['SurfAz'],
DHI=TMY.DHI,
DNI=TMY.DNI,
GHI=TMY.GHI,
HExtra=TMY.HExtra,
SunZen=TMY.SunZen,
SunAz=TMY.SunAz)
yearly = TMY[models].resample('A', how='sum').squeeze() / 1000.0 # kWh
monthly = TMY[models].resample('M', how='sum', kind='period') / 1000.0
daily = TMY[models].resample('D', how='sum') / 1000.0
ax = TMY[models].plot(title='Modelled in-plane diffuse irradiance')
# ax.legend(ncol=3, loc='upper center', title='In-plane diffuse irradiance')
ax.set_ylim(0, 800)
ylabel = ax.set_ylabel('Diffuse Irradiance [W]')
TMY[models].describe()
Perez | Hay-Davies | Isotropic | King | Klucher | Reindl | |
---|---|---|---|---|---|---|
count | 4531.000000 | 8760.000000 | 8760.000000 | 8760.000000 | 8760.000000 | 8760.000000 |
mean | 126.020357 | 58.782821 | 56.830698 | 64.360708 | 64.622931 | 59.137119 |
std | 98.944032 | 87.587632 | 85.451560 | 91.309309 | 95.026610 | 88.152787 |
min | 0.715740 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
25% | 58.800900 | 0.000000 | 0.000000 | 0.000000 | 0.000000 | 0.000000 |
50% | 96.567972 | 5.616401 | 5.598076 | 6.155411 | 5.667071 | 5.622082 |
75% | 173.452541 | 89.248955 | 84.904156 | 103.102595 | 98.899346 | 89.601183 |
max | 558.131136 | 531.675901 | 526.219164 | 533.076491 | 533.061174 | 535.251475 |
TMY[models].dropna().plot(kind='density')
<matplotlib.axes._subplots.AxesSubplot at 0x7f50b79a52d0>
ax_daily = daily[models].plot(title='Daily diffuse irradiation')
ylabel = ax_daily.set_ylabel('Irradiation [kWh]')
ax_monthly = monthly[models].plot(title='Monthly diffuse irradiation', kind='bar')
ylabel = ax_monthly.set_ylabel('Irradiation [kWh]')
yearly.plot(kind='barh')
<matplotlib.axes._subplots.AxesSubplot at 0x7f50b7712e10>
mean_yearly = yearly.mean()
yearly_mean_deviation = (yearly - mean_yearly) / yearly * 100.0
yearly_mean_deviation.plot(kind='bar')
<matplotlib.axes._subplots.AxesSubplot at 0x7f50b5a5dad0>