Spectral mismatch calculation code

[adriesse]

• [ ] Closes #1523 23
• [ ] I am familiar with the contributing guidelines
• [ ] Tests added
• [ ] Updates entries in docs/sphinx/source/reference for API changes.
• [ ] Adds description and name entries in the appropriate "what's new" file in docs/sphinx/source/whatsnew for all changes. Includes link to the GitHub Issue with :issue:`num` or this Pull Request with :pull:`num`. Includes contributor name and/or GitHub username (link with :ghuser:`user`).
• [ ] New code is fully documented. Includes numpydoc compliant docstrings, examples, and comments where necessary.
• [ ] Pull request is nearly complete and ready for detailed review.
• [ ] Maintainer: Appropriate GitHub Labels (including remote-data) and Milestone are assigned to the Pull Request and linked Issue.

[adriesse]

Empty module: mismatch. I anticipate some more broadly usable functions could find a home here over time, so a more generic name might be considered.

[adriesse]

So we will probably have some interesting discussion on the type and shape of the outputs. While working on this I noticed that spectrl2 has wavelength as the first dimension. When converted to a dataframe this make wavelength the row index. I would prefer to have wavelength as the last dimension, which I find works more naturally.

[cwhanse]

What is the envisioned scope of pvlib/spectrum? Among other purposes, it seems like a natural place to collect the functions for spectral mismatch models that are currently in pvsystem and atmosphere. And add some code to help use pvlib with SMARTS or other such tools.

[adriesse]

Yes, I think it would make some sense to move those two functions. At the moment I can't imagine spectrum growing a lot more.

[adriesse]

I've been thinking about numpy friendliness here, and I am tending toward making these functions pandas only. The wavelength values are an integral part of the data so it makes a lot of sense to me to carry them around as row or column index.

[adriesse]

Ok, I think it works now. :)

[adriesse]

I think this is ready for review. Not sure whether I cause all the automated test failure--I hope not.

[kandersolar]

Not sure whether I cause all the automated test failure--I hope not.

At least some of the failures are because xlrd is needed to parse the excel file but it's not installed as a pvlib dependency. Maybe better to convert the file to CSV? Doesn't seem like a big loss to drop the Chart1 sheet.

[adriesse]

Should I try attempt to change docs/sphinx/source/reference, or leave that to someone more familiar with this?

[kandersolar]

Certainly not opposed to housing mismatch functions in pvlib.spectrum.mismatch, but should an alias be set up in pvlib/spectrum/__init__.py so that the functions are also accessible as e.g. pvlib.spectrum.calc_spectral_mismatch (in addition to pvlib.spectrum.mismatch.calc_spectral_mismatch)? We did that for spectrl2.

I wonder if spectrl2 ought to have been pandas-only as well.

Should I try attempt to change docs/sphinx/source/reference

Please do! I think adding a couple entries at the bottom of docs/sphinx/source/reference/effects_on_pv_system_output.rst would do.

[adriesse]

Ok, I made those adjustments. I'm afraid I don't understand why github still sees a conflict in the whatsnew though.

[adriesse]

Can someone suggest how to remove docs/sphinx/source/whatsnew/v0.9.2.rst from my PR?

[kandersolar]

I think (hope?) the 0.9.2 whatsnew changes will go away by resolving the merge conflict (Resolve conflicts button here on GH, or can be done locally). When resolving the conflict I would keep the master version in each case.

[adriesse]

I feel a bit silly here, but I don't know where to click to actually resolve the conflicts it shows.

[kandersolar]

Oh! This guide might be useful: https://docs.github.com/en/pull-requests/collaborating-with-pull-requests/addressing-merge-conflicts/resolving-a-merge-conflict-on-github

tl;dr: resolving the conflicts is done mostly via keyboard rather than mouse. In an editable text box, it presents you with the two versions (this PR's version and master's version) of each conflict, demarcated by === and >>> lines. Pick the version you want to keep, delete the other, delete the demarcation lines, and commit the changes. Doing it locally is basically the same except you use your own text editor instead of github's in-browser editor.

[adriesse]

Thanks, I had no idea that box was editable!!

[markcampanelli]

@adriesse I have just started to take a look at this. I hit some snags trying to pull the branch, but I think that's on my end. Hopefully I'll have that sorted soon.

[adriesse]

Ok, I've resolved a bunch of stuff. There were some that perhaps would have been more appropriately marked as resolved by the originators, but I don't know whether that is possible or common practice.

[markcampanelli]

@adriesse I added a new comment under the resolved comment for the usable fraction calculations that may not show up as new for you.

[markcampanelli]

Using the example data for the spectral radiance of the sun and the spectral responsivity of a Si reference device, the error in the SMM when truncating both integrals appears to small, on the order of 0.1%. One does not, however, want to only truncate the incomplete sun integral.

e_ref totals to 1.7 um, 4 um, and infinity (W/m^2): 942.88 997.47 1000.0
e_sun totals to 1.7 um, 4 um, and infinity (W/m^2): 668.028852782929 708.4855034837434 710.2825182549283
no_trunc smm: 0.9913868433999364
1.7 um single truncation smm: 1.054093308636961
4 um single truncation smm: 0.9939014139772988
1.7 um double truncation smm: 0.9902678450884946
4 um double truncation smm: 0.9923965546943567
1.7 um single truncation error in smm: 6.325125822930477 %
4 um single truncation error in smm: 0.25364171353523535 %
1.7 um double truncation error in smm: -0.11287201548935144 %
4 um double truncation error in smm: 0.10184836536235586 %

Here is the code that produced these results:

import numpy as np
from pvlib import spectrum
from pvlib.tests import test_spectrum

e_ref = spectrum.get_am15g()  # 0.28 to 4 micron, scaled to "match" IEC 60904-3 Ed. 2 (2008)
e_ref_integral_to_1p7_micron = 942.88  # W/m^2, from IEC 60904-3 Ed. 2 (2008)
e_ref_integral_to_4_micron = 997.47  # W/m^2, from IEC 60904-3 Ed. 2 (2008)
e_ref_integral_to_infinity = 1000.  # W/m^2, from IEC 60904-3 Ed. 2 (2008)

# Get example solar spectrum from test fixture.
_, e_sun = test_spectrum.spectrl2_data.__pytest_wrapped__.obj()  # 0.3 to 4 micron
e_sun = e_sun.set_index('wavelength')
e_sun = e_sun.transpose().loc['specglo']

idx_to_1p7_micron = e_sun.T.index <= 1700
e_sun_integral_to_1p7_micron = np.trapz(e_sun[idx_to_1p7_micron], x=e_sun.T.index[idx_to_1p7_micron], axis=-1)
e_sun_integral_to_4_micron = np.trapz(e_sun, x=e_sun.T.index, axis=-1)  # 708.4855034837434 W/m^2
scale_factor_estimate = e_ref_integral_to_infinity / e_ref_integral_to_4_micron
e_sun_integral_to_infinity = scale_factor_estimate * e_sun_integral_to_4_micron

sr = spectrum.get_example_spectral_response()

no_trunc = spectrum.calc_spectral_mismatch_field(sr, e_sun, e_sun_integral_to_infinity, e_ref=e_ref, e_ref_tot=e_ref_integral_to_infinity)
single_trunc_1p7_micron = spectrum.calc_spectral_mismatch_field(sr, e_sun[idx_to_1p7_micron], None, e_ref=e_ref, e_ref_tot=e_ref_integral_to_infinity)
single_trunc_4_micron_ = spectrum.calc_spectral_mismatch_field(sr, e_sun, None, e_ref=e_ref, e_ref_tot=e_ref_integral_to_infinity)
double_trunc_1p7_micron = spectrum.calc_spectral_mismatch_field(sr, e_sun[idx_to_1p7_micron], None)
double_trunc_4_micron = spectrum.calc_spectral_mismatch_field(sr, e_sun, None)

print("e_ref totals to 1.7 um, 4 um, and infinity (W/m^2):", e_ref_integral_to_1p7_micron, e_ref_integral_to_4_micron, e_ref_integral_to_infinity)
print("e_sun totals to 1.7 um, 4 um, and infinity (W/m^2):", e_sun_integral_to_1p7_micron, e_sun_integral_to_4_micron, e_sun_integral_to_infinity)
print("no_trunc smm:", no_trunc)
print("1.7 um single truncation smm:", single_trunc_1p7_micron)
print("4 um single truncation smm:", single_trunc_4_micron_)
print("1.7 um double truncation smm:", double_trunc_1p7_micron)
print("4 um double truncation smm:", double_trunc_4_micron)
print("1.7 um single truncation error in smm:", 100*(single_trunc_1p7_micron / no_trunc - 1), "%")
print("4 um single truncation error in smm:", 100*(single_trunc_4_micron_ / no_trunc - 1), "%")
print("1.7 um double truncation error in smm:", 100*(double_trunc_1p7_micron / no_trunc - 1), "%")
print("4 um double truncation error in smm:", 100*(double_trunc_4_micron / no_trunc - 1), "%")

I had to change the source code in this PR in order to calculate the "no_trunc" result (doc strings omitted):

def get_am15g(wavelength=None):
    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Aug. 2022
    # Changed by Mark Campanelli for testing purposes.

    pvlib_path = pvlib.__path__[0]
    filepath = os.path.join(pvlib_path, 'data', 'astm_g173_am15g.csv')

    am15g = pd.read_csv(filepath, index_col=0).squeeze()  # 0.28 to 4 microns
    am15g_integral_to_4_micron = np.trapz(am15g, x=am15g.index, axis=-1)  # W/m^2
    iec_integral_to_4_micron = 997.47  # W/m^2, from IEC 60904-3 Ed. 2 (2008)
    am15g_scaled_to_iec = iec_integral_to_4_micron / am15g_integral_to_4_micron * am15g

    if wavelength is not None:
        interpolator = interp1d(am15g_scaled_to_iec.index, am15g_scaled_to_iec,
                                kind='linear',
                                bounds_error=False,
                                fill_value=0.0,
                                copy=False,
                                assume_sorted=True)

        am15g_scaled_to_iec = pd.Series(data=interpolator(wavelength), index=wavelength)

    am15g_scaled_to_iec.index.name = 'wavelength'
    am15g_scaled_to_iec.name = 'am15g'

    return am15g_scaled_to_iec


def calc_spectral_mismatch_field(sr, e_sun, e_sun_tot, e_ref=None, e_ref_tot=None):
    # Contributed by Anton Driesse (@adriesse), PV Performance Labs. Aug. 2022.
    # Changed by Mark Campanelli for testing purposes.

    # get the reference spectrum at wavelengths matching the measured spectra
    if e_ref is None:
        e_ref = get_am15g(wavelength=e_sun.T.index)

    # a helper function to make usable fraction calculations more readable
    def integrate(e):
        return np.trapz(e, x=e.T.index, axis=-1)

    if e_sun_tot is None:
        # Use integrated spectral irradiance, typically an approximation
        # that misses "tail" of distribution.
        e_sun_tot = integrate(e_sun)

    if e_ref_tot is None:
        # Use integrated spectral irradiance, typically an approximation
        # that misses "tail" of distribution.
        e_ref_tot = integrate(e_ref)

    # 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)

    # calculate usable fractions
    uf_sun = integrate(e_sun * sr_sun) / e_sun_tot
    uf_ref = integrate(e_ref * sr_ref) / e_ref_tot

    # 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

I am guessing that double truncation is not a bad strategy in many (most?) cases, but I will point out that an alternate strategy (ref. Carl Osterwald) is to estimate the total sun spectrum by integrating the tail of the reference spectrum after scaling it to "match" some portion of the sun spectrum. I suppose we would have to do a bigger study to fairly compare the two approaches.

I have some additional thoughts now about how we could (optionally) include the total irradiance value(s) but still try to keep the function interface from being too confusing. (This idea is different than the code changes I used above.) That said, I would also be ok if the double truncation approximation was more clearly stated in the docstring and perhaps a note saying that the code is not appropriate for computations such as primary reference cell calibrations such as ASTM E1125. (It turns out that ASTM E1125 is where the differences between ASTM G173-03 and IEC 60904-3 is called out and handled.)

[adriesse]

Thanks very much, @markcampanelli for checking this out so thoroughly and confirming at the very least that I'm not terribly mistaken! :) The double truncation is more or less imposed when one of the inputs is pre-truncated in field measurements. One additional easy test is to image the sun spectrum being equal to the reference spectrum, so you want to get a factor of 1.0. If the spectroradiometer only measures to 1700 nm, then the only way to get 1.0 is to truncate the reference spectrum to 1700 as well.

I've added a paragraph to the docstring to explain the truncation.

[adriesse]

@nicorie would you be interested in having a look?

[adriesse]

To all: if you feel that this PR is premature, don't feel like you have to merge it. It's not a lot of code and I can easily just pack it into an example notebook.

[kandersolar]

Well this PR does have three approving reviews, which is more than most pvlib PRs :) I haven't been following the discussion with @markcampanelli closely, but I have the sense that at this point whatever limitations there are in this implementation are noted in the docstring, and of course we can always make additional improvements later if we want. Anyone please correct me if I'm mistaken there. Otherwise let's plan to merge this tomorrow?

[markcampanelli]

@adriesse Would it be possible to include an example in this PR of how this code is used for a specific use case or cases?

If I were to try to use it for an IEC 61853-3 like calculation with a 1.7 micron cutoff, then I could get tripped up about which total irradiance to multiply the SMM by to get the equivalent irradiance under the standard spectrum used in the 61853-1 matrix (i.e., using the double truncation approximation).

I am curious if pvlib workflows supports any use cases. (Not that the answer to this has to be positive, the code can be useful regardless.)

[adriesse]

@markcampanelli An example will accompany the forthcoming sr collection and spectroradiometer measurements on Duramat.

[mikofski]

There are conflicts with the main branch, maybe what’s new? That must be addressed before we can merge

[adriesse]

I think this needs coordination with the other PR(s), doesn't it?

[cwhanse]

The conflict was created when #1468 was merged. I can resolve if you'd like the help.

[adriesse]

Sure @cwhanse, it seems easiest to do this at the time of merging.