Skip to content

add linke turbidity formulas module to atmosphere #278

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 22 commits into from
Feb 6, 2017

Conversation

mikofski
Copy link
Member

@mikofski mikofski commented Dec 3, 2016

Signed-off-by: Mark Mikofski [email protected]

* 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 <[email protected]>
@wholmgren wholmgren added this to the 0.4.4 milestone Dec 28, 2016
* append linke_turb_forms.py to atmosphere.py
* rm -rf atmosphere/

Signed-off-by: Mark Mikofski <[email protected]>
* import xrange from six for python 3 compatibility
* fix pyrheliometric formula is from 1980
* just use `range` instead of importing `six`
* change function to `kasten80_lt` to match `gueymard94_pw`
* start adding references
* add lots of references
* cite text quoted directly from Ineichen
* acknowledge Armel
* fix docstrings for numpydocs
@mikofski
Copy link
Member Author

I think this is ready. The Travis failure is not related to kasten96_lt()

Copy link
Member

@wholmgren wholmgren left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @mikofski.

Is there any way to simplify the API? It feels very complicated to me. Also, the various types of potential inputs are not tested.

Please also add the functions to the api.rst file and update whatsnew.


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``.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not clear to me that we need to set these limits. Garbage in, garbage out.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

okay, removed filters and put these limits as a .. warning::

than 1 or greater than 6 will return ``NaN``. Precipitable water less than
zero or greater than 5[cm] will also return ``NaN``.

Based on originam implementation by Armel Oumbe.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

original

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Linke turbidity

References
----------
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there are a handful of typos in the references

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

1-3 good references might be more useful than 8.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I fixed the references errors. Not sure I agree about the references. It's hard for me to tease out all of these pieces from just one or two references. I suppose we could just use the Ineichen references which include almost all of the other ones, although some are incorrect:

  • "A new airmass independent ..." (2002) includes: Kasten-1980, Kasten-1996, Linke-1922, Molineaux-1989
  • "Conversion function between Linke ..." (2008) includes: Berk-1989,1996, Bird-1980 (incorrect reference), Mueller and Mayer

What is your preference?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm ok with a longer references section if that's your preference.

if method == 'Molineaux':
# get aod at 700[nm] from alpha for Molineaux (1998)
delta_a = angstrom_aod_at_lambda(aod, alpha)
else:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

probably best to elif method == 'bird-hulstrom' and then else: raise ValueError('invalid method')

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Angstrom :math:`\\alpha` exponent corresponding to wavelength in ``aod``
lambda0 : numeric
desired wavelength in nanometers, defaults to 700[nm]

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need a Returns section

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

return lt


def angstrom_aod_at_lambda(aod, alpha, lambda0=700.0):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this function needs its own tests

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

atmos_path = os.path.dirname(os.path.abspath(__file__))
pvlib_path = os.path.dirname(atmos_path)
melbourne_fl = tmy.readtmy3(os.path.join(
pvlib_path, 'data', '722040TYA.CSV')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

why add a new tmy file instead of using the existing sample file?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no particular reason, I like that Melbourne-FL site, it's close to Cocoa where I have irradiance data. also although no site matches the Linke turbidity data well, the Alaska site has bigger differences, and also bigger differences between broadband approximations relative to Melbourne-FL. Maybe these are not good tests?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'd prefer to define ~2-5 element arrays that span the range of allowed input values, and also define the expected return values. I can't quickly inspect the values that are being tested and expected in this test. Maybe these values span the space well enough?

am = np.array([1, 3, 5])
pw = np.array([0, 2.5, 5])
aod700 = np.array([0, 0.1, 1])
lt_molineaux_expected = np.array('your precalculated values')
lt_bird_hulstrom_expected = np.array('your precalculated values')
...
np.assert_allclose(lt_molineaux_expected, lt_molineaux)
np.assert_allclose(lt_bird_hulstrom_expected, lt_bird_hulstrom)

Then you don't need to read a tmy file or lookup linke turbidity in the test. Seems more simple and more explicit to me. I'm not sure how much I trust the tmy and linke turbidity files anyways.

@mikofski
Copy link
Member Author

@wholmgren great feedback on all points, thanks! 👍👍

* remove alpha argument
* update docstring to explain that if molineaux method, then aod is
shape is only aod at 700[nm] but if bird-hulstrom, then aod should be
both aod380 and aod500
* add see also to angstrom_aod_at_lambda
* add warning for air mass or pwat out of range
* check for bird-hulstrom, raise type error if neither molineaux nor
* also simplify angstrom_aod_at_lambda() api, change aod to aod0,
change lambda0 to lambda1, and add lambda0, update docstring
* update test_atmosphere.test_kasten96_lt with new changes
* separate arguments so their shapes don't change depending on method,
ie: before aod could be 1, 2, N, N x 1 or N x 2 depending on whether
Molineaux or Bird-Hulstrom were used, now, there are 3 separate args,
 aod700, aod380 and aod500.
* fix misspelling of Hulstrom, no "d", note several sources also make
 this error :(
* add new methods to api.rst
* add test_angstrom tests
* add raises value error test for bad method
* add lower() to compare only lower case of method name
* fix more hulstrom errors, no "d"
* fix typos in referecnes
@mikofski
Copy link
Member Author

ready for your review:

  • simplified kasten96_lt API
  • added tests for angstrom methods
  • added test for bad method raises ValueError
  • add new methods to api.rst
  • update what's new 0.4.4

@wholmgren
Copy link
Member

wholmgren commented Jan 31, 2017 via email

@wholmgren
Copy link
Member

These are great changes, thanks! landscape is complaining about a few things related to the doc strings.

@mikofski
Copy link
Member Author

mikofski commented Feb 1, 2017

Two of the landscape docstring warnings:

Use r""" if any backslashes in a docstring

are due to math, but I just escaped the backslashes, so they are working fine as is, no need to use a "raw" string
angstrom_docstring

@mikofski
Copy link
Member Author

mikofski commented Feb 1, 2017

The other landscape errors are due to the DOI external links which are too long. Can we add doi to extlink in conf.py like you did for issues and wiki?

extlinks = {
    'doi': ('https://dx.doi.org/%s', 'doi:'),
}

@wholmgren
Copy link
Member

Sorry, didn't look carefully enough. All sounds good.

* add 'doi': ('http://dx/doi.org/%s', 'DOI: ') to extlink in conf.py
* change doi links everywhere to use :doi:`doi number` instead of links
* instead of """docstring with :math:`\\alpha`
* add Co-Action publishing for original Tellus A DOI before moving to
 Taylor & Francis
@mikofski
Copy link
Member Author

mikofski commented Feb 2, 2017

Okay, I think it's ready now. I added doi as an extlink to the conf.py and I change the docstring to use r"""docstring with :math:\alpha""" to make landscape happy. Those are both pretty neat features, so something new to learn. 😺

@mikofski
Copy link
Member Author

mikofski commented Feb 2, 2017

Okay, that was my last change. Sorry! 😉 I really think it's ready now. Let me know if there are any more comments or suggestions.

@mikofski
Copy link
Member Author

mikofski commented Feb 2, 2017

For the remaining landscape complaint:

Too many arguments (6/5)

That's the simplified API for the Kasten LT formula. I had some ideas about that but I wasn't what was best. Basically one idea would be to just have one AOD input called aod_bb for broadband AOD. And then provide one more short formulas for Bird-Hulstrom and remove all of the methods checking. That would also remove 2 references.

def bird_hulstrom80_aod_bb(aod380, aod500):
    """Calculate broadband AOD using Bird-Hulstrom approximation"""
    # put Bird-Hulstrom references here
    return 0.27583 * aod380 + 0.35 * aod500


def kasten96_lt(am, pwat, aod_bb):
    """
    Calculate Linke turbidity factor using Kasten pyrheliometric formula.

    .. warning::
        These calculations are only valid for air mass less than 5[atm] and
        precipitable water less than 5[cm].

    Parameters
    ----------
    am : numeric
        airmass, pressure corrected in atmospheres
    pwat : numeric
        precipitable water or total column water vapor in centimeters
    aod_bb : numeric
        broadband AOD or AOD measured at 700[nm] 

    Returns
    -------
        Linke turbidity

    See also
    --------
    bird_hulstrom80_aod_bb
    angstrom_aod_at_lambda
    """
    # remove Bird-Hulstrom references
    delta_cda = -0.101 + 0.235 * am ** (-0.16)  # rayleigh "clean dry atmosphere"
    delta_w = 0.112 * am ** (-0.55) * pwat ** (0.34)  # water vapor
    delta_a = aod_bb
    return -(9.4 + 0.9 * am) * np.log(np.exp(-am * (delta_cda + delta_w + delta_a))) / am

Would you prefer this approach? Let me know and I'll make the change :)

@wholmgren
Copy link
Member

Nice. I really like the simplicity of your latest proposal. I'm guessing that most users would find it easier to use.

How about this for alternative aod_bb docstring: "Broadband AOD. Note that AOD measured at 700 nm is a good approximation." Up to you. Also, we don't usually use [] around units in pvlib python. I think text reads better with only a space and that seems to be the standard in all of the journals I read (figure labels are a different story).

For consistency with existing functions, we'd want to use airmass_absolute and precipitable_water.

What did you think about simplifying the tests by using a couple of hardcoded values instead of a TMY file?

It's definitely not necessary for this PR, but it would be nice to have a brief section in clearsky.rst (or maybe a new atmosphere.rst) with a few examples and a little bit of discussion. Something for the future.

@mikofski
Copy link
Member Author

mikofski commented Feb 2, 2017

Sounds good. I'll make the changes and get back to you.

* change kasten96_lt AOD args to just a single arg called aod_bb for
broadband AOD
* remove "method" arg and if-then-else on method.lower() in
kasten96_lt()
* add note that according to Molineaux, aod at 700 nm is a good
approximation for broadband aod
* add new function bird_hulstrom80_aod_bb to atmosphere.py
* add it to api.rst
* add new test for bird_hulstrom80_aod_bb() to test_atmosphere.py
* update test_kasten96_lt() to use aod_bb
* use verbose names:
  - s/am/airmass_absolute/g
  - s/pwat/precipitable_water/g
* fix angstrom aod test to use more realistic values
* don't use brackets around units
* remove test for raises ValueError in test_kasten
@mikofski
Copy link
Member Author

mikofski commented Feb 3, 2017

almost there, just cleaning up the data for the test, and removing the TMY file, should be ready soon.

* values span pwat range 0 - 5 cm, and airmass between 1 & 2.5 atm
* remove tmy3 file
* don't try Bird-Hulstrom, not necessary
* make pylint happy, isclose(a, b) a & b should be array or iterable
@mikofski
Copy link
Member Author

mikofski commented Feb 4, 2017

It's ready. 🚀

I made a separate issue for the docs. 📚 I feel like it may require some iterations, so would rather split into a new PR, OK?

Copy link
Member

@wholmgren wholmgren left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I will merge as soon as the asserts are fixed

aod550 = 0.15
aod1240 = 0.05
alpha = atmosphere.angstrom_alpha(aod550, 550, aod1240, 1240)
np.isclose(alpha, [1.3513924317859232])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to assert np.allclose or use np.testing.assert_allclose

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Weird. np.isclose should work since these are scalars. Pylint told me to change the scalar to an intake and it passes tests fine

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

*iterable

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

changed these back to scalars, although pylint now complains, oh well.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see, it's pylint that's the trouble:

In [4]: np.isclose(2, [2])
Out[4]: array([ True], dtype=bool)

In [5]: np.isclose(2, 2)
Out[5]: True

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's that you have to assert something is true or not. np.isclose only returns an array.

alpha = atmosphere.angstrom_alpha(aod550, 550, aod1240, 1240)
np.isclose(alpha, [1.3513924317859232])
aod700 = atmosphere.angstrom_aod_at_lambda(aod550, 550, alpha)
np.isclose(aod700, [0.10828110997681031])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same assert

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ditto

"""Test Bird_Hulstrom broadband AOD."""
aod380, aod500 = 0.22072480948195175, 0.1614279181106312
bird_hulstrom = atmosphere.bird_hulstrom80_aod_bb(aod380, aod500)
np.isclose([0.09823143641608373], bird_hulstrom)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same assert

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ditto

def test_kasten96_lt():
"""Test Linke turbidity factor calculated from AOD, Pwat and AM"""
timestamps = pd.DatetimeIndex(MELBOURNE_FL[0])
timestamps = timestamps.tz_localize('UTC').tz_convert('Etc/GMT+5')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

one more that I just noticed... this line errors on the python27-min configuration because timestamps is already tz aware: https://travis-ci.org/pvlib/pvlib-python/jobs/198282511

I don't understand why the python-36 build completed, but we should still fix it.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was not tz aware on my machine, it gets converted to a naive date time index in utc, that's why I had to localize as UTC then convert to us/Eastern. Is this a pandas or pytz version issue?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is what I get with pandas

In [21]: pd.DatetimeIndex(MELBOURNE_FL[0])
Out[21]:
DatetimeIndex(['1999-01-31 17:00:00', '2000-02-20 20:00:00',
               '2000-02-22 18:00:00', '2000-02-24 20:00:00',
               '1995-03-02 19:00:00', '1995-03-11 17:00:00',
               '1995-03-12 18:00:00', '1995-03-20 16:00:00',
               '1995-03-20 19:00:00', '1995-03-22 16:00:00',
               '1995-03-22 19:00:00', '1995-04-07 14:00:00',
               '1995-04-10 14:00:00', '1995-04-21 14:00:00',
               '2004-05-01 13:00:00', '2004-05-03 13:00:00',
               '2004-05-05 18:00:00', '2004-05-16 14:00:00',
               '2004-05-21 20:00:00', '2004-05-24 16:00:00',
               '2004-05-31 14:00:00', '2002-06-04 21:00:00',
               '2002-06-16 22:00:00', '2002-06-17 13:00:00',
               '2000-07-01 13:00:00', '2000-07-05 12:00:00',
               '2000-07-06 12:00:00', '2000-07-06 22:00:00',
               '2000-07-25 20:00:00', '2001-10-14 18:00:00',
               '2001-10-15 16:00:00', '2003-11-04 19:00:00',
               '1999-12-20 18:00:00'],
              dtype='datetime64[ns]', freq=None)

In [22]: pd.__version__
Out[22]: u'0.19.1'

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the conda list command in the travis builds reports pytz 2016.10 for both builds, so maybe it's pandas. I couldn't find anything about a change in the pandas behavior here, so maybe there's a bug in the newer version. strange.

Copy link
Member Author

@mikofski mikofski Feb 4, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no change in newest pandas or pytz

In [3]: pd.__version__
Out[3]: u'0.19.2'

In [4]: pd.DatetimeIndex(MELBOURNE_FL[0])
Out[4]:
DatetimeIndex(['1999-01-31 17:00:00', '2000-02-20 20:00:00',
               '2000-02-22 18:00:00', '2000-02-24 20:00:00',
               '1995-03-02 19:00:00', '1995-03-11 17:00:00',
               '1995-03-12 18:00:00', '1995-03-20 16:00:00',
               '1995-03-20 19:00:00', '1995-03-22 16:00:00',
               '1995-03-22 19:00:00', '1995-04-07 14:00:00',
               '1995-04-10 14:00:00', '1995-04-21 14:00:00',
               '2004-05-01 13:00:00', '2004-05-03 13:00:00',
               '2004-05-05 18:00:00', '2004-05-16 14:00:00',
               '2004-05-21 20:00:00', '2004-05-24 16:00:00',
               '2004-05-31 14:00:00', '2002-06-04 21:00:00',
               '2002-06-16 22:00:00', '2002-06-17 13:00:00',
               '2000-07-01 13:00:00', '2000-07-05 12:00:00',
               '2000-07-06 12:00:00', '2000-07-06 22:00:00',
               '2000-07-25 20:00:00', '2001-10-14 18:00:00',
               '2001-10-15 16:00:00', '2003-11-04 19:00:00',
               '1999-12-20 18:00:00'],
              dtype='datetime64[ns]', freq=None)

In [6]: pytz.__version__
Out[6]: '2016.10'

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I get the same thing on pandas 0.19.2. But also this...

In [11]: pd.Timestamp(MELBOURNE_FL[0][0])
Out[11]: Timestamp('1999-01-31 12:00:00-0500', tz='pytz.FixedOffset(-300)')

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@mikofski
Copy link
Member Author

mikofski commented Feb 4, 2017

Anyway, I didn't like this test and want to change it so this gives me the chance.

* fix assert np.isclose() should be scalar only
* super simplify test_kasten96_lt() to ranges of inputs and expected
outputs, don't need to compare to actual data, just test that method
is working nominally
@mikofski
Copy link
Member Author

mikofski commented Feb 4, 2017

let's see how this goes. sorry for all the iterations. ☹️



def test_angstrom_aod():
"""Test Angstrom turbidity model functions."""
aod550 = 0.15
aod1240 = 0.05
alpha = atmosphere.angstrom_alpha(aod550, 550, aod1240, 1240)
np.isclose(alpha, [1.3513924317859232])
np.isclose(alpha, 1.3513924317859232)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sorry for not being clear. this line and others like it should read assert np.isclose(alpha, 1.3513924317859232)

@wholmgren
Copy link
Member

no worries. it would help if the "so called judge" (the test suite) actually worked reliably.

@mikofski
Copy link
Member Author

mikofski commented Feb 4, 2017

Landscape to the rescue - cleaned up unused import (clearsky). Travis reports at least one Python 2 and one Python 3 tests each pass, the others stall or error out on unrelated issues.

@wholmgren
Copy link
Member

@mikofski In case my comment was buried in the exchange on Saturday morning... the functions that use np.isclose need assert statements in order for them to actually test something: #278 (comment)

I'll also close #101 when this is merged.

@mikofski
Copy link
Member Author

mikofski commented Feb 6, 2017

Oops! Sorry about that! Hopefully good to go.

@wholmgren
Copy link
Member

Great, thanks Mark!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

calculate Linke turbidity factor from AOD, precip water and airmass
2 participants