# [Proj] Robinson projection coefficients

Ed Campbell drescampbell at googlemail.com
Thu May 2 13:37:35 EST 2013

```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/
```