Add pw and spectral correction code from pvlib matlab 1.3

[wholmgren]

This PR is a very quick port of the PVLIB Matlab 1.3 functions pvl_calcPwat and pvl_FSspeccorr. I'm implementing these quickly to add a little more substance to a pvsc abstract.

Problems:

1. The PW function will not return values less than 0.1 cm. What!?
2. The spectral correction function diverges as PW goes to 0. What!?
3. No unit tests yet.
4. Haven't checked that the doc strings compile correctly.
5. The new function and variable names are up for debate. I used pw, but precipitable_water is probably better aligned with our style.

cc @cwhanse @jforbess. ref #2.

import numpy as np
import matplotlib.pyplot as plt
from pvlib.atmosphere import calc_pw, first_solar_spectral_correction

airmasss = np.linspace(1.2, 5)
rhs = np.linspace(20, 100, 5)
pws = calc_pw(25, rhs)
for pw in pws:
    cdte = first_solar_spectral_correction(pw, airmasss, 'CdTe')
    xsi = first_solar_spectral_correction(pw, airmasss, 'xSi')
    plt.plot(airmasss, cdte, 'r-')
    plt.plot(airmasss, xsi, 'b-')

plt.xlabel('Air mass')
plt.ylabel('Spectral mismatch modifier')
plt.title('Effect of Relative Humidity on Spectral Mismatch')
plt.legend(['x-Si', 'CdTe'], loc='lower center')
plt.text(4, .995, 'RH = 100%', color='r')
plt.text(3.6, .92, 'RH = 20%', color='r')
plt.text(2.5, 1.046, 'RH = 100%', color='b')
plt.text(3.5, 1.024, 'RH = 20%', color='b')

[cwhanse]

Will,

I’m not an atmospheric scientist but the limint of 0.1cm on preciptable water seems reasonable to me.

http://www.cgd.ucar.edu/cas/adai/papers/Wang_etal_JGR07_2006JD007529.pdf

I noted the division by Pwat also. But, those functions are what First Solar contributed, so I’d rather leave them as is and attribute them.

Chris Gueymard emailed about the PW function saying there’s more recent work that is considered better. We haven’t looked into writing a companion PW function yet, nor if that AM/Pwat term in the spectral correction model is significant.

Cliff

[wholmgren]

0.1 cm is certainly drier than most places in the US ever get. Cooler sites at higher elevation can occasionally see values below 0.1 cm. Here's GPS IPW measurements from Albuquerque (P034 and ZAB1) and Flagstaff (FST5) over the last couple of months.

Looks to me like there were a few data points below 0.1 cm. Plots from http://gpsmet.noaa.gov/cgi-bin/gnuplots/rti.cgi

The SPC Sounding Climatology Page is also worth looking at.

I'm always skeptical of the extreme points, but there is some QA in the data.

Is there a downside to a limit of 0 in the PW function?

Concerning the FS function, in this case I would support adding an internal coercion to some limiting number to prevent it from blowing up.

[cwhanse]

Good points, thanks for sending the figures. I think the PW limit is there only to avoid blowup in the spectral correction model.

First Solar guys would be more than open to improvements if you have an idea of a reasonable way to prevent /0. There’s no really firm physical limit on that modifier, though, I suppose we could impose something like +/-50% and no one would question it.

[wholmgren]

Sorry for the delay in getting back to this. I think it can now be merged, though I'm definitely open to further discussion and improvements.

Here's what I've changed...

I changed the name of calc_pw to gueymard94_pw. I prefer to start off being more specific with the name, even if it means being different from the matlab library, so that we can more easily add new functions later. I skimmed the reference documentation and verified that the equations are correct.

See below for a couple of figures that show the results gueymard94_pw(temp, rh) calculation, one with the color spanning the full calculation output and one with the color spanning only 0 to 5 cm.

I left the 0.1 lower limit alone for now. I suggest adding an upper limit of ~5 cm.

I updated the FS coefficients according to an email that Josh forwarded to me.

I added some tests for each function.

A link to the docs for this branch:

http://wholmgren-pvlib-python-new.readthedocs.org/en/pwfsspec/modules.html#module-pvlib.atmosphere

[wholmgren]

Forgot to add the FS correction plots... note that the scales are not the same...

[cwhanse]

Renaming to gueymard94_pw is good, I’m sure there are other PW functions so we ought to have a structure that accommodates other models. I’ll copy that in the Matlab files.

FS is still fine tuning that model, I’d wait until after PVSC to publish it in the toolbox.

Cliff

[MLEEFS]

It appears that people have been talking about me (and making really pretty plots that I am going to try to recreate). Now that this work has been published, I am going to work with a colleague at First Solar to submit this through the appropriate GitHub channels.

I now recognize the divergence issue. I didn't realize that Pwat could go below 0.1 cm in terrestrial environments. My natural inclination is to replace Pwat < 0.1 cm with 0.1 cm (or NaN) within the spectral correction function. Some basic SMARTS simulations could be run to quantify bias introduced. M(Pwat=0,AM=1) for CdTe and cSi.

With regards to Pwat above 5 cm, I would rather not filter them out. The M(Pwat,AM) surface is fairly flat for extremely high Pwat values. This is because there is only so much infrared light that can be absorbed. Likewise, I could do a couple of SMARTS simulations to see what happens above Pwat = 5cm,

• Mitch Lee (First Solar)

[wholmgren]

Hi Mitch, thanks for commenting. Just so there's no confusion... did you click on the Files changed tab to see my implementation of the algorithm? I have no problem discarding my python code in favor of yours. Another option is that you can start with my code, make changes, and then make a new pull request. It sounded like you have a git expert on your team so I’ll leave it at that for now.

Here's the code that made those plots...

import numpy as np
import matplotlib.pyplot as plt
from pvlib.atmosphere import first_solar_spectral_correction, gueymard94_pw

temp_air = np.linspace(-10, 40, 51)
relative_humidity = np.linspace(0, 100, 51)

temps_humids = np.meshgrid(temp_air, relative_humidity)

pws = gueymard94_pw(temps_humids[0], temps_humids[1])

fig, ax = plt.subplots(figsize=(12,9), subplot_kw=dict(aspect=0.5))
im = ax.contourf(temps_humids[0], temps_humids[1], pws, 50, cmap=plt.get_cmap('viridis'), vmin=0, vmax=5)
ax.set_ylabel('relative humidity (%)')
ax.set_xlabel('temperature (C)')
fig.colorbar(im, label='pw (cm)')

fig, ax = plt.subplots(figsize=(12,9), subplot_kw=dict(aspect=0.5))
im = ax.contourf(temps_humids[0], temps_humids[1], pws, 50, cmap=plt.get_cmap('viridis'))
ax.set_ylabel('relative humidity (%)')
ax.set_xlabel('temperature (C)')
fig.colorbar(im, label='pw (cm)')

ams = np.linspace(1, 5, 51)
pws = np.linspace(0, 5, 51)
ams, pws = np.meshgrid(ams, pws)

cmap = plt.get_cmap('viridis')
n = 15
vmin = .94
vmax = 1.05

for module_type in ['cdte', 'xsi', 'polysi']:
    mods = first_solar_spectral_correction(pws, ams, module_type=module_type)

    fig, ax = plt.subplots(figsize=(10,7))

    im = ax.contour(ams, pws, mods, n, cmap=cmap, vmin=vmin, vmax=vmax)
    im2 = ax.contourf(ams, pws, mods, n, cmap=cmap, vmin=vmin, vmax=vmax)
    ax.set_ylabel('preciptable water (cm)')
    ax.set_xlabel('airmass')
    ax.clabel(im, colors='k', fmt='%.2f')
    fig.colorbar(im2, label='modifier')
    ax.set_title(module_type)

[MLEEFS]

Hi Will,

I would like to work off of what you have already done. I'm guessing the easiest way for me to do that is for you to merge the work that you have already done. Then, I can make the necessary changes, such that the PVlib python, matches my IEEE PVSC paper (and the PVlib matlab function I submitted today)

sandialabs/MATLAB_PV_LIB#3

• Mitch

[wholmgren]

Yes, I think that's a good way to do it in this case. Merging...