[Proj] Robinson projection coefficients

[Ed Campbell]

A trac ticket on the Robinson projection has recently been reopened
(http://trac.osgeo.org/proj/ticket/113). The issue is that an applied
correction to the constants that define the Robinson projection was
only a partial fix. The higher-order interpolation coefficients also
needed to be updated. Below is a short python script that I have
written to determine the coefficients:

#!/usr/bin/env python2.7
"""
Short script to calculate the interpolation coefficients
used by proj.4 Robinson projection.

.. note::

    The scipy interpolate method is not a natural cubic
    spline (curvature zero at the endpoints).

"""
import sys

import numpy
import scipy.interpolate


def spline_coeffs(x, y):
    """
    Calculates the 5 deg interpolation coefficients based on
    input arrays x and y where y = f(x).

    """
    s = scipy.interpolate.InterpolatedUnivariateSpline(x, y)
    coeffs = []
    for t in  numpy.linspace(0.0, 90.0, 19):
        c0 = s(t)
        c1 = s(t, 1)
        c2 = s(t, 2) / 2.0
        c3 = s(t, 3) / 6.0
        coeffs.append([numpy.float64(c) for c in (c0, c1, c2, c3)])
    return coeffs


def print_coeffs(coeffs, var_name, file_out=None):
    """Output the coefficients formatted as a C code."""
    if file_out is None:
       file_out = sys.stdout
    file_out.write('{}[] = {{\n'.format(var_name))
    for i, c in enumerate(coeffs):
        if i != 0:
            file_out.write(',\n')
        file_out.write('    {{{:.6g}, {:.6g}, {:.6g}, {:.6g}}}'.format(*c))
    file_out.write(' };\n')


def main():
    # 0 to 90 deg values
    x = numpy.array([1.0, 0.9986, 0.9954, 0.99, 0.9822, 0.973, 0.96,
         0.9427, 0.9216, 0.8962, 0.8679, 0.835, 0.7986,
         0.7597, 0.7186, 0.6732, 0.6213, 0.5722, 0.5322])

    y = numpy.array([0.0, 0.062, 0.124, 0.186, 0.248, 0.31, 0.372,
         0.434, 0.4958, 0.5571,0.6176, 0.6769, 0.7346,
         0.7903, 0.8435, 0.8936, 0.9394, 0.9761, 1.0])

    # Produce -90 to  90 arrays
    phi = numpy.linspace(-90.0, 90.0, 37)
    x = numpy.append(x[::-1], x[1:])
    y = numpy.append(-y[::-1], y[1:])

    # Input values
    print('Input values:')
    print('phi = {}'.format(phi))
    print('x = {}'.format(x))
    print('y = {}'.format(y))

    print('=====================================')
    xcoeffs = spline_coeffs(phi, x)
    print_coeffs(xcoeffs, 'X')

    ycoeffs = spline_coeffs(phi, y)
    print_coeffs(ycoeffs, 'Y')
    print('=====================================')

    return xcoeffs, ycoeffs

if __name__ == '__main__':
    main()

This code produces the following:

X[] = {
    {1, 2.2199e-17, -7.15515e-05, 3.1103e-06},
    {0.9986, -0.000482243, -2.4897e-05, -1.3309e-06},
    {0.9954, -0.00083103, -4.48605e-05, -9.86701e-07},
    {0.99, -0.00135364, -5.9661e-05, 3.6777e-06},
    {0.9822, -0.00167442, -4.49547e-06, -5.72411e-06},
    {0.973, -0.00214868, -9.03571e-05, 1.8736e-08},
    {0.96, -0.00305085, -9.00761e-05, 1.64917e-06},
    {0.9427, -0.00382792, -6.53386e-05, -2.6154e-06},
    {0.9216, -0.00467746, -0.00010457, 4.81243e-06},
    {0.8962, -0.00536223, -3.23831e-05, -5.43432e-06},
    {0.8679, -0.00609363, -0.000113898, 3.32484e-06},
    {0.835, -0.00698325, -6.40253e-05, 9.34959e-07},
    {0.7986, -0.00755338, -5.00009e-05, 9.35324e-07},
    {0.7597, -0.00798324, -3.5971e-05, -2.27626e-06},
    {0.7186, -0.00851367, -7.01149e-05, -8.6303e-06},
    {0.6732, -0.00986209, -0.000199569, 1.91974e-05},
    {0.6213, -0.010418, 8.83923e-05, 6.24051e-06},
    {0.5722, -0.00906601, 0.000182, 6.24051e-06},
    {0.5322, -0.00677797, 0.000275608, 6.24051e-06} };
Y[] = {
    {-5.20417e-18, 0.0124, 1.21431e-18, -8.45284e-11},
    {0.062, 0.0124, -1.26793e-09, 4.22642e-10},
    {0.124, 0.0124, 5.07171e-09, -1.60604e-09},
    {0.186, 0.0123999, -1.90189e-08, 6.00152e-09},
    {0.248, 0.0124002, 7.10039e-08, -2.24e-08},
    {0.31, 0.0123992, -2.64997e-07, 8.35986e-08},
    {0.372, 0.0124029, 9.88983e-07, -3.11994e-07},
    {0.434, 0.0123893, -3.69093e-06, -4.35621e-07},
    {0.4958, 0.0123198, -1.02252e-05, -3.45523e-07},
    {0.5571, 0.0121916, -1.54081e-05, -5.82288e-07},
    {0.6176, 0.0119938, -2.41424e-05, -5.25327e-07},
    {0.6769, 0.011713, -3.20223e-05, -5.16405e-07},
    {0.7346, 0.0113541, -3.97684e-05, -6.09052e-07},
    {0.7903, 0.0109107, -4.89042e-05, -1.04739e-06},
    {0.8435, 0.0103431, -6.4615e-05, -1.40374e-09},
    {0.8936, 0.00969686, -6.4636e-05, -8.547e-06},
    {0.9394, 0.00840947, -0.000192841, -4.2106e-06},
    {0.9761, 0.00616527, -0.000256, -4.2106e-06},
    {1, 0.00328947, -0.000319159, -4.2106e-06} };

These can be compared directly with the arrays in PJ_robin.c. Note
that the values at the extremes differ slightly, possibly because of a
difference in the interpolation algorithm and/or the boundary
conditions. I'd appreciate it if someone could post details on the
exact method used to calculate the coefficients in the proj.4 source
code.


Ed Campbell PhD
UK Met Office
Iris & Cartopy - http://scitools.org.uk/