Warn when SPA is reloaded

pvlib/forecast_tools.py

@@ -0,0 +1,189 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Aug 23 15:03:04 2018
+
+@author: cwhanse
+"""
+
+import pandas as pd
+import statsmodels.api as sm
+from datetime import timedelta
+
+# Forecast functions
+
+def get_num_intervals(end, start, deltat):
+    """
+    Calculate number of intervals of length deltat from end working back to
+    start, to include start.
+
+    Parameters
+    ------------
+    end : datetime
+
+    start : datetime
+
+    deltat : timedelta
+
+    """
+    return int((end - start).total_seconds() / deltat.seconds) + 1
+
+
+def _get_data_for_ARMA_forecast(start, end, deltat, history,
+                                dataWindowLength):
+
+    """
+    Select data from history for fitting the forecast model.
+
+    Parameters
+    -----------
+
+    start : datetime
+        the first time in the forecast
+
+    end : datetime
+        the last time in be forecast
+
+    deltat : timedelta
+        the time interval for the forecast
+
+    history : pandas Series
+        historical values of the quantity to be forecast
+
+    Returns
+    ----------
+    fitdata : pandas Series
+        data from history aligned to be in phase with requested forecast
+
+    f_intervals : integer
+        number of time steps in the forecast
+    """
+
+    # find number of deltat intervals between start and last time in history
+    # assumes that start > max(history.index)
+    K = get_num_intervals(start, max(history.index), deltat)
+
+    # find number of deltat intervals covering dataWindowLength in history
+    N = get_num_intervals(start - K*deltat,
+                          start - K*deltat - dataWindowLength,
+                          deltat)
+
+    # start time of resampled history in phase with forecast start
+    resample_start = start - (K + N - 1) * deltat
+    resample_end = start - K * deltat
+
+    # calculate base time for resample, minutes out of phase with midnight
+    midnight = resample_start.replace(hour=0, minute=0, second=0)
+    first_after_midnight = resample_start - \
+                    int((resample_start - midnight) / deltat) * deltat
+    base = (first_after_midnight - midnight).total_seconds() / 60
+
+    # resample history
+    idata = history.resample(rule=pd.to_timedelta(deltat),
+                             closed='left',
+                             label='left',
+                             base=base).mean()
+
+    # extract window from resampled history
+    return idata.loc[(idata.index>=resample_start) &
+                        (idata.index<=resample_end)].copy()
+
+
+def forecast_ARMA(start, end, history, deltat,
+                  dataWindowLength=timedelta(hours=1),
+                  order=None,
+                  start_params=None):
+
+    """
+    Generate forecast from start to end with steps deltat
+    using an ARMA model fit to data in history.
+
+    Parameters
+    -----------
+
+    start : datetime
+        the first time to be forecast
+
+    end : datetime
+        the last time to be forecast
+
+    deltat : timedelta
+        the time step for the forecast
+
+    history : pandas Series
+        historical values from which the forecast is made.
+
+    dataWindowLenth : timedelta, default one hour
+        time interval for data in history to be considered for model fitting
+
+    order : tuple of three integers
+        autoregressive, difference, and moving average orders for an ARMA
+        forecast model
+
+    start_params : list of float
+        initial guess of model parameters, one value for each autoregressive
+        and moving average coefficient, followed by the value for the variance
+        term
+
+
+    """
+
+    # TODO: input validation
+
+    # create datetime index for forecast
+    fdr = pd.DatetimeIndex(start=start, end=end, freq=pd.to_timedelta(deltat))
+
+    # get time-averaged data from history over specified data window and that
+    # is in phase with forecast times
+    fitdata = _get_data_for_ARMA_forecast(start, end,
+                                          deltat, history,
+                                          dataWindowLength)
+
+    # TODO: model identification logic
+
+    # fit model of order (p, d, q)
+    # defaults to first order differencing to help with non-stationarity
+    if not order:
+        if deltat.total_seconds()>=15*60:
+            # use MA process for time intervals 15 min or longer
+            p = 0
+            d = 1
+            q = 1
+            start_params = [1, 0]
+        else:
+            # use AR process for time intervals < 15 min
+            p = 1
+            d = 1
+            q = 0
+            start_params = [0, 1]
+        order = (p, d, q)
+
+    model = sm.tsa.statespace.SARIMAX(fitdata,
+                                      trend='n',
+                                      order=order)
+
+    # if not provided, generate guess of model parameters, helps overcome 
+    # non-stationarity
+    if not start_params:
+        # generate a list with one entry '0' for each AR or MA parameter 
+        # plus an entry '1' for the variance parameter
+        start_params = []
+        for i in range(0, order[0]+order[2]):
+            start_params.append(0)
+        start_params.append(1)
+
+    # generate the ARMA model object
+    results = model.fit(start_params=start_params)
+
+    # total intervals to forecast from end of data to end of forecast
+    idr = pd.DatetimeIndex(start=max(fitdata.index),
+                           end=end,
+                           freq=pd.to_timedelta(deltat))
+
+    f_intervals = len(idr-1) # first time in idr is last data point
+
+    # forecast
+    f = results.forecast(f_intervals)
+    # return the requested forecast times
+
+    return f[fdr]
+

pvlib/irradiance.py

@@ -15,6 +15,7 @@
 
 from pvlib import atmosphere, solarposition, tools
 from pvlib._deprecation import deprecated
+from tools import nanmaximum
 
 # see References section of grounddiffuse function
 SURFACE_ALBEDOS = {'urban': 0.18,
@@ -1147,24 +1148,24 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     F2c = np.vstack((F2c, nans))
 
     F1 = (F1c[ebin, 0] + F1c[ebin, 1] * delta + F1c[ebin, 2] * z)
-    F1 = np.maximum(F1, 0)
+    F1 = nanmaximum(F1, 0)
 
     F2 = (F2c[ebin, 0] + F2c[ebin, 1] * delta + F2c[ebin, 2] * z)
-    F2 = np.maximum(F2, 0)
+    F2 = nanmaximum(F2, 0)
 
     A = aoi_projection(surface_tilt, surface_azimuth,
                        solar_zenith, solar_azimuth)
-    A = np.maximum(A, 0)
+    A = nanmaximum(A, 0)
 
     B = tools.cosd(solar_zenith)
-    B = np.maximum(B, tools.cosd(85))
+    B = nanmaximum(B, tools.cosd(85))
 
     # Calculate Diffuse POA from sky dome
     term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))
     term2 = F1 * A / B
     term3 = F2 * tools.sind(surface_tilt)
 
-    sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)
+    sky_diffuse = nanmaximum(dhi * (term1 + term2 + term3), 0)
 
     # we've preserved the input type until now, so don't ruin it!
     if isinstance(sky_diffuse, pd.Series):

pvlib/solarposition.py

@@ -21,6 +21,7 @@
 
 import numpy as np
 import pandas as pd
+import warnings
 
 from pvlib import atmosphere
 from pvlib.tools import datetime_to_djd, djd_to_datetime
@@ -242,12 +243,14 @@ def _spa_python_import(how):
         # reload the module without compiling
         # the PVLIB_USE_NUMBA env variable is used to tell the module
         # to not compile with numba
+        warnings.warn('Reloading spa to use numpy')
         os.environ['PVLIB_USE_NUMBA'] = '0'
         spa = reload(spa)
         del os.environ['PVLIB_USE_NUMBA']
     elif how == 'numba' and not using_numba:
         # The spa module was not compiled to numba code, so set
         # PVLIB_USE_NUMBA so it does compile to numba on reload.
+        warnings.warn('Reloading spa to use numba')
         os.environ['PVLIB_USE_NUMBA'] = '1'
         spa = reload(spa)
         del os.environ['PVLIB_USE_NUMBA']