Linke Turbidity factor uses nearest neighbor instead of bilinear interpolation

[mikofski]

lookup_linke_turbidity() uses "nearest neighbor" and not bilinear interpolation, so it is an approximation. The maximum error is over 530[km] which seems large to me:

I'll send a patch for your consideration.

[wholmgren]

The data is pretty rough to begin with, but go for it. I only added the time interpolation because I didn't like the month-to-month discontinuities in my modeling results, not because I thought it would make them any more accurate.

That maximum error does seem large. Maybe we'd need to think harder about the map projection.

[mikofski]

I am working on this, but I've still some bugs - need to send pull request, but here is patch:

"""
Bilinear interpolation for Linke turbidity factor.
"""

import scipy.io as scio
import numpy as np
import os, sys
from pvlib.clearsky import lookup_linke_turbidity
import pandas as pd
from functools import partial
import logging

logging.basicConfig()
LOGGER = logging.getLogger(__name__)
LOGGER.setLevel(logging.DEBUG)

LTSHAPE = (2160, 4320)  # get this from mat file?
DX, DY = 180.0 / LTSHAPE[0], 360.0 / LTSHAPE[1]
LTFILE = os.path.join('data', 'LinkeTurbidities.mat')


class LinkeConfig:
    ltfile = LTFILE
    data = scio.loadmat(ltfile)
    linke_turbidity_table = data['LinkeTurbidity'] / 20.0


def wikipedia_interp2d(x, y, f, ix, iy, dx=DX, dy=DY):
    """
    Bilinear interpolation for (x, y) bounded in the rectangle between nodes
    (x1, y1) and (x1 + dx, y1 + dx).

    :param x: desired x coordinate
    :param y: desired y coordinate
    :param f: the output field
    :param ix: lower bounding nodes in x-direction
    :param iy: lower bounding nodes in y-direction
    :param dx: spacing between x nodes, default = 1
    :param dy: spacing between y nodes, default = 1
    https://en.wikipedia.org/wiki/Bilinear_interpolation
    """
    x1, y1 = float(ix), float(iy)
    x2, y2 = x1 + dx, y1 + dy  # upper bounding nodes
    return np.dot(np.dot([x2 - x, x - x1], f).T, [[y2 - y], [y - y1]]) / dx / dy


def get_linke_turbidity_at_latlon(latitude, longitude):
    """
    Get the Linke Turbidity from the ``LinkeTurbidities.mat`` data file.
    :param latitude: latitude [degrees]
    :type latitude: float
    :param longitude: longitude
    :type longitude: float
    :return: Linke turbdity factor
    :rtype: np.array
    """
    lat, lon = latitude, longitude
    # int() always rounds towards zero
    ilat_0 = int((90.0 - lat) / DX)  # lower bound
    # avoid out of bounds indices
    if ilat_0 < 0 or ilat_0 > LTSHAPE[0]:
        raise IndexError('latitude should be between -90 and 90 degrees')
    if ilat_0 == LTSHAPE[0]:
        ilat_0 = LTSHAPE[0] - 1  # need room to interpolate
    ilat_1 = ilat_0 + 1  # upper bound
    lon = lon % 360.0  # adjust longitude range from -180 - 180 to 0 - 360
    ilon_0 = int(lon / DY) % LTSHAPE[1]  # lower bound
    ilon_1 = (ilon_0 + 1) % LTSHAPE[1]  # upper bound
    f11 = LinkeConfig.linke_turbidity_table[ilat_0][ilon_0]
    f12 = LinkeConfig.linke_turbidity_table[ilat_0][ilon_1]
    f21 = LinkeConfig.linke_turbidity_table[ilat_1][ilon_0]
    f22 = LinkeConfig.linke_turbidity_table[ilat_1][ilon_1]
    f = np.asarray([[f11, f12], [f21, f22]])
    return wikipedia_interp2d(latitude, longitude, f, 90.0 - ilat_0 * DX, ilon_0 * DY)


TIME = pd.DatetimeIndex(['%d/1/2016' % m for m in xrange(1,13)])


def test():
    lookup_linke_turbidity_test = partial(
        lookup_linke_turbidity, filepath=LTFILE, interp_turbidity=False
    )
    test1 = [79, 179]
    lt_test1 = get_linke_turbidity_at_latlon(*test1)
    lt_bench1 = lookup_linke_turbidity_test(time, *test1)
    LOGGER.debug('test1:\n%r', lt_test1)
    LOGGER.debug('bench1:\n%r', lt_bench1)
    assert lt_test == lookup_linke_turbidity_test(time, *test1)


if __name__ == '__main__':
    print get_linke_turbidity_at_latlon(*[float(ll) for ll in sys.argv[1:]])