Repository: pvlib/pvlib-python (python)

Issue:
calculate Linke turbidity factor from AOD, precip water and airmass
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)][2] 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][2])

        delta_a = aod700

    or using (Bird-Hulstrom 1980)

        delta_a = 0.27583*aod380 + 0.35*aod500

6. Using the [Kasten pyrheliometric formula (1980)][1], 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

[1]: http://www.asrc.albany.edu/people/faculty/perez/publications/Solar%20Resource%20Assessment%20and%20Modeling/Papers%20on%20Resource%20Assessment%20and%20Satellites/new%20airmass%20for%20Linke-02.pdf
[2]: http://www.cuepe.ch/html/biblio/pdf/molineaux-ineichen%201998%20%20-%20Equivalence%20of%20pyrheliometric%20and%20aerosol%20optical%20depth%20(ao).pdf

references:
1. [_A new airmass independent formulation for the Linke turbidity coefficient_, Ineichen, Perez](http://www.asrc.albany.edu/people/faculty/perez/publications/Solar%20Resource%20Assessment%20and%20Modeling/Papers%20on%20Resource%20Assessment%20and%20Satellites/new%20airmass%20for%20Linke-02.pdf)
2. [**Equivalence of pyrheliometric and monochromatic aerosol optical depth at a single key wavelength**, Mlineaux, Ineichen and O'Neill](http://www.cuepe.ch/html/biblio/pdf/molineaux-ineichen%201998%20%20-%20Equivalence%20of%20pyrheliometric%20and%20aerosol%20optical%20depth%20(ao).pdf)

Pull request:
add linke turbidity formulas module to atmosphere
* closes #269 
* add kasten 1996 pyrheliometric formula for linke turbidity factor
* move atmosphere into package, just for fun, why not?
* make __init__ with same api as before

Signed-off-by: Mark Mikofski <mark.mikofski@sunpowercorp.com>

Patch:
diff --git a/.gitignore b/.gitignore
index 658fe2d..f18eba0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -22,6 +22,7 @@ var/
 Drafts/
 */build/
 build/
+venv/
 
 *.egg-info/
 .installed.cfg
@@ -58,6 +59,8 @@ coverage.xml
 .pydevproject
 .spyderproject
 .idea/
+*.sublime-project
+*.sublime-workspace
 
 # Rope
 .ropeproject
diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst
index 332b1b7..ff34b7c 100644
--- a/docs/sphinx/source/api.rst
+++ b/docs/sphinx/source/api.rst
@@ -87,6 +87,10 @@ Airmass and atmospheric models
    atmosphere.alt2pres
    atmosphere.gueymard94_pw
    atmosphere.first_solar_spectral_correction
+   atmosphere.bird_hulstrom80_aod_bb
+   atmosphere.kasten96_lt
+   atmosphere.angstrom_aod_at_lambda
+   atmosphere.angstrom_alpha
 
 
 Irradiance
diff --git a/docs/sphinx/source/conf.py b/docs/sphinx/source/conf.py
index 546e61c..21d3cd5 100644
--- a/docs/sphinx/source/conf.py
+++ b/docs/sphinx/source/conf.py
@@ -270,7 +270,8 @@ latex_documents = [
 extlinks = {'issue': ('https://github.com/pvlib/pvlib-python/issues/%s',
                       'GH'),
             'wiki': ('https://github.com/pvlib/pvlib-python/wiki/%s',
-                     'wiki ')}
+                     'wiki '),
+            'doi': ('http://dx.doi.org/%s', 'DOI: ')}
 
 # -- Options for manual page output ---------------------------------------
 
diff --git a/docs/sphinx/source/whatsnew/v0.4.4.txt b/docs/sphinx/source/whatsnew/v0.4.4.txt
index 48ff7da..e2be6a7 100644
--- a/docs/sphinx/source/whatsnew/v0.4.4.txt
+++ b/docs/sphinx/source/whatsnew/v0.4.4.txt
@@ -6,6 +6,9 @@ v0.4.4 (January xx, 2017)
 Enhancements
 ~~~~~~~~~~~~
 
+* Added Kasten pyrheliometric formula to calculate Linke turbidity factors with
+  improvements by Ineichen and Perez to extend range of air mass (:issue:`278`)
+
 
 Documentation
 ~~~~~~~~~~~~~
@@ -19,3 +22,4 @@ Contributors
 ~~~~~~~~~~~~
 
 * Will Holmgren
+* Mark Mikofski
diff --git a/pvlib/atmosphere.py b/pvlib/atmosphere.py
index 7ea4b5e..0b90094 100644
--- a/pvlib/atmosphere.py
+++ b/pvlib/atmosphere.py
@@ -292,7 +292,7 @@ def gueymard94_pw(temp_air, relative_humidity):
     ----------
     .. [1] W. M. Keogh and A. W. Blakers, Accurate Measurement, Using Natural
        Sunlight, of Silicon Solar Cells, Prog. in Photovoltaics: Res.
-       and Appl. 2004, vol 12, pp. 1-19 (DOI: 10.1002/pip.517)
+       and Appl. 2004, vol 12, pp. 1-19 (:doi:`10.1002/pip.517`)
 
     .. [2] C. Gueymard, Analysis of Monthly Average Atmospheric Precipitable
        Water and Turbidity in Canada and Northern United States,
@@ -467,3 +467,193 @@ def first_solar_spectral_correction(pw, airmass_absolute, module_type=None,
         coeff[4]*np.sqrt(pw) + coeff[5]*ama/np.sqrt(pw))
 
     return modifier
+
+
+def bird_hulstrom80_aod_bb(aod380, aod500):
+    """
+    Approximate broadband aerosol optical depth.
+
+    Bird and Hulstrom developed a correlation for broadband aerosol optical
+    depth (AOD) using two wavelengths, 380 nm and 500 nm.
+
+    Parameters
+    ----------
+    aod380 : numeric
+        AOD measured at 380 nm
+    aod500 : numeric
+        AOD measured at 500 nm
+
+    Returns
+    -------
+    aod_bb : numeric
+        broadband AOD
+
+    See also
+    --------
+    kasten96_lt
+
+    References
+    ----------
+    [1] Bird and Hulstrom, "Direct Insolation Models" (1980)
+    `SERI/TR-335-344 <http://www.nrel.gov/docs/legosti/old/344.pdf>`_
+
+    [2] R. E. Bird and R. L. Hulstrom, "Review, Evaluation, and Improvement of
+    Direct Irradiance Models", Journal of Solar Energy Engineering 103(3),
+    pp. 182-192 (1981)
+    :doi:`10.1115/1.3266239`
+    """
+    # approximate broadband AOD using (Bird-Hulstrom 1980)
+    return 0.27583 * aod380 + 0.35 * aod500
+
+
+def kasten96_lt(airmass_absolute, precipitable_water, aod_bb):
+    """
+    Calculate Linke turbidity factor using Kasten pyrheliometric formula.
+
+    Note that broadband aerosol optical depth (AOD) can be approximated by AOD
+    measured at 700 nm according to Molineaux [4] . Bird and Hulstrom offer an
+    alternate approximation using AOD measured at 380 nm and 500 nm.
+
+    Based on original implementation by Armel Oumbe.
+
+    .. warning::
+        These calculations are only valid for air mass less than 5 atm and
+        precipitable water less than 5 cm.
+
+    Parameters
+    ----------
+    airmass_absolute : numeric
+        airmass, pressure corrected in atmospheres
+    precipitable_water : numeric
+        precipitable water or total column water vapor in centimeters
+    aod_bb : numeric
+        broadband AOD
+
+    Returns
+    -------
+    lt : numeric
+        Linke turbidity
+
+    See also
+    --------
+    bird_hulstrom80_aod_bb
+    angstrom_aod_at_lambda
+
+    References
+    ----------
+    [1] F. Linke, "Transmissions-Koeffizient und Trubungsfaktor", Beitrage
+    zur Physik der Atmosphare, Vol 10, pp. 91-103 (1922)
+
+    [2] F. Kasten, "A simple parameterization of the pyrheliometric formula for
+    determining the Linke turbidity factor", Meteorologische Rundschau 33,
+    pp. 124-127 (1980)
+
+    [3] Kasten, "The Linke turbidity factor based on improved values of the
+    integral Rayleigh optical thickness", Solar Energy, Vol. 56, No. 3,
+    pp. 239-244 (1996)
+    :doi:`10.1016/0038-092X(95)00114-7`
+
+    [4] B. Molineaux, P. Ineichen, N. O'Neill, "Equivalence of pyrheliometric
+    and monochromatic aerosol optical depths at a single key wavelength",
+    Applied Optics Vol. 37, issue 10, 7008-7018 (1998)
+    :doi:`10.1364/AO.37.007008`
+
+    [5] P. Ineichen, "Conversion function between the Linke turbidity and the
+    atmospheric water vapor and aerosol content", Solar Energy 82,
+    pp. 1095-1097 (2008)
+    :doi:`10.1016/j.solener.2008.04.010`
+
+    [6] P. Ineichen and R. Perez, "A new airmass independent formulation for
+    the Linke Turbidity coefficient", Solar Energy, Vol. 73, no. 3, pp. 151-157
+    (2002)
+    :doi:`10.1016/S0038-092X(02)00045-2`
+    """
+    # "From numerically integrated spectral simulations done with Modtran
+    # (Berk, 1989), Molineaux (1998) obtained for the broadband optical depth
+    # of a clean and dry atmospshere (fictitious atmosphere that comprises only
+    # the effects of Rayleigh scattering and absorption by the atmosphere gases
+    # other than the water vapor) the following expression"
+    # - P. Ineichen (2008)
+    delta_cda = -0.101 + 0.235 * airmass_absolute ** (-0.16)
+    # "and the broadband water vapor optical depth where pwat is the integrated
+    # precipitable water vapor content of the atmosphere expressed in cm and am
+    # the optical air mass. The precision of these fits is better than 1% when
+    # compared with Modtran simulations in the range 1 < am < 5 and
+    # 0 < pwat < 5 cm at sea level" - P. Ineichen (2008)
+    delta_w = 0.112 * airmass_absolute ** (-0.55) * precipitable_water ** 0.34
+    # broadband AOD
+    delta_a = aod_bb
+    # "Then using the Kasten pyrheliometric formula (1980, 1996), the Linke
+    # turbidity at am = 2 can be written. The extension of the Linke turbidity
+    # coefficient to other values of air mass was published by Ineichen and
+    # Perez (2002)" - P. Ineichen (2008)
+    lt = -(9.4 + 0.9 * airmass_absolute) * np.log(
+        np.exp(-airmass_absolute * (delta_cda + delta_w + delta_a))
+    ) / airmass_absolute
+    # filter out of extrapolated values
+    return lt
+
+
+def angstrom_aod_at_lambda(aod0, lambda0, alpha, lambda1=700.0):
+    r"""
+    Get AOD at specified wavelength using Angstrom turbidity model.
+
+    Parameters
+    ----------
+    aod0 : numeric
+        aerosol optical depth (AOD) measured at known wavelength
+    lambda0 : numeric
+        wavelength in nanometers corresponding to ``aod0``
+    alpha : numeric
+        Angstrom :math:`\alpha` exponent corresponding to ``aod0``
+    lambda1 : numeric
+        desired wavelength in nanometers, defaults to 700 nm
+
+    Returns
+    -------
+    aod1 : numeric
+        AOD at desired wavelength, ``lambda1``
+
+    See also
+    --------
+    angstrom_alpha
+
+    References
+    ----------
+    [1] Anders Angstrom, "On the Atmospheric Transmission of Sun Radiation and
+    On Dust in the Air", Geografiska Annaler Vol. 11, pp. 156-166 (1929) JSTOR
+    :doi:`10.2307/519399`
+
+    [2] Anders Angstrom, "Techniques of Determining the Turbidity of the
+    Atmosphere", Tellus 13:2, pp. 214-223 (1961) Taylor & Francis
+    :doi:`10.3402/tellusa.v13i2.9493` and Co-Action Publishing
+    :doi:`10.1111/j.2153-3490.1961.tb00078.x`
+    """
+    return aod0 * ((lambda1 / lambda0) ** (-alpha))
+
+
+def angstrom_alpha(aod1, lambda1, aod2, lambda2):
+    r"""
+    Calculate Angstrom alpha exponent.
+
+    Parameters
+    ----------
+    aod1 : numeric
+        first aerosol optical depth
+    lambda1 : numeric
+        wavelength in nanometers corresponding to ``aod1``
+    aod2 : numeric
+        second aerosol optical depth
+    lambda2 : numeric
+        wavelength in nanometers corresponding to ``aod2``
+
+    Returns
+    -------
+    alpha : numeric
+        Angstrom :math:`\alpha` exponent for AOD in ``(lambda1, lambda2)``
+
+    See also
+    --------
+    angstrom_aod_at_lambda
+    """
+    return - np.log(aod1 / aod2) / np.log(lambda1 / lambda2)


Test patch:
diff --git a/pvlib/test/test_atmosphere.py b/pvlib/test/test_atmosphere.py
index 5522fa3..40e58c7 100644
--- a/pvlib/test/test_atmosphere.py
+++ b/pvlib/test/test_atmosphere.py
@@ -1,4 +1,3 @@
-import datetime
 import itertools
 
 import numpy as np
@@ -7,10 +6,8 @@ import pandas as pd
 import pytest
 from numpy.testing import assert_allclose
 
-from pvlib.location import Location
-from pvlib import solarposition
 from pvlib import atmosphere
-
+from pvlib import solarposition
 
 latitude, longitude, tz, altitude = 32.2, -111, 'US/Arizona', 700
 
@@ -113,3 +110,43 @@ def test_first_solar_spectral_correction_supplied():
 def test_first_solar_spectral_correction_ambiguous():
     with pytest.raises(TypeError):
         atmosphere.first_solar_spectral_correction(1, 1)
+
+
+def test_kasten96_lt():
+    """Test Linke turbidity factor calculated from AOD, Pwat and AM"""
+    amp = np.array([1, 3, 5])
+    pwat = np.array([0, 2.5, 5])
+    aod_bb = np.array([0, 0.1, 1])
+    lt_expected = np.array(
+        [[[1.3802, 2.4102, 11.6802],
+          [1.16303976, 2.37303976, 13.26303976],
+          [1.12101907, 2.51101907, 15.02101907]],
+
+         [[2.95546945, 3.98546945, 13.25546945],
+          [2.17435443, 3.38435443, 14.27435443],
+          [1.99821967, 3.38821967, 15.89821967]],
+
+         [[3.37410769, 4.40410769, 13.67410769],
+          [2.44311797, 3.65311797, 14.54311797],
+          [2.23134152, 3.62134152, 16.13134152]]]
+    )
+    lt = atmosphere.kasten96_lt(*np.meshgrid(amp, pwat, aod_bb))
+    assert np.allclose(lt, lt_expected, 1e-3)
+    return lt
+
+
+def test_angstrom_aod():
+    """Test Angstrom turbidity model functions."""
+    aod550 = 0.15
+    aod1240 = 0.05
+    alpha = atmosphere.angstrom_alpha(aod550, 550, aod1240, 1240)
+    assert np.isclose(alpha, 1.3513924317859232)
+    aod700 = atmosphere.angstrom_aod_at_lambda(aod550, 550, alpha)
+    assert np.isclose(aod700, 0.10828110997681031)
+
+
+def test_bird_hulstrom80_aod_bb():
+    """Test Bird_Hulstrom broadband AOD."""
+    aod380, aod500 = 0.22072480948195175, 0.1614279181106312
+    bird_hulstrom = atmosphere.bird_hulstrom80_aod_bb(aod380, aod500)
+    assert np.isclose(0.11738229553812768, bird_hulstrom)