Repository: pvlib/pvanalytics

Issue:
Detect sunny days
The PVFleets QA project does this by fitting curves (quadratic for fixed systems, quartic for tracking systems) on each day. 

Discussion:

Patch:
diff --git a/.travis.yml b/.travis.yml
index dd6c7fb..ecc1892 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -7,8 +7,8 @@ python:
   - "3.8"
 
 install:
+  - pip install pytest==5.4.3 pytest-cov coveralls
   - pip install -r requirements.txt
-  - pip install pytest pytest-cov coveralls
 
 script:
   - pytest --cov=pvanalytics --cov-report term-missing
diff --git a/docs/api.rst b/docs/api.rst
index e6b92fa..a1b6a72 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -51,23 +51,7 @@ Identify gaps in the data.
    :toctree: generated/
 
    quality.gaps.interpolation_diff
-
-Data sometimes contains sequences of values that are "stale" or
-"stuck." These are contiguous spans of data where the value does not
-change within the precision given. The functions below
-can be used to detect stale values.
-
-.. note::
-
-   If the data has been altered in some way (i.e. temperature that has
-   been rounded to an integer value) before being passed to these
-   functions you may see unexpectedly large amounts of stale data.
-
-.. autosummary::
-   :toctree: generated/
-
    quality.gaps.stale_values_diff
-   quality.gaps.stale_values_round
 
 The following functions identify days with incomplete data.
 
@@ -134,6 +118,12 @@ Quality checks for weather data.
    quality.weather.temperature_limits
    quality.weather.wind_limits
 
+.. rubric:: References
+
+.. [1]  C. N. Long and Y. Shi, An Automated Quality Assessment and Control
+        Algorithm for Surface Radiation Measurements, The Open Atmospheric
+        Science Journal 2, pp. 23-37, 2008.
+
 Features
 ========
 
@@ -148,7 +138,6 @@ Functions for identifying inverter clipping
    :toctree: generated/
 
    features.clipping.levels
-   features.clipping.threshold
 
 Clearsky
 --------
@@ -158,8 +147,25 @@ Clearsky
 
    features.clearsky.reno
 
-.. rubric:: References
+Orientation
+-----------
 
-.. [1]  C. N. Long and Y. Shi, An Automated Quality Assessment and Control
-        Algorithm for Surface Radiation Measurements, The Open Atmospheric
-        Science Journal 2, pp. 23-37, 2008.
+System orientation refers to mounting type (fixed or tracker) and the
+azimuth and tilt of the mounting. A system's orientation can be
+determined by examining power or POA irradiance on days that are
+relatively sunny.
+
+This module provides functions that operate on power or POA irradiance
+to identify system orientation on a daily basis. These functions can
+tell you whether a day's profile matches that of a fixed system or
+system with a single-axis tracker.
+
+Care should be taken when interpreting function output since
+other factors such as malfunctioning trackers can interfere with
+identification.
+
+.. autosummary::
+   :toctree: generated/
+
+   features.orientation.fixed_nrel
+   features.orientation.tracking_nrel
diff --git a/pvanalytics/features/clipping.py b/pvanalytics/features/clipping.py
index d589781..7bb14ac 100644
--- a/pvanalytics/features/clipping.py
+++ b/pvanalytics/features/clipping.py
@@ -88,138 +88,3 @@ def levels(ac_power, window=4, fraction_in_window=0.75,
         flags = flags | _label_clipping(temp, window=window,
                                         frac=fraction_in_window)
     return flags
-
-
-def _daytime_powercurve(ac_power, power_quantile):
-    # return the `power_quantile` quantile of power at each minute of
-    # the day.
-    positive_power = ac_power[ac_power >= 0]
-    return positive_power.groupby(
-        positive_power.index.hour * 60 + positive_power.index.minute
-    ).quantile(power_quantile)
-
-
-def _clipped(power, slope, power_min, slope_max):
-    # return a mask that is true where `power` is greater than
-    # `power_min` and the absolute value of `slope` is less
-    # than or equal to `slope_max`
-    return (np.abs(slope) <= slope_max) & (power > power_min)
-
-
-def _clipping_power(ac_power, slope_max, power_min,
-                    power_quantile, freq=None):
-    # Calculate a power threshold, above which the power is being
-    # clipped.
-    #
-    # - The daytime power curve is calculated using
-    #   _daytime_powercurve(). This function groups `ac_power` by
-    #   minute of the day and returns the `power_quantile`-quantile of
-    #   power for each minute.
-    #
-    # - Each timestamp in the daytime power curve that satisfies the
-    #   clipping criteria[*] is flagged.
-    #
-    # - The clipping threshold is calculated as the mean power during the
-    #   longest flagged period in the daytime power curve.
-    #
-    # [*] clipping criteria: a timestamp satisfies the clipping
-    # criteria if the absolute value of the slope of the daytime power curve
-    # is less than `slope_max` and the value of the daytime
-    # power curve is greater than `power_min` times the median of the
-    # daytime power curve.
-    #
-    # Based on the PVFleets QA Analysis project
-    if not freq:
-        freq = pd.Timedelta(pd.infer_freq(ac_power.index)).seconds / 60
-    elif isinstance(freq, str):
-        freq = pd.Timedelta(freq).seconds / 60
-
-    # Use the slope of the 99.5% quantile of daytime power at
-    # each minute to identify clipping.
-    powercurve = _daytime_powercurve(
-        ac_power, power_quantile
-    )
-    normalized_power = powercurve / powercurve.max()
-    power_slope = (normalized_power.diff()
-                   / normalized_power.index.to_series().diff()) * freq
-
-    clipped_times = _clipped(powercurve, power_slope,
-                             powercurve.median() * power_min,
-                             slope_max)
-    clipping_cumsum = (~clipped_times).cumsum()
-    # get the value of the cumulative sum of the longest True span
-    longest_clipped = clipping_cumsum.value_counts().idxmax()
-    # select the longest span that satisfies the clipping criteria
-    longest = powercurve[clipping_cumsum == longest_clipped]
-
-    if longest.index.max() - longest.index.min() >= 60:
-        # if the period of clipping is at least 60 minutes then we
-        # have enough evidence to determine the clipping threshold.
-        return longest.mean()
-    return None
-
-
-def threshold(ac_power, slope_max=0.0035, power_min=0.75,
-              power_quantile=0.995, freq=None):
-    """Detect clipping based on a maximum power threshold.
-
-    This is a two-step process. First a clipping threshold is
-    identified, then any values in `ac_power` greater than or equal to
-    that threshold are flagged.
-
-    The clipping threshold is determined by computing a 'daily power
-    curve' which is the `power_quantile` quantile of all values in
-    `ac_power` at each minute of the day. This gives a rough estimate
-    of the maximum power produced at each minute of the day.
-
-    The daily power curve is normalized by its maximum and the minutes
-    of the day are identified where the normalized curve's slope is
-    less than `slope_max`. If there is a continuous period of time
-    spanning at least one hour where the slope is less than
-    `slope_max` and the value of the normalized daily power curve is
-    greater than `power_min` times the median of the normalized daily
-    power curve then the data has clipping in it. If no sufficiently
-    long period with both a low slope and high power exists then there
-    is no clipping in the data.  The average of the daily power curve
-    (not normalized) during the longest period that satisfies the
-    criteria above is the clipping threshold.
-
-    Parameters
-    ----------
-    ac_power : Series
-        DatetimeIndexed series of AC power data.
-    slope_max : float, default 0.0035
-        Maximum absolute value of slope of AC power quantile for
-        clipping to be indicated. The default value has been derived
-        empirically to prevent false positives for tracking PV
-        systems.
-    power_min: float, default 0.75
-        The power during periods with slope less than `slope_max` must
-        be greater than `power_min` times the median normalized
-        daytime power.
-    power_quantile : float, default 0.995
-        Quantile used to calculate the daily power curve.
-    freq : string, default None
-        A pandas string offset giving the frequency of data in
-        `ac_power`. If None then the frequency is inferred from the
-        series index.
-
-    Returns
-    -------
-    Series
-        True when `ac_power` is greater than or equal to the clipping
-        threshold.
-
-    Notes
-    -----
-    This function is based on the pvfleets_qa_analysis project.
-
-    """
-    threshold = _clipping_power(
-        ac_power,
-        slope_max,
-        power_min,
-        power_quantile,
-        freq=freq
-    )
-    return ac_power >= threshold
diff --git a/pvanalytics/features/orientation.py b/pvanalytics/features/orientation.py
new file mode 100644
index 0000000..f209534
--- /dev/null
+++ b/pvanalytics/features/orientation.py
@@ -0,0 +1,191 @@
+"""Functions for identifying system orientation."""
+import pandas as pd
+from pvanalytics.util import _fit, _group
+
+
+def _conditional_fit(day, fitfunc, freq, default=0.0, min_hours=0.0,
+                     peak_min=None):
+    # If there are at least `min_hours` of data in `day` and the
+    # maximum for the day is greater than `peak_min` then `fitfunc` is
+    # applied to fit a curve to the data. `fitfunc` must be a function
+    # that takes a Series and returns the :math:`r^2` for a curve fit.
+    # If the two conditions are not met then `default` is returned and
+    # no curve fitting is performed.
+    high_enough = True
+    if peak_min is not None:
+        high_enough = day.max() > peak_min
+    if (_hours(day, freq) > min_hours) and high_enough:
+        return fitfunc(day)
+    return default
+
+
+def _freqstr_to_hours(freq):
+    # Convert pandas freqstr to hours (as a float)
+    if freq.isalpha():
+        freq = '1' + freq
+    return pd.to_timedelta(freq).seconds / 60
+
+
+def _hours(data, freq):
+    # Return the number of hours in `data` with timestamp
+    # spacing given by `freq`.
+    return data.count() * _freqstr_to_hours(freq)
+
+
+def tracking_nrel(power_or_irradiance, daytime, r2_min=0.915,
+                  r2_fixed_max=0.96, min_hours=5, peak_min=None,
+                  quadratic_mask=None):
+    """Flag days that match the profile of a single-axis tracking PV system
+    on a sunny day.
+
+    This algorithm relies on the observation that the power profile of
+    a single-axis tracking PV system tends to resemble a quartic
+    polynomial on a sunny day, I.e., two peaks are observed, one
+    before and one after the sun crosses the tracker azimuth. By
+    contrast, the power profile for a fixed tilt PV system often
+    resembles a quadratic polynomial on a sunny day, with a single
+    peak when the sun is near the system azimuth.
+
+    The algorithm fits both a quartic and a quadratic polynomial to
+    each day's data.  A day is marked True if the quartic fit has a
+    sufficiently high :math:`r^2` and the quadratic fit has a
+    sufficiently low :math:`r^2`.  Specifically, a day is marked True
+    when three conditions are met:
+
+    1. a restricted quartic [#]_ must fit the data with :math:`r^2`
+       greater than `r2_min`
+    2. the :math:`r^2` for the restricted quartic fit must be greater
+       than the :math:`r^2` for a quadratic fit
+    3. the :math:`r^2` for a quadratic fit must be less than
+       `r2_fixed_max`
+
+    Values on days where any one of these conditions is not met are
+    marked False.
+
+    .. [#] The specific quartic used for this fit is centered within
+       70 minutes of 12:00, the y-value at the center must be within
+       15% of the median for the day, and it must open downwards.
+
+    Parameters
+    ----------
+    power_or_irradiance : Series
+        Timezone localized series of power or irradiance measurements.
+    daytime : Series
+        Boolean series with True for times that are during the
+        day. For best results this mask should exclude early morning
+        and late afternoon as well as night. Data at these times may have
+        problems with shadows that interfere with curve fitting.
+    r2_min : float, default 0.915
+        Minimum :math:`r^2` of a quartic fit for a day to be marked True.
+    r2_fixed_max : float, default 0.96
+        If the :math:`r^2` of a quadratic fit exceeds
+        `r2_fixed_max`, then tracking/fixed cannot be distinguished
+        and the day is marked False.
+    min_hours : float, default 5.0
+        Minimum number of hours with data to attempt a fit on a day.
+    peak_min : float, default None
+        The maximum `power_or_irradiance` value for a day must be
+        greater than `peak_min` for a fit to be attempted. If the
+        maximum for a day is less than `peak_min` then the day will be
+        marked False.
+    quadratic_mask : Series, default None
+        If None then `daytime` is used. This Series is used to remove
+        morning and afternoon times from the data before applying a
+        quadratic fit. The mask should
+        typically exclude more data than `daytime` in order to
+        eliminate long tails in the morning or afternoon that can
+        appear if a tracker is stuck in a West or East orientation.
+
+    Returns
+    -------
+    Series
+        Boolean series with True for every value on a day that has a
+        tracking profile (see criteria above).
+
+    Notes
+    -----
+    This algorithm is based on the PVFleets QA Analysis
+    project. Copyright (c) 2020 Alliance for Sustainable Energy, LLC.
+
+    """
+    if quadratic_mask is None:
+        quadratic_mask = daytime
+    freq = pd.infer_freq(power_or_irradiance.index)
+    daily_data = _group.by_day(power_or_irradiance[daytime])
+    tracking_days = daily_data.apply(
+        _conditional_fit,
+        _fit.quartic_restricted,
+        freq=freq,
+        min_hours=min_hours,
+        peak_min=peak_min
+    )
+    fixed_days = _group.by_day(power_or_irradiance[quadratic_mask]).apply(
+        _conditional_fit,
+        _fit.quadratic,
+        freq=freq,
+        min_hours=min_hours,
+        peak_min=peak_min
+    )
+    return (
+        (tracking_days > r2_min)
+        & (tracking_days > fixed_days)
+        & (fixed_days < r2_fixed_max)
+    ).reindex(power_or_irradiance.index, method='pad', fill_value=False)
+
+
+def fixed_nrel(power_or_irradiance, daytime, r2_min=0.94,
+               min_hours=5, peak_min=None):
+    """Flag days that match the profile of a fixed PV system on a sunny day.
+
+    This algorithm relies on the observation that the power profile of a
+    fixed tilt PV system often resembles a quadratic polynomial on a
+    sunny day, with a single peak when the sun is near the system azimuth.
+
+    A day is marked True when the :math:`r^2` for a quadratic fit to the
+    power data is greater than `r2_min`.
+
+    Parameters
+    ----------
+    power_or_irradiance : Series
+        Timezone localized series of power or irradiance measurements.
+    daytime : Series
+        Boolean series with True for times that are during the
+        day. For best results this mask should exclude early morning
+        and evening as well as night. Data at these times may have
+        problems with shadows that interfere with curve fitting.
+    r2_min : float, default 0.94
+        Minimum :math:`r^2` of a quadratic fit for a day to be marked True.
+    min_hours : float, default 5.0
+        Minimum number of hours with data to attempt a fit on a day.
+    peak_min : float, default None
+        The maximum `power_or_irradiance` value for a day must be
+        greater than `peak_min` for a fit to be attempted. If the
+        maximum for a day is less than `peak_min` then the day will be
+        marked False.
+
+    Returns
+    -------
+    Series
+        True for values on days where `power_or_irradiance` matches
+        the expected parabolic profile for a fixed PV system on a sunny day.
+
+    Notes
+    -----
+    This algorithm is based on the PVFleets QA Analysis
+    project. Copyright (c) 2020 Alliance for Sustainable Energy, LLC.
+
+    """
+    freq = pd.infer_freq(power_or_irradiance.index)
+    daily_data = _group.by_day(
+        power_or_irradiance[daytime]
+    )
+    fixed_days = daily_data.apply(
+        _conditional_fit,
+        _fit.quadratic,
+        freq=freq,
+        min_hours=min_hours,
+        peak_min=peak_min
+    )
+    return (
+        fixed_days > r2_min
+    ).reindex(power_or_irradiance.index, method='pad', fill_value=False)
diff --git a/pvanalytics/quality/gaps.py b/pvanalytics/quality/gaps.py
index 171d8ca..ff872da 100644
--- a/pvanalytics/quality/gaps.py
+++ b/pvanalytics/quality/gaps.py
@@ -51,17 +51,7 @@ def _backfill_window(endpoints, window):
     return flags
 
 
-def _mark(flags, window, mark):
-    if mark not in ['all', 'tail', 'end']:
-        raise ValueError("mark must be one of'tail', 'end' or 'all'")
-    if mark == 'all':
-        return _backfill_window(flags, window)
-    if mark == 'tail':
-        return _backfill_window(flags, window - 1)
-    return flags
-
-
-def stale_values_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
+def stale_values_diff(x, window=6, rtol=1e-5, atol=1e-8, label_all=False):
     """Identify stale values in the data.
 
     For a window of length N, the last value (index N-1) is considered
@@ -82,18 +72,9 @@ def stale_values_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
         relative tolerance for detecting a change in data values
     atol : float, default 1e-8
         absolute tolerance for detecting a change in data values
-    mark : str, default 'tail'
-        How much of the window to mark ``True`` when a sequence of
-        stale values is detected. Can one be of 'tail', 'end', or
-        'all'.
-
-        - If 'tail' (the default) then every point in the window
-          *except* the first point is marked ``True``.
-        - If 'end' then the first `window - 1` values in a stale
-          sequence sequence are marked ``False`` and all subsequent
-          values in the sequence are marked ``True``.
-        - If 'all' then every point in the window *including* the
-          first point is marked ``True``.
+    label_all : bool, default False
+        Whether to label all values in the window. If False, then only
+        the right endpoint of the window is labeled.
 
     Returns
     -------
@@ -103,8 +84,7 @@ def stale_values_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
     Raises
     ------
     ValueError
-        If `window < 2` or `mark` is not one of 'tail', 'end', or
-        'all'.
+        If window < 2.
 
     Notes
     -----
@@ -123,62 +103,12 @@ def stale_values_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
         raw=True,
         kwargs={'rtol': rtol, 'atol': atol}
     ).fillna(False).astype(bool)
-    return _mark(flags, window, mark)
-
-
-def stale_values_round(x, window=6, decimals=3, mark='tail'):
-    """Identify stale values by rounding.
-
-    A value is considered stale if it is part of a sequence of length
-    `window` of values that are identical when rounded to `decimals`
-    decimal places.
-
-    Parameters
-    ----------
-    x : Series
-        Data to be processed.
-    window : int, default 6
-        Number of consecutive identical values for a data point to be
-        considered stale.
-    decimals : int, default 3
-        Number of decimal places to round to.
-    mark : str, default 'tail'
-        How much of the window to mark ``True`` when a sequence of
-        stale values is detected. Can be one of 'tail', 'end', or
-        'all'.
-
-        - If 'tail' (the default) then every point in the window
-          *except* the first point is marked ``True``.
-        - If 'end' then the first `window - 1` values in a stale
-          sequence sequence are marked ``False`` and all subsequent
-          values in the sequence are marked ``True``.
-        - If 'all' then every point in the window *including* the
-          first point is marked ``True``.
-
-    Returns
-    -------
-    Series
-        True for each value that is part of a stale sequence of data.
-
-    Raises
-    ------
-    ValueError
-        If `mark` is not one of 'tail', 'end', or 'all'.
-
-    Notes
-    -----
-        Based on code from the pvfleets_qa_analysis project. Copyright
-        (c) 2020 Alliance for Sustainable Energy, LLC.
-
-    """
-    rounded_diff = x.round(decimals=decimals).diff()
-    endpoints = rounded_diff.rolling(window=window-1).apply(
-        lambda xs: len(xs[xs == 0]) == window-1
-    ).fillna(False).astype(bool)
-    return _mark(endpoints, window, mark)
+    if label_all:
+        return _backfill_window(flags, window)
+    return flags
 
 
-def interpolation_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
+def interpolation_diff(x, window=6, rtol=1e-5, atol=1e-8, label_all=False):
     """Identify sequences which appear to be linear.
 
     Sequences are linear if the first difference appears to be
@@ -199,18 +129,9 @@ def interpolation_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
         tolerance relative to max(abs(x.diff()) for detecting a change
     atol : float, default 1e-8
         absolute tolerance for detecting a change in first difference
-    mark : str, default 'tail'
-        How much of the window to mark ``True`` when a sequence of
-        interpolated values is detected. Can be one of 'tail', 'end',
-        or 'all'.
-
-        - If 'tail' (the default) then every point in the window
-          *except* the first point is marked ``True``.
-        - If 'end' then the first `window - 1` values in an
-          interpolated sequence are marked ``False`` and all
-          subsequent values in the sequence are marked ``True``.
-        - If 'all' then every point in the window *including* the
-          first point is marked ``True``.
+    label_all : bool, default False
+        Whether to label all values in the window. If False, then only the
+        right endpoint of the window is labeled.
 
     Returns
     -------
@@ -220,8 +141,7 @@ def interpolation_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
     Raises
     ------
     ValueError
-        If `window < 3` or `mark` is not one of 'tail', 'end', or
-        'all'.
+        If window < 3.
 
     Notes
     -----
@@ -240,10 +160,11 @@ def interpolation_diff(x, window=6, rtol=1e-5, atol=1e-8, mark='tail'):
         x.diff(periods=1),
         window=window-1,
         rtol=rtol,
-        atol=atol,
-        mark='end'
+        atol=atol
     )
-    return _mark(flags, window, mark)
+    if label_all:
+        return _backfill_window(flags, window)
+    return flags
 
 
 def _freq_to_seconds(freq):
diff --git a/pvanalytics/util/__init__.py b/pvanalytics/util/__init__.py
new file mode 100644
index 0000000..47cf6a3
--- /dev/null
+++ b/pvanalytics/util/__init__.py
@@ -0,0 +1,2 @@
+from pvanalytics.util import _fit    # noqa: F401
+from pvanalytics.util import _group  # noqa: F401
diff --git a/pvanalytics/util/_fit.py b/pvanalytics/util/_fit.py
new file mode 100644
index 0000000..183173b
--- /dev/null
+++ b/pvanalytics/util/_fit.py
@@ -0,0 +1,89 @@
+"""Internal module for curve fitting functions."""
+import numpy as np
+import scipy.optimize
+import pandas as pd
+
+
+def _to_minute_of_day(index):
+    # Transform the index into minutes of the day. If `index` is a
+    # DatetimeIndex then it is converted to an Int64Index with values
+    # equal to the minute of the day since midnight. Any other type
+    # and a ValueError is raised.
+    if isinstance(index, pd.DatetimeIndex):
+        return index.hour * 60 + index.minute
+    raise ValueError("cannot convert index to minutes since midnight")
+
+
+def quadratic(data):
+    """Fit a quadratic to the data.
+
+    Parameters
+    ----------
+    data : Series
+        Series of power or irradiance measurements.
+
+    Reurns
+    ------
+    float
+        The :math:`R^2` value for the fit.
+
+    Notes
+    -----
+    Based on the PVFleets QA Analysis project. Copyright (c) 2020
+    Alliance for Sustainable Energy, LLC.
+
+    """
+    minute_of_day = _to_minute_of_day(data.index)
+    # Fit a quadratic to `data` returning R^2 for the fit.
+    coefficients = np.polyfit(minute_of_day, data, 2)
+    quadratic = np.poly1d(coefficients)
+    # Calculate the R^2 for the fit
+    _, _, correlation, _, _ = scipy.stats.linregress(
+        data, quadratic(minute_of_day)
+    )
+    return correlation**2
+
+
+def quartic_restricted(data, noon=720):
+    """Fit a restricted quartic to the data.
+
+    The quartic is restricted to match the expected shape for a
+    tracking pv system under clearsky conditions. The quartic must:
+
+    - open downwards
+    - be centered within 70 minutes of the expected solar noon (see
+      `noon` parameter)
+    - the y-value at the center must be within 15% of the median of
+      `data`
+
+    Parameters
+    ----------
+    data : Series
+        power or irradiance data.
+    noon : int, default 720
+       The minute for solar noon. Defaults to the clock-noon.
+
+    Returns
+    -------
+    rsquared : float
+        The :math:`R^2` value for the fit.
+
+    Notes
+    -----
+    Based on the PVFleets QA Analysis project. Copyright (c) 2020
+    Alliance for Sustainable Energy, LLC.
+
+    """
+    def _quartic(x, a, b, c, e):
+        return a * (x - e)**4 + b * (x - e)**2 + c
+    minute_of_day = _to_minute_of_day(data.index)
+    median = data.median()
+    params, _ = scipy.optimize.curve_fit(
+        _quartic,
+        minute_of_day, data,
+        bounds=((-1e-05, 0, median * 0.85, noon - 70),
+                (-1e-10, median * 3e-05, median * 1.15, noon + 70))
+    )
+    model = _quartic(minute_of_day, params[0], params[1], params[2], params[3])
+    residuals = data - model
+    return 1 - (np.sum(residuals**2) / np.sum((data - np.mean(data))**2))
diff --git a/pvanalytics/util/_group.py b/pvanalytics/util/_group.py
new file mode 100644
index 0000000..4f3ebfd
--- /dev/null
+++ b/pvanalytics/util/_group.py
@@ -0,0 +1,27 @@
+"""Functions for grouping data"""
+import pandas as pd
+
+
+def by_day(data):
+    """Group data by day, preserving timezone in each group.
+
+    This function can be used in place of ``pd.Series.resample('D')``
+    when it is necessary to preserve the timezone in the index of each
+    group (e.g. when the function being applied to the Resampler
+    requires a localized time series).
+
+    Parameters
+    ----------
+    data : Series
+        DatetimeIndexed series.
+
+    Returns
+    -------
+    GroupBy
+        Data grouped by day, with a DatetimeIndex with daily frequency
+        and the same timezone as `data`.
+
+    """
+    return data.groupby(
+        pd.to_datetime(data.index.date).tz_localize(data.index.tz)
+    )


Test patch:
diff --git a/pvanalytics/tests/features/conftest.py b/pvanalytics/tests/conftest.py
similarity index 100%
rename from pvanalytics/tests/features/conftest.py
rename to pvanalytics/tests/conftest.py
diff --git a/pvanalytics/tests/features/test_clipping.py b/pvanalytics/tests/features/test_clipping.py
index 347ff40..add65af 100644
--- a/pvanalytics/tests/features/test_clipping.py
+++ b/pvanalytics/tests/features/test_clipping.py
@@ -2,7 +2,6 @@
 import pytest
 from pandas.util.testing import assert_series_equal
 import numpy as np
-import pandas as pd
 from pvanalytics.features import clipping
 
 
@@ -69,166 +68,3 @@ def test_levels_two_periods(quadratic, quadratic_clipped):
     assert clipped[35:40].all()
     assert not clipped[0:10].any()
     assert not clipped[50:].any()
-
-
-def test_threshold_no_clipping(quadratic):
-    """In a data set with a single quadratic there is no clipping."""
-    quadratic.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    assert not clipping.threshold(quadratic).any()
-
-
-def test_threshold_no_clipping_with_night(quadratic):
-    """In a data set with a single quadratic surrounded by zeros there is
-    no clipping."""
-    quadratic.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    full_day = quadratic.reindex(
-        pd.date_range(
-            start='01/01/2020 00:00',
-            end='01/01/2020 23:50',
-            freq='10T')
-    )
-    full_day.fillna(0)
-    assert not clipping.threshold(quadratic).any()
-
-
-def test_threshold_clipping(quadratic_clipped):
-    """In a data set with a single clipped quadratic clipping is
-    indicated."""
-    quadratic_clipped.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    assert not clipping.threshold(quadratic_clipped).all()
-    assert clipping.threshold(quadratic_clipped).any()
-
-
-def test_threshold_clipping_with_night(quadratic_clipped):
-    """Clipping is identified in the daytime with periods of zero power
-    before and after simulating night time conditions."""
-    quadratic_clipped.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    full_day = quadratic_clipped.reindex(
-        pd.date_range(
-            start='01/01/2020 00:00',
-            end='01/01/2020 23:50',
-            freq='10T')
-    )
-    full_day.fillna(0)
-    assert not clipping.threshold(full_day).all()
-    assert clipping.threshold(full_day)[quadratic_clipped.index].any()
-
-
-def test_threshold_clipping_with_freq(quadratic_clipped):
-    """Passing the frequency gives same result as infered frequency."""
-    quadratic_clipped.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    assert_series_equal(
-        clipping.threshold(quadratic_clipped),
-        clipping.threshold(quadratic_clipped, freq='10T')
-    )
-
-
-def test_threshold_clipping_with_interruption(quadratic_clipped):
-    """Test threshold clipping with period of no clipping mid-day."""
-    quadratic_clipped.loc[28:31] = [750, 725, 700, 650]
-    quadratic_clipped.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    clipped = clipping.threshold(quadratic_clipped)
-
-    assert not clipped.iloc[0:10].any()
-    assert not clipped.iloc[28:31].any()
-    assert not clipped.iloc[50:].any()
-    assert clipped.iloc[17:27].all()
-    assert clipped.iloc[32:40].all()
-
-
-def test_threshold_clipping_four_days(quadratic, quadratic_clipped):
-    """Clipping is identified in the first of four days."""
-    quadratic.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    quadratic_clipped.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    full_day_clipped = quadratic_clipped.reindex(
-        pd.date_range(
-            start='01/01/2020 00:00',
-            end='01/01/2020 23:50',
-            freq='10T')
-    )
-    full_day = quadratic.reindex(
-        pd.date_range(
-            start='01/01/2020 00:00',
-            end='01/01/2020 23:50',
-            freq='10T')
-    )
-    full_day_clipped.fillna(0)
-    full_day.fillna(0)
-
-    # scale the rest of the days below the clipping threshold
-    full_day *= 0.75
-
-    power = full_day_clipped
-    power.append(full_day)
-    power.append(full_day)
-    power.append(full_day)
-
-    power.index = pd.date_range(
-        start='01/01/2020 00:00', freq='10T', periods=len(power)
-    )
-
-    clipped = clipping.threshold(power)
-
-    assert clipped['01/01/2020'].any()
-    assert not clipped['01/02/2020':].any()
-
-
-def test_threshold_no_clipping_four_days(quadratic):
-    """Four days with no clipping."""
-    quadratic.index = pd.date_range(
-        start='01/01/2020 07:30',
-        freq='10T',
-        periods=61
-    )
-    full_day = quadratic.reindex(
-        pd.date_range(
-            start='01/01/2020 00:00',
-            end='01/01/2020 23:50',
-            freq='10T')
-    )
-    full_day.fillna(0)
-
-    power = full_day
-    power.append(full_day * 1.3)
-    power.append(full_day * 1.2)
-    power.append(full_day * 1.1)
-
-    power.index = pd.date_range(
-        start='01/01/2020 00:00', freq='10T', periods=len(power)
-    )
-
-    clipped = clipping.threshold(power)
-
-    assert not clipped.any()
diff --git a/pvanalytics/tests/features/test_orientation.py b/pvanalytics/tests/features/test_orientation.py
new file mode 100644
index 0000000..90bbcf5
--- /dev/null
+++ b/pvanalytics/tests/features/test_orientation.py
@@ -0,0 +1,159 @@
+import pytest
+from pandas.util.testing import assert_series_equal
+import pandas as pd
+from pvlib import location, pvsystem, tracking, modelchain, irradiance
+from pvlib.temperature import TEMPERATURE_MODEL_PARAMETERS
+from pvanalytics.features import orientation
+
+
+@pytest.fixture(scope='module')
+def times():
+    """Three days with one hour timestamp spacing in Etc/GMT+7"""
+    return pd.date_range(
+        start='03/01/2020',
+        end='03/04/2020',
+        closed='left',
+        freq='H',
+        tz='Etc/GMT+7'
+    )
+
+
+@pytest.fixture(scope='module')
+def albuquerque():
+    """pvlib Location for Albuquerque, NM."""
+    return location.Location(
+        35.0844,
+        -106.6504,
+        name='Albuquerque',
+        altitude=1500,
+        tx='Etc/GMT+7'
+    )
+
+
+@pytest.fixture(scope='module')
+def system_parameters():
+    """System parameters for generating simulated power data."""
+    sandia_modules = pvsystem.retrieve_sam('SandiaMod')
+    sapm_inverters = pvsystem.retrieve_sam('cecinverter')
+    module = sandia_modules['Canadian_Solar_CS5P_220M___2009_']
+    inverter = sapm_inverters['ABB__MICRO_0_25_I_OUTD_US_208__208V_']
+    temperature_model_parameters = (
+        TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
+    )
+    return {
+        'module_parameters': module,
+        'inverter_parameters': inverter,
+        'temperature_model_parameters': temperature_model_parameters
+    }
+
+
+@pytest.fixture(scope='module')
+def clearsky(times, albuquerque):
+    """Clearsky at `times` in `albuquerque`."""
+    return albuquerque.get_clearsky(times, model='simplified_solis')
+
+
+@pytest.fixture(scope='module')
+def solarposition(times, albuquerque):
+    """Solar position at `times` in `albuquerque`."""
+    return albuquerque.get_solarposition(times)
+
+
+def test_clearsky_ghi_fixed(clearsky, solarposition):
+    """Identify every day as fixed, since clearsky GHI is sunny."""
+    assert orientation.fixed_nrel(
+        clearsky['ghi'],
+        solarposition['zenith'] < 87,
+        r2_min=0.94,
+    ).all()
+
+
+@pytest.mark.filterwarnings("ignore:invalid value encountered in",
+                            "ignore:divide by zero encountered in")
+def test_perturbed_ghi_fixed(clearsky, solarposition):
+    """If the clearsky for one day is perturbed then that day is not sunny."""
+    ghi = clearsky['ghi']
+    ghi.iloc[0:24] = 1
+    expected = pd.Series(True, ghi.index)
+    expected[0:24] = False
+    assert_series_equal(
+        expected,
+        orientation.fixed_nrel(
+            ghi,
+            solarposition['zenith'] < 87
+        ),
+        check_names=False
+    )
+
+
+def test_ghi_not_tracking(clearsky, solarposition):
+    """If we pass GHI measurements and tracking=True then no days are sunny."""
+    assert (~orientation.tracking_nrel(
+        clearsky['ghi'], solarposition['zenith'] < 87
+    )).all()
+
+
+@pytest.fixture
+def power_tracking(clearsky, albuquerque, system_parameters):
+    """Simulated power for a pvlib SingleAxisTracker PVSystem in Albuquerque"""
+    system = tracking.SingleAxisTracker(**system_parameters)
+    mc = modelchain.ModelChain(
+        system,
+        albuquerque,
+        orientation_strategy='south_at_latitude_tilt'
+    )
+    mc.run_model(clearsky)
+    return mc.ac
+
+
+def test_power_tracking(power_tracking, solarposition):
+    """simulated power from a single axis tracker is identified as sunny
+    with tracking=True"""
+    assert orientation.tracking_nrel(
+        power_tracking,
+        solarposition['zenith'] < 87
+    ).all()
+
+
+def test_power_tracking_perturbed(power_tracking, solarposition):
+    """A day with perturbed values is not marked as tracking."""
+    power_tracking.iloc[6:18] = 10
+    expected = pd.Series(True, index=power_tracking.index)
+    expected.iloc[0:24] = False
+    assert_series_equal(
+        expected,
+        orientation.tracking_nrel(
+            power_tracking,
+            solarposition['zenith'] < 87
+        )
+    )
+    assert_series_equal(
+        expected,
+        orientation.tracking_nrel(
+            power_tracking,
+            solarposition['zenith'] < 87,
+            peak_min=100
+        )
+    )
+
+
+def test_stuck_tracker_profile(solarposition, clearsky):
+    """Test POA irradiance at a awkward orientation (high tilt and
+    oriented West)."""
+    poa = irradiance.get_total_irradiance(
+        surface_tilt=45,
+        surface_azimuth=270,
+        **clearsky,
+        solar_zenith=solarposition['apparent_zenith'],
+        solar_azimuth=solarposition['azimuth']
+    )
+    assert not orientation.tracking_nrel(
+        poa['poa_global'],
+        solarposition['zenith'] < 87,
+    ).any()
+    # by restricting the data to the middle of the day (lower zenith
+    # angles) we should classify the system as fixed
+    assert orientation.fixed_nrel(
+        poa['poa_global'],
+        solarposition['zenith'] < 70
+    ).all()
diff --git a/pvanalytics/tests/quality/test_gaps.py b/pvanalytics/tests/quality/test_gaps.py
index 6d008c7..6c27558 100644
--- a/pvanalytics/tests/quality/test_gaps.py
+++ b/pvanalytics/tests/quality/test_gaps.py
@@ -52,14 +52,14 @@ def test_stale_values_diff(stale_data):
     for more information.
 
     """
-    res0 = gaps.stale_values_diff(stale_data, mark='end')
-    res1 = gaps.stale_values_diff(stale_data, window=3, mark='end')
-    res2 = gaps.stale_values_diff(stale_data, rtol=1e-8, window=2, mark='end')
-    res3 = gaps.stale_values_diff(stale_data, window=7, mark='end')
-    res4 = gaps.stale_values_diff(stale_data, window=8, mark='end')
-    res5 = gaps.stale_values_diff(stale_data, rtol=1e-8, window=4, mark='end')
-    res6 = gaps.stale_values_diff(stale_data[1:], window=3, mark='end')
-    res7 = gaps.stale_values_diff(stale_data[1:8], window=3, mark='end')
+    res0 = gaps.stale_values_diff(stale_data)
+    res1 = gaps.stale_values_diff(stale_data, window=3)
+    res2 = gaps.stale_values_diff(stale_data, rtol=1e-8, window=2)
+    res3 = gaps.stale_values_diff(stale_data, window=7)
+    res4 = gaps.stale_values_diff(stale_data, window=8)
+    res5 = gaps.stale_values_diff(stale_data, rtol=1e-8, window=4)
+    res6 = gaps.stale_values_diff(stale_data[1:], window=3)
+    res7 = gaps.stale_values_diff(stale_data[1:8], window=3)
     assert_series_equal(res0, pd.Series([False, False, False, False, False,
                                          False, True, True, False, False]))
     assert_series_equal(res1, pd.Series([False, False, False, True, True, True,
@@ -91,22 +91,16 @@ def test_stale_values_diff_handles_negatives(data_with_negatives):
     for more information.
 
     """
-    res = gaps.stale_values_diff(data_with_negatives, window=3, mark='end')
+    res = gaps.stale_values_diff(data_with_negatives, window=3)
     assert_series_equal(res, pd.Series([False, False, True, True, False, False,
                                         False]))
-    res = gaps.stale_values_diff(
-        data_with_negatives, window=3, atol=1e-3, mark='end'
-    )
+    res = gaps.stale_values_diff(data_with_negatives, window=3, atol=1e-3)
     assert_series_equal(res, pd.Series([False, False, True, True, True, True,
                                         True]))
-    res = gaps.stale_values_diff(
-        data_with_negatives, window=3, atol=1e-5, mark='end'
-    )
+    res = gaps.stale_values_diff(data_with_negatives, window=3, atol=1e-5)
     assert_series_equal(res, pd.Series([False, False, True, True, True, False,
                                         False]))
-    res = gaps.stale_values_diff(
-        data_with_negatives, window=3, atol=2e-5, mark='end'
-    )
+    res = gaps.stale_values_diff(data_with_negatives, window=3, atol=2e-5)
     assert_series_equal(res, pd.Series([False, False, True, True, True, True,
                                         True]))
 
@@ -127,28 +121,14 @@ def test_stale_values_diff_raises_error(stale_data):
         gaps.stale_values_diff(stale_data, window=1)
 
 
-def test_stale_values_diff_raises_error_for_bad_mark(stale_data):
-    """Passing mark not in ['all', 'end', 'tail'] raises a ValueError."""
-    with pytest.raises(ValueError):
-        gaps.stale_values_diff(stale_data, mark='head')
-
-
-def test_stale_values_diff_mark_all(stale_data):
-    """When mark='all' the full window is marked stale"""
+def test_stale_values_diff_label_all(stale_data):
+    """When label_all is True the full window is marked stale"""
     assert_series_equal(
         pd.Series([False, True, True, True, True,
                    True, True, True, False, False]),
-        gaps.stale_values_diff(stale_data, window=4, mark='all')
-    )
-
-
-def test_stale_values_diff_mark_tail(stale_data):
-    """When mark='tail' (the default), every point in the window except
-    the first is marked stale."""
-    assert_series_equal(
-        pd.Series([False, False, True, True, True,
-                   True, True, True, False, False]),
-        gaps.stale_values_diff(stale_data, window=4)
+        gaps.stale_values_diff(
+            stale_data, window=4, label_all=True
+        )
     )
 
 
@@ -170,10 +150,10 @@ def interpolated_data():
     return pd.Series(data=data)
 
 
-def test_interpolation_diff_mark_all(interpolated_data):
-    """When mark='all' the full window is marked interpoated"""
+def test_interpolation_diff_label_all(interpolated_data):
+    """When label_all is True the full window is marked interpoated"""
     assert_series_equal(
-        gaps.interpolation_diff(interpolated_data, window=3, mark='all'),
+        gaps.interpolation_diff(interpolated_data, window=3, label_all=True),
         pd.Series([False, False, False, False, False,
                    True, True, True, False, False,
                    False, True, True, True, True, True,
@@ -181,18 +161,6 @@ def test_interpolation_diff_mark_all(interpolated_data):
     )
 
 
-def test_interpolation_diff_mark_tail(interpolated_data):
-    """When mark='tail' (the default), all but the first point an the
-    window is marked interpolated."""
-    assert_series_equal(
-        gaps.interpolation_diff(interpolated_data, window=3),
-        pd.Series([False, False, False, False, False,
-                   False, True, True, False, False,
-                   False, False, True, True, True, True,
-                   False])
-    )
-
-
 def test_interpolation_diff(interpolated_data):
     """Interpolation is detected correclty.
 
@@ -205,31 +173,27 @@ def test_interpolation_diff(interpolated_data):
     for more information.
 
     """
-    res0 = gaps.interpolation_diff(interpolated_data, mark='end')
+    res0 = gaps.interpolation_diff(interpolated_data)
     assert_series_equal(res0, pd.Series([False, False, False, False, False,
                                          False, False, False, False, False,
                                          False, False, False, False, False,
                                          False, False]))
-    res1 = gaps.interpolation_diff(interpolated_data, window=3, mark='end')
+    res1 = gaps.interpolation_diff(interpolated_data, window=3)
     assert_series_equal(res1, pd.Series([False, False, False, False, False,
                                          False, False, True, False, False,
                                          False, False, False, True, True, True,
                                          False]))
-    res2 = gaps.interpolation_diff(
-        interpolated_data, window=3, rtol=1e-2, mark='end'
-    )
+    res2 = gaps.interpolation_diff(interpolated_data, window=3, rtol=1e-2)
     assert_series_equal(res2, pd.Series([False, False, True, True, True,
                                          False, False, True, False, False,
                                          False, False, False, True, True, True,
                                          False]))
-    res3 = gaps.interpolation_diff(interpolated_data, window=5, mark='end')
+    res3 = gaps.interpolation_diff(interpolated_data, window=5)
     assert_series_equal(res3, pd.Series([False, False, False, False, False,
                                          False, False, False, False, False,
                                          False, False, False, False, False,
                                          True, False]))
-    res4 = gaps.interpolation_diff(
-        interpolated_data, window=3, atol=1e-2, mark='end'
-    )
+    res4 = gaps.interpolation_diff(interpolated_data, window=3, atol=1e-2)
     assert_series_equal(res4, pd.Series([False, False, True, True, True,
                                          True, True, True, False, False,
                                          False, False, False, True, True, True,
@@ -248,14 +212,10 @@ def test_interpolation_diff_handles_negatives(data_with_negatives):
     for more information.
 
     """
-    res = gaps.interpolation_diff(
-        data_with_negatives, window=3, atol=1e-5, mark='end'
-    )
+    res = gaps.interpolation_diff(data_with_negatives, window=3, atol=1e-5)
     assert_series_equal(res, pd.Series([False, False, True, True, True, True,
                                         False]))
-    res = gaps.stale_values_diff(
-        data_with_negatives, window=3, atol=1e-4, mark='end'
-    )
+    res = gaps.stale_values_diff(data_with_negatives, window=3, atol=1e-4)
     assert_series_equal(res, pd.Series([False, False, True, True, True, True,
                                         True]))
 
@@ -273,13 +233,7 @@ def test_interpolation_diff_raises_error(interpolated_data):
 
     """
     with pytest.raises(ValueError):
-        gaps.interpolation_diff(interpolated_data, window=2, mark='end')
-
-
-def test_interpolation_diff_raises_error_for_bad_mark(interpolated_data):
-    """Passing mark not in ['all', 'end', 'tail'] raises a ValueError."""
-    with pytest.raises(ValueError):
-        gaps.interpolation_diff(interpolated_data, mark='bad')
+        gaps.interpolation_diff(interpolated_data, window=2)
 
 
 def test_start_stop_dates_all_true():
@@ -651,94 +605,3 @@ def test_complete():
         pd.Series(True, index=data.index),
         gaps.complete(data, minimum_completeness=0.2)
     )
-
-
-def test_stale_values_round_no_stale():
-    """No stale values in a monotonically increasing sequence."""
-    data = pd.Series(np.linspace(0, 10))
-    assert not gaps.stale_values_round(data, mark='all').any()
-
-
-def test_stale_values_round_all_same():
-    """If all data is identical, then all values are stale."""
-    data = pd.Series(1, index=range(0, 10))
-    assert gaps.stale_values_round(data, mark='all').all()
-
-
-def test_stale_values_round_noisy():
-    """If all values are the same +/- 0.0005"""
-    data = pd.Series(
-        [1.555, 1.5551, 1.5549, 1.555, 1.555, 1.5548, 1.5553]
-    )
-    assert gaps.stale_values_round(data, decimals=3, mark='all').all()
-
-
-def test_stale_values_round_span_in_middle():
-    """A span of stale values between not-stale data."""
-    data = pd.Series(
-        [1.0, 1.1, 1.2, 1.5, 1.5, 1.5, 1.5, 1.9, 2.0, 2.2]
-    )
-    assert_series_equal(
-        gaps.stale_values_round(data, window=4, mark='all'),
-        pd.Series([False, False, False,
-                   True, True, True, True,
-                   False, False, False], dtype='bool')
-    )
-
-
-def test_stale_values_larger_window():
-    """Increasing the window size excludes short spans of repeated
-    values."""
-    data = pd.Series(
-        [1, 2, 2, 2, 2, 3, 4, 4, 4, 4, 4, 6]
-    )
-    assert_series_equal(
-        gaps.stale_values_round(data, window=4, mark='all'),
-        (data == 2) | (data == 4)
-    )
-    assert_series_equal(
-        gaps.stale_values_round(data, window=5, mark='all'),
-        (data == 4)
-    )
-
-
-def test_stale_values_round_bad_mark():
-    """passing an invalid value for `mark` raises a ValueError."""
-    data = pd.Series(1, index=range(1, 10))
-    with pytest.raises(ValueError):
-        gaps.stale_values_round(data, mark='bad')
-
-
-def test_stale_values_round_mark():
-    """Test that different values for `mark` have the correct semantics."""
-    data = pd.Series(1, index=range(0, 10))
-    expected = pd.Series(True, index=range(0, 10))
-    assert_series_equal(
-        expected,
-        gaps.stale_values_round(data, mark='all')
-    )
-    expected.iloc[0] = False
-    assert_series_equal(
-        expected,
-        gaps.stale_values_round(data)
-    )
-    assert_series_equal(
-        expected,
-        gaps.stale_values_round(data, mark='tail')
-    )
-    expected.iloc[1] = False
-    assert_series_equal(
-        expected,
-        gaps.stale_values_round(data, window=3, mark='end')
-    )
-
-
-def test_stale_values_round_smaller_window():
-    """Decreasing window size includes shorter spans of repeated values."""
-    data = pd.Series(
-        [1, 2, 2, 2, 2, 3, 3, 4, 4, 4, 5, 6]
-    )
-    assert_series_equal(
-        gaps.stale_values_round(data, window=3, mark='all'),
-        (data == 2) | (data == 4)
-    )
diff --git a/pvanalytics/tests/util/test_fit.py b/pvanalytics/tests/util/test_fit.py
new file mode 100644
index 0000000..3a7dacf
--- /dev/null
+++ b/pvanalytics/tests/util/test_fit.py
@@ -0,0 +1,9 @@
+"""Tests for curve fitting functions"""
+import pytest
+from pvanalytics.util import _fit
+
+
+def test_fit_not_datetime_raises_error(quadratic):
+    """Trying to fit to an integer index raises a ValueError."""
+    with pytest.raises(ValueError):
+        _fit.quadratic(quadratic)