Spectral mismatch calculation code

docs/sphinx/source/whatsnew/v0.9.2.rst

@@ -17,6 +17,9 @@ Enhancements
 * Improve error message about uneven time intervals for
   :py:func:`~pvlib.clearsky.detect_clearsky` and :py:func:`~pvlib.temperature.prilliman`
   (:issue:`1476`, :pull:`1490`)
+* New module and function to calculate spectral mismatch
+  :py:func:`~pvlib.spetrum.mismatch.calc_spectral_mismatch`
+  (:issue:`1523`, :pull:`1524`)
 
 Bug fixes
 ~~~~~~~~~
@@ -59,9 +62,10 @@ Contributors
 * Adam R. Jensen (:ghuser:`AdamRJensen`)
 * Naman Priyadarshi (:ghuser:`Naman-Priyadarshi`)
 * Chencheng Luo (:ghuser:`roger-lcc`)
-* Prajwal Borkar (:ghuser:`PrajwalBorkar`) 
+* Prajwal Borkar (:ghuser:`PrajwalBorkar`)
 * Kevin Anderson (:ghuser:`kanderso-nrel`)
 * Cliff Hansen (:ghuser:`cwhanse`)
 * Jules Chéron (:ghuser:`jules-ch`)
 * Kurt Rhee (:ghuser:`kurt-rhee`)
 * Will Hobbs (:ghuser:`williamhobbs`)
+* Anton Driesse (:ghuser:`adriesse`)

pvlib/spectrum/mismatch.py

@@ -0,0 +1,164 @@
+"""
+The ``mismatch`` module provides functions for spectral mismatch calculations.
+"""
+
+import pvlib
+import numpy as np
+import pandas as pd
+from scipy.interpolate import interp1d
+import os
+
+
+def get_sample_sr(wavelength=None):
+    '''
+    Generate a generic smooth c-si spectral response for tests and experiments.
+
+    Parameters
+    ----------
+    wavelength: 1-D sequence of numeric, optional
+        The wavelengths at which the sample sr should be interpolated. [nm]
+        By default the wavelengths are from 280 to 1200 in 5 nm intervals.
+
+    Notes
+    -----
+    This sr is based on measurements taken on a reference cell.
+    The measured data points have been adjusted by PV Performance Labs so that
+    standard cubic spline interpolation produces a curve without oscillations.
+    '''
+    SR_DATA = np.array([[ 290, 0.00],
+                        [ 350, 0.27],
+                        [ 400, 0.37],
+                        [ 500, 0.52],
+                        [ 650, 0.71],
+                        [ 800, 0.88],
+                        [ 900, 0.97],
+                        [ 957, 1.00],
+                        [1000, 0.93],
+                        [1050, 0.58],
+                        [1100, 0.21],
+                        [1150, 0.05],
+                        [1190, 0.00]]).transpose()
+
+    if wavelength is None:
+        resolution = 5.0
+        wavelength = np.arange(280, 1200 + resolution, resolution)
+
+    interpolator = interp1d(SR_DATA[0], SR_DATA[1],
+                            kind='cubic',
+                            bounds_error=False,
+                            fill_value=0.0,
+                            copy=False,
+                            assume_sorted=True)
+
+    sr = pd.Series(data=interpolator(wavelength), index=wavelength)
+
+    sr.index.name = 'wavelength'
+    sr.name = 'sr'
+
+    return sr
+
+
+def get_am15g(wavelength=None):
+    '''
+    Read the ASTM G173-03 AM1.5 global tilted spectrum, optionally interpolated
+    to the specified wavelength(s).
+
+    Parameters
+    ----------
+    wavelength: 1-D sequence of numeric, optional
+        The wavelengths at which the spectrum should be interpolated. [nm]
+        By default the 2002 wavelengths of the standard are returned.
+
+    Notes
+    -----
+    This function uses linear interpolation.  If the reference spectrum is too
+    coarsely interpolated, its integral may deviate from 1000.37 W/m2.
+
+    More information about reference spectra is found here:
+    https://www.nrel.gov/grid/solar-resource/spectra-am1.5.html
+
+    The original file containing the reference spectra can be found here:
+    https://www.nrel.gov/grid/solar-resource/assets/data/astmg173.xls
+    '''
+
+    pvlib_path = pvlib.__path__[0]
+    filepath = os.path.join(pvlib_path, 'data', 'astmg173.xls')
+
+    g173 = pd.read_excel(filepath, index_col=0, skiprows=1)
+    am15g = g173['Global tilt  W*m-2*nm-1']
+
+    if wavelength is not None:
+        interpolator = interp1d(am15g.index, am15g,
+                                kind='linear',
+                                bounds_error=False,
+                                fill_value=0.0,
+                                copy=False,
+                                assume_sorted=True)
+
+        am15g = pd.Series(data=interpolator(wavelength), index=wavelength)
+
+    am15g.index.name = 'wavelength'
+    am15g.name = 'am15g'
+
+    return am15g
+
+
+def calc_spectral_mismatch(sr, e_sun, e_ref=None):
+    """
+    Calculate the spectral mismatch under a given measured spectrum with
+    respect to a reference spectrum.
+
+    Parameters
+    ----------
+    sr: pandas.Series
+        The spectral response of one (photovoltaic) device.
+        The index of the Series must contain wavelength values in nm.
+    e_sun: pandas.DataFrame or pandase.Series
+        One or more measured irradiance spectra in a pandas.DataFrame
+        having wavelength in nm as column index.  A single spectrum may be
+        be given as a pandas.Series having wavelength in nm as index.
+    e_ref: pandas.Series, optional
+        The reference spectrum to use for the mismatch calculation.
+        The index of the Series must contain wavelength values in nm.
+        The default is the ASTM G173-03 global tilted spectrum.
+
+    Returns
+    -------
+    smm: pandas.Series or float if a single measured spectrum is provided.
+
+    Notes
+    -----
+    If the default reference spectrum is used it is linearly interpolated
+    to the wavelengths of the measured spectrum.  To achieve alternate
+    behavior e_ref can be transformed before calling this function and
+    provided as an argument.
+
+    The spectral response is linearly interpolated to the wavelengths of the
+    spectrum with which is it multiplied internally (e_sun and e_ref). To
+    achieve alternate behavior the sr can be transformed before calling
+    this function.
+    """
+    # get the reference spectrum at wavelengths matching the measured spectra
+    if e_ref is None:
+        e_ref = get_am15g(wavelength=e_sun.T.index)
+
+    # interpolate the sr at the wavelengths of the spectra
+    # reference spectrum wavelengths may differ if e_ref is from caller
+    sr_sun = np.interp(e_sun.T.index, sr.index, sr, left=0.0, right=0.0)
+    sr_ref = np.interp(e_ref.T.index, sr.index, sr, left=0.0, right=0.0)
+
+    # a helper function to make usable fraction calculations more readable
+    def integrate(e):
+        return np.trapz(e, x=e.T.index, axis=-1)
+
+    # calculate usable fractions
+    uf_sun = integrate(e_sun * sr_sun) / integrate(e_sun)
+    uf_ref = integrate(e_ref * sr_ref) / integrate(e_ref)
+
+    # mismatch is the ratio or quotient of the usable fractions
+    smm = uf_sun / uf_ref
+
+    if isinstance(e_sun, pd.DataFrame):
+        smm = pd.Series(smm, index=e_sun.index)
+
+    return smm

pvlib/tests/test_spectrum.py

@@ -1,13 +1,13 @@
 import pytest
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_approx_equal, assert_equal
 import pandas as pd
 import numpy as np
 from pvlib import spectrum
+from pvlib.spectrum import mismatch
 from .conftest import DATA_DIR
 
 SPECTRL2_TEST_DATA = DATA_DIR / 'spectrl2_example_spectra.csv'
 
-
 @pytest.fixture
 def spectrl2_data():
     # reference spectra generated with solar_utils==0.3
@@ -34,9 +34,9 @@ def spectrl2_data():
         'aerosol_turbidity_500nm': 0.1,
         'dayofyear': 75
     }
-    df = pd.read_csv(SPECTRL2_TEST_DATA)
+    df = pd.read_csv(SPECTRL2_TEST_DATA, index_col=0)
     # convert um to nm
-    df['wavelength'] *= 1000
+    df['wavelength'] = np.round(df['wavelength'] * 1000, 1)
     df[['specdif', 'specdir', 'specetr', 'specglo']] /= 1000
     return kwargs, df
 
@@ -106,3 +106,68 @@ def test_aoi_gt_90(spectrl2_data):
     for key in ['poa_direct', 'poa_global']:
         message = f'{key} contains negative values for aoi>90'
         assert np.all(spectra[key] >= 0), message
+
+
+def test_get_sample_sr():
+    # test that the sample sr is read and interpolated correctly
+    sr = mismatch.get_sample_sr()
+    assert_equal(len(sr), 185)
+    assert_equal(np.sum(sr.index), 136900)
+    assert_approx_equal(np.sum(sr), 107.6027)
+
+    wavelength = [270, 850, 957, 1200, 4001]
+    expected = [0.0, 0.92732, 1.0, 0.0, 0.0]
+
+    sr = mismatch.get_sample_sr(wavelength)
+    assert_equal(len(sr), len(wavelength))
+    assert_allclose(sr, expected, rtol=1e-5)
+
+
+def test_get_am15g():
+    # test that the reference spectrum is read and interpolated correctly
+    e = mismatch.get_am15g()
+    assert_equal(len(e), 2002)
+    assert_equal(np.sum(e.index), 2761442)
+    assert_approx_equal(np.sum(e), 1002.88, significant=6)
+
+    wavelength = [270, 850, 957, 1200, 4001]
+    expected = [0.0, 0.893720, 0.270670, 0.448250, 0.0]
+
+    e = mismatch.get_am15g(wavelength)
+    assert_equal(len(e), len(wavelength))
+    assert_allclose(e, expected, rtol=1e-6)
+
+
+def test_calc_mismatch(spectrl2_data):
+    # test that the mismatch is calculated correctly with
+    # - default and custom reference sepctrum
+    # - single or multiple sun spectra
+
+    # sample data
+    _, e_sun = spectrl2_data
+    e_sun = e_sun.set_index('wavelength')
+    e_sun = e_sun.transpose()
+
+    e_ref = mismatch.get_am15g()
+    sr = mismatch.get_sample_sr()
+
+    # test with single sun spectrum, same as ref spectrum
+    mm = mismatch.calc_spectral_mismatch(sr, e_sun=e_ref)
+    assert_approx_equal(mm, 1.0, significant=6)
+
+    # test with single sun spectrum
+    mm = mismatch.calc_spectral_mismatch(sr, e_sun=e_sun.loc['specglo'])
+    assert_approx_equal(mm, 0.992393, significant=6)
+
+    # test with single sun spectrum, also used as reference spectrum
+    mm = mismatch.calc_spectral_mismatch(sr,
+                                         e_sun=e_sun.loc['specglo'],
+                                         e_ref=e_sun.loc['specglo'])
+    assert_approx_equal(mm, 1.0, significant=6)
+
+    # test with multiple sun spectra
+    expected = [0.973016, 0.995571, 0.899784, 0.992393]
+
+    mm = mismatch.calc_spectral_mismatch(sr, e_sun=e_sun)
+    assert mm.index is e_sun.index
+    assert_allclose(mm, expected, rtol=1e-6)