diff --git a/pvlib/forecast_tools.py b/pvlib/forecast_tools.py new file mode 100644 index 0000000000..b431f931b1 --- /dev/null +++ b/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] + diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py index eb63371b81..e012059940 100644 --- a/pvlib/irradiance.py +++ b/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): diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index ad67a9efa9..e4d95fcf3e 100644 --- a/pvlib/solarposition.py +++ b/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']