calculate Linke turbidity factor from AOD, precip water and airmass

[mikofski]

Given aerosol optical depth (AOD), Angstrom alpha and precipitable water (AKA total column water vapor) calculate Linke turbidity factor using Molineaux (1998), Bird-Hulstrom (1980), Berk (1996) and Kasten (1980).

1. calculate Angstrom turbidity alpha exponent if not known, from AOD at two wavelengths, lambda1 and lambda2:

   alpha0 = -log(aod1 / aod2) / log(lambda1 / lambda2)

   Example with lambda1 = 1240nm and lambda2 = 550nm

   alpha0 = -log(aod1240nm / aod550nm) / log(1240 / 550)

2. get aod at 700[nm] from alpha

   aod700 = aod550* ((700 / 550) ^ (-alpha0))

3. From numerically integrated spectral simulations done with Modtran (Berk, 1996), Molineaux (1998) obtained for the broadband optical depth of a clean and dry atmopshere (fictious atmosphere that comprises only the effects of Rayleigh scattering and absorption by the atmosphere gases other than the water vapor) the following expression where M is airmass.

   delta_cda = - 0.101 + 0.235 * M ^ (-0.16)

4. The broadband water vapor optical depth where pwat is the precipitable water vapor content of the atmosphere in [cm]. The precision of these fits is better than 1% when compared with Modtran simulations in the range 1 < M < 6 and 0 < w < 5 cm.

   delta_w = 0.112 * M ^ (-0.55) * pwat ^ (0.34);

5. Aerosol

   either using (Molineaux 1998)

   delta_a = aod700

   or using (Bird-Hulstrom 1980)

   delta_a = 0.27583*aod380 + 0.35*aod500

6. Using the Kasten pyrheliometric formula (1980), the Linke turbidity at M

       TL = -(9.4 + 0.9*M) * log(exp(-M * (delta_cda + delta_w + delta_a))) / M

This derivation was developed in collaboration with Armel Oumbe @aoumbe

references:

1. A new airmass independent formulation for the Linke turbidity coefficient, Ineichen, Perez
2. Equivalence of pyrheliometric and monochromatic aerosol optical depth at a single key wavelength, Mlineaux, Ineichen and O'Neill

[mikofski]

References: molineaux-ineichen 1998 - Equivalence of pyrheliometric and aerosol optical depth (ao).pdf new airmass for Linke-02.pdf

[mikofski]

@wholmgren there seem to be several derivations of Linke turbidity. Is this worth submitting?

[wholmgren]

Sorry, I don't fully understand your proposal. Are you proposing to add functions to implement some of the equations above? I actually meant to add Bird-Hulstrom 1980 to atmosphere.py when I added the simplified solis model a few months ago.

[mikofski]

@wholmgren these formulas together allow you to convert AOD and Pwat measurements directly to Linke turbidity for use in Ineichen clearsky, which is evidently the only place it's used. Maybe it's related to #101.

For example if you had AOD and Pwat from AERONET and MERRA (or any source), then you could calculate Linke turbidity factors and use them in the clearsky model.

[wholmgren]

I do think that could be useful. My understanding is that using your AOD and PW data with simplified solis is more accurate, however, your proposed functionality would be necessary to actually prove that in pvlib.

[adriesse]

I think it would be really nice to have these cleanly implemented so that various comparisons can be done.

[cwhanse]

@mikofski I'd like to see this function added.

[mikofski]

I struggled with testing, so I compared the broadband AOD measurements from TMY3 from Melbourne, FL with the lookup Linke turbidity function and the output was within an order of magnitude. The figure below shows that Linke turbidity calculated with lambda set to 500[nm] and the daily averages differ by about 10%.

But the TMY3 manual says that the reported AOD is broadband so if I use lambda of 700[nm] then the calculated values are much closer, but the trends deviated in autumn.

I think another test would be to compare it with ECMWF AOD/Pwat data in which the wavelength of the AOD measurement is also reported, and there are measurements at five different bands.

Here's a patch. @wholmgren and @cwhanse, not sure where you want me to put it.

linke_turb_from_aod.py.txt

"""Linke turbidity factor calculated from AOD, Pwat and AM"""

from datetime import datetime
import pandas as pd
import numpy as np
import pvlib
from solar_utils import *
from matplotlib import pyplot as plt
import seaborn as sns

plt.ion()
sns.set_context(rc={'figure.figsize': (12, 8)})

def kasten_80lt(aod, am, pwat, alpha0=1.14, method='Molineaux'):
    """
    Calculate Linke turbidity factor using Kasten pyrheliometric formula (1980).

    :param aod: aerosol optical depth table or value at 500
    :param am: airmass, pressure corrected in atmospheres
    :param pwat: precipitable water or total column water vapor in centimeters
    :param alpha0: Angstrom turbidity alpha exponent, default is 1.14
    :param method: Molineaux (default) or Bird-Huldstrom
    :return: Linke turbidity

    Aerosol optical depth can be given as a list of tuples with wavelength in
    nanometers as the first item in each tuple and values as AOD as the second
    item. The list must be in ascending order by wavelength. If ``aod`` is given
    as a sequence of floats or as a float then a wavelength of 500[nm] will be
    used and alpha will default to 1.14, unless alpha is also given. Otherwise
    alpha is calculated from the given wavelength and aod.

    Method can be either ``'Molineaux'`` or ``'Bird-Huldstrom'``. Airmass less
    than 1 or greater than 6 will return ``NaN``. Precipitable water less than zero
    or greater than 5[cm] will also return ``NaN``.
    """
    # calculate Angstrom turbidity alpha exponent if not known, from AOD at two
    # wavelengths, lambda1 and lambda2
    alpha = []
    try:
        # xrange(0) means iterate zero times, xrange(negative) == xrange(0)
        for idx in xrange(len(aod) - 1):
            lambda1, aod1 = aod[idx]
            lambda2, aod2 = aod[idx + 1]
            alpha.append(-np.log(aod1 / aod2) / np.log(lambda1 / lambda2))
    except TypeError:
        # case 1: aod is a float, so len(aod) raises TypeError
        # case 2: aod is an array of float, so (lambda1, aod1) = aod[idx] raises
        # TypeError
        aod = [(500.0, aod)]
    else:
        # case 3: len(aod) == 1, then alpha == []
        if len(alpha) > 1:
            alpha0 = alpha
    # make sure alpha can be indexed
    try:
        alpha = list(alpha0)
    except TypeError:
        alpha = [alpha0]
    # make sure aod has lambda
    try:
        # case 3: len(aod) == 1 and aod == [aod]
        lambda1, aod1 = zip(*aod)
    except TypeError:
        aod = [(500.0, aod)]
    # From numerically integrated spectral simulations done with Modtran (Berk,
    # 1996), Molineaux (1998) obtained for the broadband optical depth of a
    # clean and dry atmopshere (fictious atmosphere that comprises only the
    # effects of Rayleigh scattering and absorption by the atmosphere gases
    # other than the water vapor) the following expression where am is airmass.
    delta_cda = -0.101 + 0.235 * am ** (-0.16)
    # The broadband water vapor optical depth where pwat is the precipitable
    # water vapor content of the atmosphere in [cm]. The precision of these fits
    # is better than 1% when compared with Modtran simulations in the range
    # 1 < am < 6 and 0 < pwat < 5 cm.
    delta_w = 0.112 * am ** (-0.55) * pwat ** (0.34)
    if method == 'Molineaux':
        # get aod at 700[nm] from alpha for Molineaux (1998)
        delta_a = get_aod_at_lambda(aod, alpha)
    else:
        # using (Bird-Hulstrom 1980)
        aod380 = get_aod_at_lambda(aod, alpha, 380.0)
        aod500 = get_aod_at_lambda(aod, alpha, 500.0)
        delta_a = 0.27583 * aod380 + 0.35 * aod500
    # the Linke turbidity at am using the Kasten pyrheliometric formula (1980)
    lt = -(9.4 + 0.9 * am) * np.log(
        np.exp(-am * (delta_cda + delta_w + delta_a))
    ) / am
    # filter out of extrapolated values
    filter = (amp < 1.0) | (amp > 6.0) | (pwat_cm < 0) | (pwat_cm > 5.0)
    lt[filter] = np.nan  # set out of range to NaN
    return lt

def get_aod_at_lambda(aod, alpha, lambda0=700.0):
    """
    Get AOD at specified wavelenth.

    :param aod: sequence of (wavelength, aod) in ascending order by wavelength
    :param alpha: sequence of Angstrom alpha corresponding to wavelength in aod
    :param lambda0: desired wavelength in nanometers, defaults to 700[nm]
    """
    lambda1, aod = zip(*aod)
    # lambda0 is between (idx - 1) and idx
    idx = np.searchsorted(lambda1, lambda0)
    # unless idx is zero
    if idx == 0:
        idx = 1
    return aod[idx - 1] * ((lambda0 / lambda1[idx - 1]) ** (-alpha[idx - 1]))

if __name__ == '__main__':
    melbourne_fl = pvlib.tmy.readtmy3('NREL/TMY3/722040TYA.CSV')
    aod500 = melbourne_fl[0]['AOD']
    pwat_cm = melbourne_fl[0]['Pwat']
    press_mbar = melbourne_fl[0]['Pressure']
    dry_temp = melbourne_fl[0]['DryBulb']
    timestamps = melbourne_fl[0].index
    location = (melbourne_fl[1]['latitude'],
                melbourne_fl[1]['longitude'],
                melbourne_fl[1]['TZ'])
    _, airmass = zip(*[solposAM(
        location, d.timetuple()[:6], (press_mbar.loc[d], dry_temp.loc[d])
    ) for d in timestamps])
    _, amp = zip(*airmass)
    amp = np.array(amp)
    filter = amp < 0
    amp[filter] = np.nan
    t = pd.DatetimeIndex([datetime.replace(d, year=2016) for d in timestamps])
    lt_molineaux = kasten_80lt(aod=aod500, am=amp, pwat=pwat_cm)
    lt_molineaux.index = t
    lt_bird_huldstrom = kasten_80lt(aod=aod500, am=amp, pwat=pwat_cm,
                                    method='Bird-Huldstrom')
    lt_bird_huldstrom.index = t
    pvlib.clearsky.lookup_linke_turbidity(t, *location[:2]).plot()
    lt_molineaux.resample('D').mean().plot()
    lt_bird_huldstrom.resample('D').mean().plot()
    title = ['Linke turbidity factor comparison at Melbourne, FL (722040TYA),',
             'calculated using Kasten Pyrheliometric formula']
    plt.title(' '.join(title))
    plt.ylabel('Linke turbidity factor')
    plt.legend(['Linke Turbidity', 'Molineaux', 'Bird-Huldstrom'])

To get the 700[nm] plot I change L129:

lt_molineaux = kasten_80lt(aod=[(700.0, aod500)], am=amp, pwat=pwat_cm)

and L131-132:

lt_bird_huldstrom = kasten_80lt(aod=[(700.0, aod500)], am=amp, pwat=pwat_cm,
                                method='Bird-Huldstrom')

Let me know what other ideas you have.

[mikofski]

Another reference from Ineichen that outlines these steps pretty well:

• Conversion function between the Linke turbidity and the atmospheric water vapor and aerosol content, P. Ineichecn, Solar Energy 82 (2008), p 1095-1097

Ineichen_SolarEnergy(82)_Conversion-of-Linke-turbidity_2008.pdf

Ineichen also outlines an improved method by Mueller (2004) and Mayer (1997) derived from the simplified solis scheme:

where p is the atmospheric pressure at the considered altitude, po the sea level atmospheric pressure and aod550 the aerosol optical depth at 550 nm

Should be pretty easy to add this as an option to the existing kasten_80lt() method.

[wholmgren]

Nice. Can you really use a wavelength of 700 as a substitute for broadband, or was that just to prove that the wavelength matters? Good idea to test with ECMWF data, though I don't know if that's necessary to do before adding this to the library.

I can see arguments for putting this kind of thing in either atmosphere.py or clearsky.py. I don't have a recommendation at this point.

+1 on the second TL formulation.

It might be best to write separate (possibly private) functions for each of the steps you initially outlined, then another function or two that ties them all together. We can also defer discussion of implementation details to a PR.

[cwhanse]

https://www.ncbi.nlm.nih.gov/pubmed/18301519

There’s an observed equivalence between broadband AOD and AOD at 700nm, so I wouldn’t be surprised if 700nm is also a good surrogate for TL.

We’ve seen a similar relationship between spectral irradiance and short-circuit current. Generally speaking, you don’t have to measure the whole spectrum to get an appropriate spectral mismatch factor, only a few wavelengths, which look to be fairly independent of cell technology (Matt Lave’s PVSC42 papers). I think what we’re seeing is the same phenomenon as in the Molineaux paper, but interpreted through a PV module rather than an irradiance instrument.

Cliff

[wholmgren]

Thanks, Cliff. That all makes sense and once again I should have read Mark's original post more carefully.

[mikofski]

@BenBourne can you take a look at this?

[mikofski]

@wholmgren @cwhanse this is done maybe? Now rush to do the abstract, :frowning_face: sorry so slow