Merge pull request #43 from NREL/combined_update

Combine update

.gitignore

@@ -10,3 +10,4 @@
 .DS_Store
 
 docs/.ipynb_checkpoints/degradation_example-checkpoint.ipynb
+*.ipynb_checkpoints*

MANIFEST.in

@@ -1,2 +1,3 @@
 include versioneer.py
 include rdtools/_version.py
+include rdtools/data/*

rdtools/__init__.py

@@ -1,9 +1,15 @@
 from normalization import normalize_with_sapm
 from normalization import normalize_with_pvwatts
+from normalization import irradiance_rescale
 from degradation import degradation_ols
 from degradation import degradation_classical_decomposition
 from degradation import degradation_year_on_year
 from aggregation import aggregation_insol
+from clearsky_temperature import get_clearsky_tamb
+from filtering import csi_filter
+from filtering import poa_filter
+from filtering import tcell_filter
+from filtering import clip_filter
 
 from ._version import get_versions
 __version__ = get_versions()['version']

rdtools/clearsky_temperature.py

@@ -0,0 +1,110 @@
+import h5py
+from numpy import arange
+from datetime import timedelta
+import math
+import pandas as pd
+import pkg_resources
+import numpy as np
+
+
+def get_clearsky_tamb(times, latitude, longitude, window_size=40, gauss_std=20):
+    '''
+    :param times:       DateTimeIndex in local time
+    :param latitude:    float degrees
+    :param longitude:   float degrees
+    :return:            pandas Series of cell sky ambient temperature
+    '''
+
+    filepath = pkg_resources.resource_filename('rdtools', 'data/temperature.hdf5')
+
+    buffer = timedelta(days=window_size)
+    freq_actual = pd.infer_freq(times)
+    dt_daily = pd.date_range(times[0] - buffer, times[-1] + buffer, freq='D')
+
+    f = h5py.File(filepath, "r")
+
+    a = f['temperature']['day']
+    b = f['temperature']['night']
+
+    lons = len(a[:, 0, 0])
+    lats = len(a[0, :, 0])
+
+    lon_temp = longitude - 180
+    if lon_temp < 0:
+        lon_temp += 360
+    lon_index = round(float(lons) * float(lon_temp) / 360.0)
+    lat_index = round(float(lats) * (90.0 - float(latitude)) / 180.0)
+
+    df = pd.DataFrame(index=dt_daily)
+    df['month'] = df.index.month
+
+    ave_day = []
+    ave_night = []
+
+    radius = 0
+    for k in range(12):
+
+        day = _get_pixel_value(a, lon_index, lat_index, k, radius)
+        night = _get_pixel_value(b, lon_index, lat_index, k, radius)
+
+        if day == float("NaN"):
+            day = a[:, lat_index, k]
+        if night == float("NaN"):
+            night = a[:, lat_index, k]
+
+        ave_day.append(day)
+        ave_night.append(night)
+
+    for i in range(12):
+        df.loc[df['month'] == i + 1, 'day'] = ave_day[i]
+        df.loc[df['month'] == i + 1, 'night'] = ave_night[i]
+
+    df = df.rolling(window=window_size, win_type='gaussian', min_periods=1, center=True).mean(std=gauss_std)
+
+    df = df.resample(freq_actual).interpolate(method='linear')
+    df['month'] = df.index.month
+
+    df = df.reindex(times, method='nearest')
+
+    utc_offsets = [y.utcoffset().total_seconds() / 3600.0 for y in df.index]
+
+    def solar_noon_offset(utc_offset):
+        return longitude / 180.0 * 12.0 - utc_offset
+
+    df['solar_noon_offset'] = solar_noon_offset(np.array(utc_offsets))
+
+    df['hour_of_day'] = df.index.hour + df.index.minute / 60.0
+    df['Clear Sky Temperature (C)'] = _get_temperature(df['hour_of_day'].values, df['night'].values,\
+                                                       df['day'].values, df['solar_noon_offset'].values)
+    return df['Clear Sky Temperature (C)']
+
+
+def _get_pixel_value(data, i, j, k, radius):
+    list = []
+    for x in arange(i - radius, i + radius + 1):
+        if x < 0 or x >= len(data[:, 0, 0]):
+            continue
+        for y in arange(j - radius, j + radius + 1):
+            if y < 0 or y >= len(data[0, :, 0]):
+                continue
+
+            value = data[x, y, k]
+            if value == float("NaN"):
+                continue
+
+            list.append(value)
+
+    if len(list) == 0:
+        return float("NaN")
+
+    return pd.Series(list).median()
+
+
+def _get_temperature(hour_of_day, night_temp, day_temp, solar_noon_offset):
+    hour_offset = 8.0 + solar_noon_offset
+    temp_scaler = 0.7
+    t_diff = day_temp - night_temp
+    t_ave = (day_temp + night_temp) / 2.0
+    v = np.cos((hour_of_day + hour_offset) / 24.0 * 2.0 * np.pi)
+    t = t_diff * 0.5 * v * temp_scaler + t_ave
+    return t

rdtools/degradation.py

@@ -168,7 +168,7 @@ def degradation_classical_decomposition(normalized_energy):
     return (Rd_pct, Rd_CI, calc_info)
 
 
-def degradation_year_on_year(normalized_energy, recenter=True):
+def degradation_year_on_year(normalized_energy, recenter=True, exceedance_prob=95):
     '''
     Description
     -----------
@@ -180,6 +180,7 @@ def degradation_year_on_year(normalized_energy, recenter=True):
         corrected performance ratio timeseries index in monthly format
     recenter:  bool, default value True
         specify whether data is centered to normalized yield of 1 based on first year
+    exceedance_prob (float): the probability level to use for exceedance value calculation
 
     Returns
     -------
@@ -191,6 +192,8 @@ def degradation_year_on_year(normalized_energy, recenter=True):
         calc_info:  dict
             ('YoY_values') pandas series of right-labeled year on year slopes
             ('renormalizing_factor') float of value used to recenter data
+            ('exceedance_level') the degradation rate that was ouperformed with
+            probability of exceedance_prob
     '''
 
     # Ensure the data is in order
@@ -199,7 +202,7 @@ def degradation_year_on_year(normalized_energy, recenter=True):
     normalized_energy.index.name = 'dt'
 
     # Detect sub-daily data:
-    if min(np.diff(normalized_energy.index.values, n=1)) < np.timedelta64(24, 'h'):
+    if min(np.diff(normalized_energy.index.values, n=1)) < np.timedelta64(23, 'h'):
         raise ValueError('normalized_energy must not be more frequent than daily')
 
     # Detect less than 2 years of data
@@ -233,19 +236,26 @@ def degradation_year_on_year(normalized_energy, recenter=True):
 
     yoy_result = df.yoy.dropna()
 
+    calc_info = {
+        'YoY_values': yoy_result,
+        'renormalizing_factor': renorm
+    }
+
+    if not len(yoy_result):
+        raise ValueError('no year-over-year aggregated data pairs found')
+
     Rd_pct = yoy_result.median()
 
-    # bootstrap to determine 68% CI
+    # bootstrap to determine 68% CI and exceedance probability
     n1 = len(yoy_result)
     reps = 10000
     xb1 = np.random.choice(yoy_result, (n1, reps), replace=True)
     mb1 = np.median(xb1, axis=0)
     Rd_CI = np.percentile(mb1, [15.9, 84.1])
 
-    calc_info = {
-        'YoY_values': yoy_result,
-        'renormalizing_factor': renorm
-    }
+    P_level = np.percentile(mb1, 100 - exceedance_prob)
+
+    calc_info['exceedance_level'] = P_level
 
     return (Rd_pct, Rd_CI, calc_info)
 

rdtools/filtering.py

@@ -0,0 +1,58 @@
+import pandas as pd
+
+
+def poa_filter(poa, low_irradiance_cutoff=200, high_irradiance_cutoff=1200):
+    # simple filter based on irradiance sensors
+    return (poa > low_irradiance_cutoff) & (poa < high_irradiance_cutoff)
+
+
+def tcell_filter(tcell, low_tcell_cutoff=-50, high_tcell_cutoff=110):
+    # simple filter based on temperature sensors
+    return (tcell > low_tcell_cutoff) & (tcell < high_tcell_cutoff)
+
+
+def clip_filter(power, quant=0.98, low_power_cutoff=0.01):
+    '''
+    Filter data points likely to be affected by clipping
+    with power greater than or equal to 99% of the 'quant'
+    quantile and less than 'low_power_cutoff'
+
+    Parameters
+    ----------
+    power: Pandas series (numeric)
+        AC power
+    quant: float
+        threshold for quantile
+    low_power_cutoff
+
+    Returns
+    -------
+    Pandas Series (boolean)
+        mask to exclude points equal to and
+        above 99% of the percentile threshold
+    '''
+    v = power.quantile(quant)
+    return (power < v * 0.99) & (power > low_power_cutoff)
+
+
+def csi_filter(measured_poa, clearsky_poa, threshold=0.15):
+    '''
+    Filtering based on clear sky index (csi)
+
+    Parameters
+    ----------
+    measured_poa: Pandas series (numeric)
+        Plane of array irradiance based on measurments
+    clearsky_poa: Pandas series (numeric)
+        Plane of array irradiance based on a clear sky model
+    threshold: float
+        threshold for filter
+
+    Returns
+    -------
+    Pandas Series (boolean)
+        mask to exclude points below the threshold
+    '''
+
+    csi = measured_poa / clearsky_poa
+    return (csi >= 1.0 - threshold) & (csi <= 1.0 + threshold)

rdtools/normalization.py

@@ -7,6 +7,7 @@
 import pandas as pd
 import pvlib
 import numpy as np
+from scipy.optimize import minimize
 
 
 def pvwatts_dc_power(poa_global, P_ref, T_cell=None, G_ref=1000, T_ref=25, gamma_pdc=None):
@@ -191,7 +192,6 @@ def sapm_dc_power(pvlib_pvsystem, met_data):
         .pvwatts_dc(g_poa_effective=effective_poa,
                     temp_cell=temp_cell['temp_cell'])
 
-
     return dc_power, effective_poa
 
 
@@ -271,3 +271,31 @@ def delta_index(series):
 
     deltas = (series.index - series.index.shift(-1)).astype('int64') / (10.0**9 * 3600.0)
     return deltas, np.mean(deltas)
+
+
+def irradiance_rescale(irrad, modeled_irrad):
+    '''
+    Attempts to rescale modeled irradiance to match measured irradiance on clear days
+    Parameters
+    ----------
+    irrad: Pandas Series (numeric)
+        measured irradiance time series
+    modeled_irrad: Pandas Series (numeric)
+        modeled irradiance time series
+    Returns
+    -------
+    Pandas Series (numeric): resacaled modeled irradaince time series
+    '''
+    def _rmse(fact):
+        rescaled_modeled_irrad = fact * modeled_irrad
+        csi = irrad / rescaled_modeled_irrad
+        filt = (csi >= 0.8) & (csi <= 1.2)
+        rmse = np.sqrt(((rescaled_modeled_irrad[filt] - irrad[filt]) ** 2.0).mean())
+        return rmse
+
+    guess = np.percentile(irrad.dropna(), 90) / np.percentile(modeled_irrad.dropna(), 90)
+    min_result = minimize(_rmse, guess, method='Nelder-Mead')
+    factor = min_result['x'][0]
+
+    out_irrad = factor * modeled_irrad
+    return out_irrad

requirements.txt

@@ -1,6 +1,6 @@
 nose==1.3.7
 numpy==1.11.2
-pandas==0.19.1
+pandas==0.19.2
 patsy==0.4.1
 pvlib==0.4.1
 python-dateutil==2.6.0

tests/clearsky_temperature_test.py

@@ -0,0 +1,41 @@
+import unittest
+import pandas as pd
+from datetime import datetime
+
+from rdtools.clearsky_temperature import get_clearsky_tamb
+
+
+class ClearSkyTemperatureTestCase(unittest.TestCase):
+    '''Unit tests for clearsky_temperature module'''
+
+    def setUp(self):
+
+        dt = pd.date_range(datetime(2015,1,1), datetime(2015,2,1), freq='15min', tz = 'Asia/Shanghai')
+
+        self.china_west = get_clearsky_tamb(dt, 37.951721, 80.609843)
+        self.china_east = get_clearsky_tamb(dt, 36.693692, 117.699686)
+
+        #self.china_west.to_csv("west.csv")
+        #self.china_east.to_csv("east.csv")
+
+
+    # Test for shifting temperature peak with longitude for same timezone
+    def test_hour_offset(self):
+
+        df = pd.DataFrame(index = self.china_west.index)
+        df['west'] = self.china_west
+        df['east'] = self.china_east
+        df['hour'] = df.index.hour
+
+        west_hottest_hour = df.sort_values(by='west', ascending=False)['hour'].iloc[0]
+        east_hottest_hour = df.sort_values(by='east', ascending=False)['hour'].iloc[0]
+
+        #print west_hottest_hour , east_hottest_hour
+
+        self.assertTrue(west_hottest_hour > 12)
+        self.assertTrue(east_hottest_hour > 12)
+        self.assertTrue(west_hottest_hour > east_hottest_hour)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/run_tests

@@ -3,3 +3,4 @@
 python -m tests.degradation_test
 python -m tests.normalization_pvwatts_test
 python -m tests.normalization_sapm_test
+python -m tests.clearsky_temperature_test