diff --git a/lib/matplotlib/projections/geo.py b/lib/matplotlib/projections/geo.py
index 8a20f5c2768d..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
@@ -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