From 1391ad05bed6ed6217a1a1563858d0ad3c04fd69 Mon Sep 17 00:00:00 2001 From: Eric Firing Date: Sun, 2 Oct 2011 11:48:39 -1000 Subject: [PATCH 1/2] geo.py: prevent division by zero in Mollweide transform The Newton-Raphson iteration converges more slowly near the poles, and triggers a divide-by-zero at the poles. Although the previous transform was handling this correctly in the sense that it was still giving the right answer, the divide-by-zero warning was non-optimal. The present version switches to a Taylor-series approximation for latitudes within 5 degrees of the poles. This avoids the divide-by-zero, and also runs faster on arrays with a full range of latitudes. --- lib/matplotlib/projections/geo.py | 36 ++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/lib/matplotlib/projections/geo.py b/lib/matplotlib/projections/geo.py index 8a20f5c2768d..182266ca0a62 100644 --- a/lib/matplotlib/projections/geo.py +++ b/lib/matplotlib/projections/geo.py @@ -436,23 +436,35 @@ def __init__(self, resolution): def transform(self, ll): def d(theta): delta = -(theta + np.sin(theta) - pi_sin_l) / (1 + np.cos(theta)) - return delta, abs(delta) > 0.001 + return delta, np.abs(delta) > 0.001 - longitude = ll[:, 0:1] - latitude = ll[:, 1:2] + longitude = ll[:, 0] + latitude = ll[:, 1] + + clat = np.pi/2 - np.abs(latitude) + ihigh = clat < 0.087 # within 5 degrees of the poles + ilow = ~ihigh + aux = np.empty(latitude.shape, dtype=np.float) - pi_sin_l = np.pi * np.sin(latitude) - theta = 2.0 * latitude - delta, large_delta = d(theta) - while np.any(large_delta): - theta += np.where(large_delta, delta, 0) + if ilow.any(): # Newton-Raphson iteration + pi_sin_l = np.pi * np.sin(latitude[ilow]) + theta = 2.0 * latitude[ilow] delta, large_delta = d(theta) - aux = theta / 2 + while np.any(large_delta): + theta[large_delta] += delta[large_delta] + delta, large_delta = d(theta) + aux[ilow] = theta / 2 - x = (2.0 * np.sqrt(2.0) * longitude * np.cos(aux)) / np.pi - y = (np.sqrt(2.0) * np.sin(aux)) + if ihigh.any(): # Taylor series-based approx. solution + e = clat[ihigh] + d = 0.5 * (3 * np.pi * e**2) ** (1.0/3) + aux[ihigh] = (np.pi/2 - d) * np.sign(latitude[ihigh]) - return np.concatenate((x, y), 1) + xy = np.empty(ll.shape, dtype=np.float) + xy[:,0] = (2.0 * np.sqrt(2.0) / np.pi) * longitude * np.cos(aux) + xy[:,1] = np.sqrt(2.0) * np.sin(aux) + + return xy transform.__doc__ = Transform.transform.__doc__ transform_non_affine = transform From 533e5fba0bee341bec31bb2987c19e86b05cb631 Mon Sep 17 00:00:00 2001 From: Eric Firing Date: Sun, 2 Oct 2011 12:55:28 -1000 Subject: [PATCH 2/2] geo.py: fix divide-by-zero warning in Aitoff transform An attempt to do this had been make by masking zero values in the denominator; but at least with current numpy, this is not enough, because the __div__ method of the first argument is used. The solution is to ensure that the first argument is also a masked array. --- lib/matplotlib/projections/geo.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/matplotlib/projections/geo.py b/lib/matplotlib/projections/geo.py index 182266ca0a62..a5d8e3a50a75 100644 --- a/lib/matplotlib/projections/geo.py +++ b/lib/matplotlib/projections/geo.py @@ -271,16 +271,16 @@ def transform(self, ll): cos_latitude = np.cos(latitude) alpha = np.arccos(cos_latitude * np.cos(half_long)) - # Mask this array, or we'll get divide-by-zero errors + # Mask this array or we'll get divide-by-zero errors alpha = ma.masked_where(alpha == 0.0, alpha) + # The numerators also need to be masked so that masked + # division will be invoked. # We want unnormalized sinc. numpy.sinc gives us normalized sinc_alpha = ma.sin(alpha) / alpha - x = (cos_latitude * np.sin(half_long)) / sinc_alpha - y = (np.sin(latitude) / sinc_alpha) - x.set_fill_value(0.0) - y.set_fill_value(0.0) - return np.concatenate((x.filled(), y.filled()), 1) + x = (cos_latitude * ma.sin(half_long)) / sinc_alpha + y = (ma.sin(latitude) / sinc_alpha) + return np.concatenate((x.filled(0), y.filled(0)), 1) transform.__doc__ = Transform.transform.__doc__ transform_non_affine = transform