Calculate total alkalinity according to Lee et al. (2006)

alkalinity_Lee2006.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
This script is written to calculate total alkalinity according to Lee et al. (2006).
Please see function help for more details. 
Please acknowledge use.
"""

__author__ = "Luke Gregor"
__email__ = "luke.e.gregor@gmail.com"
__date__ = "2017-05-19"

import numpy as np


def calc_TA_lee2006(lat, lon, salt, tempC, return_details=False):
    """
    Calculates total alkalinity according to:
    Lee, K., Tong, L. T., Millero, F. J., Sabine, C. L., Dickson, A. G.,
        Goyet, C., … Key, R. M. (2006). Global relationships of total alkalinity
        with salinity and temperature in surface waters of the world’s oceans.
        Geophysical Research Letters, 33(19), L19605.
        https://doi.org/10.1029/2006GL027207

    Accepts a single location or a numpy array of locations with associated measurements.
    lat   =  latitude
    lon   =  longitude
    salt  =  salinity
    tempC =  temperature in degrees C
    return_details = True (default), False

    if return_details is set to True a pandas.DataFrame will be returned
        with original input and alkalinity, zone number and name, and standard deviation

    Please acknowledge usage of this script.
    """

    def is_eqPac(yA, xA):
        # north south / east west boundaries
        p0 = (abs(yA) < 20) & (xA > -140) & (xA < -75)

        # southern diagonal boundary
        x1, y1 = -140, -10
        x2, y2 = -110, -20
        # create two vectors and then calculate the cross product
        v1 = np.array([x2 - x1, y2 - y1])
        v2 = np.array([x2 - xA, y2 - yA])
        # if p1 is negative it is in EQ Pac
        p1 = (v1[0] * v2[1] - v1[1] * v2[0]) < 0

        # northern diagonal boundary
        x1, y1 = -140, 10
        x2, y2 = -110, 20
        # create two vectors and then calculate the cross product
        v1 = np.array([x2 - x1, y2 - y1])
        v2 = np.array([x2 - xA, y2 - yA])
        # if p2 is positive it is in EQ Pac
        p2 = (v1[0] * v2[1] - v1[1] * v2[0]) > 0

        return (p0 * p1 * p2)

    y = lat
    x = lon
    s = salt
    t = tempC

    zones = {
        1: {
            "name": "(Sub)tropics",
            "condition": lambda y, x, s, t: (y <= 30.) & (y >= -30) & (t > 20.) & (s > 31.) & (s < 38.) & (~is_eqPac(y, x)),
            "function": lambda y, x, s, t: (2305. + 58.66 * (s - 35.) + 2.32 * (s - 35.)**2 - 1.41 * (t - 20.) + 0.040 * (t - 20.)**2),
            "std_reported": 8.6,
        },
        2: {
            "name": "Pacific equatorial",
            "condition": lambda y, x, s, t: is_eqPac(y, x) & (s > 31) & (s < 36.5) & (t > 18),
            "function": lambda y, x, s, t: (2294. + 64.88 * (s - 35.) + 0.39 * (s - 35.)**2 - 4.52 * (t - 29.) + 0.232 * (t - 29)**2),
            "std_reported": 7.5,
        },
        3: {
            "name": "North Atlantic",
            "condition": lambda y, x, s, t: (y > 30) & (x > -90) & (x < 75) & (t > 0) & (t < 20) & (s > 31) & (s < 37),
            "function": lambda y, x, s, t: (2305. + 53.97 * (s - 35.) + 2.74 * (s - 35.)**2 - 1.16 * (t - 20.) + 0.040 * (t - 20)**2),
            "std_reported": 6.4,
        },
        4: {
            "name": "North Pacific",
            "condition": lambda y, x, s, t: (y > 30) & (y < -80) & (x > 120) & (x < -105) & (t > 0) & (t < 20) & (s > 31) & (s < 35),
            "function": lambda y, x, s, t: (2305. + 53.97 * (s - 35.) + 2.74 * (s - 35.)**2 - 1.16 * (t - 20.) + 0.040 * (t - 20)**2),
            "std_reported": 8.7,
        },
        5: {
            "name": "Southern Ocean",
            "condition": lambda y, x, s, t: (y < -30.) & (y > -70) & (t < 20.0) & (s > 33.0) & (s < 36.),
            "function": lambda y, x, s, t: (2305. + 52.48 * (s - 35.) + 2.85 * (s - 35.)**2 - 0.49 * (t - 20.) + 0.086 * (t - 20)**2),
            "std_reported": 8.4,
        }
    }

    alkalinity = (
        (zones[1]["function"](y, x, s, t) * zones[1]["condition"](y, x, s, t)) +
        (zones[2]["function"](y, x, s, t) * zones[2]["condition"](y, x, s, t)) +
        (zones[3]["function"](y, x, s, t) * zones[3]["condition"](y, x, s, t)) +
        (zones[4]["function"](y, x, s, t) * zones[4]["condition"](y, x, s, t)) +
        (zones[5]["function"](y, x, s, t) * zones[5]["condition"](y, x, s, t))
    )

    if not return_details:
        return alkalinity

    iszone = (
        (1 * zones[1]["condition"](y, x, s, t)) +
        (2 * zones[2]["condition"](y, x, s, t)) +
        (3 * zones[3]["condition"](y, x, s, t)) +
        (4 * zones[4]["condition"](y, x, s, t)) +
        (5 * zones[5]["condition"](y, x, s, t))
    )

    dct = {
        "latitude": lat,
        "longitude": lon,
        "salinity": salt,
        "temperatureC": tempC,
        "alkalinity": alkalinity,
        "zone_number": iszone,
        "zone_name": [zones[z]["name"] for z in iszone],
        "std_deviation": [zones[z]["std_reported"] for z in iszone],
    }
    try:
        from pandas import DataFrame
        return DataFrame.from_dict(dct)
    except ImportError:
        return dct


def test():
      
    lat = np.array([0,  10, 19])
    lon = np.array([-110, -130, -120])
    sss = np.array([34.5, 34.5, 35.5])
    sst = np.array([20.,  28., 25])

    print (calc_TA_lee2006(lat, lon, sss, sst, return_details=True))


if __name__ == '__main__':
    test()