diff --git a/Changelog b/Changelog
index ecfe79ecf..dfb78096b 100644
--- a/Changelog
+++ b/Changelog
@@ -1,4 +1,5 @@
 version 1.0.3 (not yet released)
+           * added hexbin method.
            * drawmapboundary now uses axes bgcolor as default fill_color. If
              no color fill wanted, set fill_color='none' (a string).
            * clip coastlines for nplaea,npaeqd,splaea,spaeqd in stereographic
diff --git a/examples/hexbin_demo.py b/examples/hexbin_demo.py
index a7f1fd854..f47324d45 100644
--- a/examples/hexbin_demo.py
+++ b/examples/hexbin_demo.py
@@ -6,8 +6,9 @@
 
 # create north polar stereographic basemap
 m = Basemap(lon_0=270, boundinglat=20, projection='npstere',round=True)
+#m = Basemap(lon_0=-105,lat_0=40,projection='ortho')
 
-# number of points to plot.
+# number of points, bins to plot.
 npts = 10000
 bins = 40
 # generate random points on a sphere,
@@ -16,13 +17,16 @@
 # http://mathworld.wolfram.com/SpherePointPicking.html
 u = uniform(0.,1.,size=npts)
 v = uniform(0.,1.,size=npts)
-lons1 = 360.*u
-lats1 = (180./np.pi)*np.arccos(2*v-1) - 90.
+lons = 360.*u
+lats = (180./np.pi)*np.arccos(2*v-1) - 90.
 # toss points outside of map region.
-lats = np.compress(lats1 > 20, lats1)
-lons = np.compress(lats1 > 20, lons1)
+lats = np.compress(lats > 20, lats)
+lons = np.compress(lats > 20, lons)
 # convert to map projection coordinates.
-x, y = m(lons, lats)
+x1, y1 = m(lons, lats)
+# remove points outside projection limb.
+x = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), x1)
+y = np.compress(np.logical_or(x1 < 1.e20,y1 < 1.e20), y1)
 # function to plot at those points.
 xscaled = 4.*(x-0.5*(m.xmax-m.xmin))/m.xmax
 yscaled = 4.*(y-0.5*(m.ymax-m.ymin))/m.ymax
@@ -31,7 +35,7 @@
 # make plot using hexbin
 fig = plt.figure(figsize=(12,5))
 ax = fig.add_subplot(121)
-CS = plt.hexbin(x,y,C=z,gridsize=bins,cmap=plt.cm.jet)
+CS = m.hexbin(x,y,C=z,gridsize=bins,cmap=plt.cm.jet)
 # draw coastlines, lat/lon lines.
 m.drawcoastlines()
 m.drawparallels(np.arange(0,81,20))
@@ -41,16 +45,16 @@
 
 # use histogram2d instead of hexbin.
 ax = fig.add_subplot(122)
-#m = Basemap(lon_0=270, boundinglat=20, projection='npstere',round=True)
+# remove points outside projection limb.
 bincount, xedges, yedges = np.histogram2d(x, y, bins=bins)
 mask = bincount == 0
 # reset zero values to one to avoid divide-by-zero
 bincount = np.where(bincount == 0, 1, bincount)
 H, xedges, yedges = np.histogram2d(x, y, bins=bins, weights=z)
 H = np.ma.masked_where(mask, H/bincount)
-# set color of masked values to white (hexbin does this by default)
+# set color of masked values to axes background (hexbin does this by default)
 palette = plt.cm.jet
-palette.set_bad('white', 1.0)
+palette.set_bad(ax.get_axis_bgcolor(), 1.0)
 CS = m.pcolormesh(xedges,yedges,H.T,shading='flat',cmap=palette)
 # draw coastlines, lat/lon lines.
 m.drawcoastlines()
diff --git a/lib/mpl_toolkits/basemap/__init__.py b/lib/mpl_toolkits/basemap/__init__.py
index 3a6fee864..b8224d608 100644
--- a/lib/mpl_toolkits/basemap/__init__.py
+++ b/lib/mpl_toolkits/basemap/__init__.py
@@ -1513,16 +1513,22 @@ def fillcontinents(self,color='0.8',lake_color=None,ax=None,zorder=None):
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         # clip continent polygons for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            for poly in polys:
-                poly.set_clip_path(c)
+        if self.round: polys = self._clipcircle(ax,polys)
         return polys
 
+    def _clipcircle(self,ax,coll):
+        c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
+            radius=0.5*(self.xmax-self.xmin),fc='none')
+        if c not in ax.patches:
+            p = ax.add_patch(c)
+            p.set_clip_on(False)
+        try:
+            coll.set_clip_path(c)
+        except:
+            for item in coll:
+                item.set_clip_path(c)
+        return coll
+
     def drawcoastlines(self,linewidth=1.,color='k',antialiased=1,ax=None,zorder=None):
         """
         Draw coastlines.
@@ -1554,13 +1560,7 @@ def drawcoastlines(self,linewidth=1.,color='k',antialiased=1,ax=None,zorder=None
         if zorder is not None:
             coastlines.set_zorder(zorder)
         # clip coastlines for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            coastlines.set_clip_path(c)
+        if self.round: coastlines = self._clipcircle(ax,coastlines)
         ax.add_collection(coastlines)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
@@ -1603,13 +1603,7 @@ def drawcountries(self,linewidth=0.5,color='k',antialiased=1,ax=None,zorder=None
             countries.set_zorder(zorder)
         ax.add_collection(countries)
         # clip countries for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            countries.set_clip_path(c)
+        if self.round: countries = self._clipcircle(ax,countries)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return countries
@@ -1651,13 +1645,7 @@ def drawstates(self,linewidth=0.5,color='k',antialiased=1,ax=None,zorder=None):
             states.set_zorder(zorder)
         ax.add_collection(states)
         # clip states for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            states.set_clip_path(c)
+        if self.round: states = self._clipcircle(ax,states)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return states
@@ -1699,13 +1687,7 @@ def drawrivers(self,linewidth=0.5,color='k',antialiased=1,ax=None,zorder=None):
             rivers.set_zorder(zorder)
         ax.add_collection(rivers)
         # clip rivers for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            rivers.set_clip_path(c)
+        if self.round: rivers = self._clipcircle(ax,rivers)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return rivers
@@ -1869,13 +1851,7 @@ def readshapefile(self,shapefile,name,drawbounds=True,zorder=None,
                lines.set_zorder(zorder)
             ax.add_collection(lines)
             # clip boundaries for round polar plots.
-            if self.round:
-                c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                    radius=0.5*(self.xmax-self.xmin),fc='none')
-                if c not in ax.patches:
-                    p = ax.add_patch(c)
-                    p.set_clip_on(False)
-                lines.set_clip_path(c)
+            if self.round: lines = self._clipcircle(ax,lines)
             # set axes limits to fit map region.
             self.set_axes_limits(ax=ax)
             info = info + (lines,)
@@ -2516,13 +2492,7 @@ def tissot(self,lon_0,lat_0,radius_deg,npts,ax=None,**kwargs):
         poly = Polygon(seg,**kwargs)
         ax.add_patch(poly)
         # clip polygons for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            poly.set_clip_path(c)
+        if self.round: poly = self._clipcircle(ax,poly)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return poly
@@ -2903,13 +2873,7 @@ def scatter(self, *args, **kwargs):
         if plt:
             plt.sci(ret)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            ret.set_clip_path(c)
+        if self.round: ret = self._clipcircle(ax,ret)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return ret
@@ -2936,13 +2900,7 @@ def plot(self, *args, **kwargs):
             raise
         ax.hold(b)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            ret.set_clip_path(c)
+        if self.round: ret = self._clipcircle(ax,ret)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return ret
@@ -2981,13 +2939,7 @@ def imshow(self, *args, **kwargs):
         if plt:
             plt.sci(ret)
         # clip image for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            ret.set_clip_path(c)
+        if self.round: ret = self._clipcircle(ax,ret)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return ret
@@ -3055,13 +3007,7 @@ def pcolor(self,x,y,data,tri=False,**kwargs):
         if plt:
             plt.sci(ret)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            ret.set_clip_path(c)
+        if self.round: ret = self._clipcircle(ax,ret)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         if self.round:
@@ -3094,13 +3040,7 @@ def pcolormesh(self,x,y,data,**kwargs):
         if plt:
             plt.sci(ret)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            ret.set_clip_path(c)
+        if self.round: ret = self._clipcircle(ax,ret)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         if self.round:
@@ -3108,6 +3048,53 @@ def pcolormesh(self,x,y,data,**kwargs):
             ax.set_frame_on(False)
         return ret
 
+    def hexbin(self,x,y,**kwargs):
+        """
+        Make a hexagonal binning plot of x versus y, where x, y are 1-D
+        sequences of the same length, N. If C is None (the default), this is a
+        histogram of the number of occurences of the observations at
+        (x[i],y[i]).
+
+        If C is specified, it specifies values at the coordinate (x[i],y[i]).
+        These values are accumulated for each hexagonal bin and then reduced
+        according to reduce_C_function, which defaults to numpy?s mean function
+        (np.mean). (If C is specified, it must also be a 1-D sequence of the
+        same length as x and y.)
+
+        x, y and/or C may be masked arrays, in which case only unmasked points
+        will be plotted.
+
+        (see matplotlib.pyplot.hexbin documentation).
+
+        Extra keyword ``ax`` can be used to override the default axis instance.
+
+        Other \**kwargs passed on to matplotlib.pyplot.hexbin
+        """
+        ax, plt = self._ax_plt_from_kw(kwargs)
+        # allow callers to override the hold state by passing hold=True|False
+        b = ax.ishold()
+        h = kwargs.pop('hold',None)
+        if h is not None:
+            ax.hold(h)
+        try:
+            # make x,y masked arrays
+            # (masked where data is outside of projection limb)
+            x = ma.masked_values(np.where(x > 1.e20,1.e20,x), 1.e20)
+            y = ma.masked_values(np.where(y > 1.e20,1.e20,y), 1.e20)
+            ret = ax.hexbin(x,y,**kwargs)
+        except:
+            ax.hold(b)
+            raise
+        ax.hold(b)
+        # reset current active image (only if pyplot is imported).
+        if plt:
+            plt.sci(ret)
+        # clip for round polar plots.
+        if self.round: ret = self._clipcircle(ax,ret)
+        # set axes limits to fit map region.
+        self.set_axes_limits(ax=ax)
+        return ret
+
     def contour(self,x,y,data,*args,**kwargs):
         """
         Make a contour plot over the map
@@ -3187,14 +3174,7 @@ def contour(self,x,y,data,*args,**kwargs):
         if plt and CS.get_array() is not None:
             plt.sci(CS)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            for cntr in CS.collections:
-                cntr.set_clip_path(c)
+        if self.round: CS.collections = self._clipcircle(ax,CS.collections)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return CS
@@ -3291,14 +3271,7 @@ def contourf(self,x,y,data,*args,**kwargs):
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            for cntr in CS.collections:
-                cntr.set_clip_path(c)
+        if self.round: CS.collections = self._clipcircle(ax,CS.collections)
         return CS
 
     def quiver(self, x, y, u, v, *args, **kwargs):
@@ -3325,13 +3298,7 @@ def quiver(self, x, y, u, v, *args, **kwargs):
         if plt is not None and ret.get_array() is not None:
             plt.sci(ret)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            ret.set_clip_path(c)
+        if self.round: ret = self._clipcircle(ax,ret)
         # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return ret
@@ -3376,15 +3343,11 @@ def barbs(self, x, y, u, v, *args, **kwargs):
         # we can't set the current image...
         #if plt is not None and ret.get_array() is not None:
         #    plt.sci(retnh)
-        # set axes limits to fit map region.
         # clip for round polar plots.
         if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            ret.set_clip_path(c)
+            retnh = self._clipcircle(ax,retnh)
+            retsh = self._clipcircle(ax,retsh)
+        # set axes limits to fit map region.
         self.set_axes_limits(ax=ax)
         return retnh,retsh
 
@@ -3524,13 +3487,7 @@ def drawlsmask(self,land_color="0.8",ocean_color="w",lsmask=None,
         # plot mask as rgba image.
         im = self.imshow(rgba,interpolation='nearest',ax=ax,**kwargs)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            im.set_clip_path(c)
+        if self.round: im = self._clipcircle(ax,im)
         return im
 
     def bluemarble(self,ax=None,scale=None,**kwargs):
@@ -3743,13 +3700,7 @@ def warpimage(self,image="bluemarble",scale=None,**kwargs):
             # bmproj True, no interpolation necessary.
             im = self.imshow(self._bm_rgba,ax=ax,**kwargs)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            im.set_clip_path(c)
+        if self.round: im = self._clipcircle(ax,im)
         return im
 
     def drawmapscale(self,lon,lat,lon0,lat0,length,barstyle='simple',\
@@ -4025,14 +3976,7 @@ def nightshade(self,date,color="k",delta=0.25,alpha=0.5,ax=None,zorder=2):
         for c in CS.collections:
             c.set_zorder(zorder)
         # clip for round polar plots.
-        if self.round:
-            c = Circle((0.5*(self.xmax+self.xmin),0.5*(self.ymax+self.ymin)),
-                radius=0.5*(self.xmax-self.xmin),fc='none')
-            if c not in ax.patches:
-                p = ax.add_patch(c)
-                p.set_clip_on(False)
-            for cntr in CS.collections:
-                cntr.set_clip_path(c)
+        if self.round: CS.collections = self._clipcircle(ax,CS.collections)
         return CS
 
     def _check_ax(self):
diff --git a/lib/mpl_toolkits/basemap/geodesic.py b/lib/mpl_toolkits/basemap/geodesic.py
index e4ca0666d..832db15f7 100644
--- a/lib/mpl_toolkits/basemap/geodesic.py
+++ b/lib/mpl_toolkits/basemap/geodesic.py
@@ -122,14 +122,18 @@ def SinCosSeries(sinp, sinx, cosx,  c, n):
       # Unroll loop x 2, so accumulators return to their original role
       k -= 1; y1 = ar * y0 - y1 + c[k]
       k -= 1; y0 = ar * y1 - y0 + c[k]
-    return ( 2 * sinx * cosx * y0 if sinp # sin(2 * x) * y0
-             else cosx * (y0 - y1) )      # cos(x) * (y0 - y1)
+    if sinp: # sin(2 * x) * y0
+      return 2 * sinx * cosx * y0
+    else:    # cos(x) * (y0 - y1)
+      return cosx * (y0 - y1)
   SinCosSeries = staticmethod(SinCosSeries)
 
   def AngNormalize(x):
-    # Place angle in [-180, 180).  Assumes x is in [-540, 540).
-    return (x - 360 if x >= 180 else
-            (x + 360 if x < -180 else x))
+    # Place angle in [-180, 180); multiples of 360
+    if x < -180 or x >= 180:
+      return ((x + 180) % 360) - 180
+    else:
+      return x
   AngNormalize = staticmethod(AngNormalize)
 
   def AngRound(x):
@@ -141,8 +145,9 @@ def AngRound(x):
     z = 0.0625                  # 1/16
     y = abs(x)
     # The compiler mustn't "simplify" z - (z - y) to y
-    y =  z - (z - y) if y < z else y
-    return -y if x < 0 else y
+    if y < z:
+      y =  z - (z - y)
+    return cmp(x, 0) * y
   AngRound = staticmethod(AngRound)
 
   def SinCosNorm(sinx, cosx):
@@ -171,11 +176,13 @@ def Astroid(x, y):
         # Pick the sign on the sqrt to maximize abs(T3).  This minimizes loss
         # of precision due to cancellation.  The result is unchanged because
         # of the way the T is used in definition of u.
-        T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
+        T3 += cmp(T3, 0) * math.sqrt(disc) # T3 = (r * t)^3
         # N.B. cbrt always returns the real root.  cbrt(-8) = -2.
         T = Math.cbrt(T3)       # T = r * t
         # T can be zero; but then r2 / T -> 0.
-        u += T + (r2 / T if T != 0 else 0)
+        u += T
+        if T != 0:
+          u += (r2 / T)
       else:
         # T is complex, but the way u is defined the result is real.
         ang = math.atan2(math.sqrt(-disc), -(S + r3))
@@ -184,7 +191,11 @@ def Astroid(x, y):
         u += 2 * r * math.cos(ang / 3)
       v = math.sqrt(Math.sq(u) + q)  # guaranteed positive
       # Avoid loss of accuracy when u < 0.
-      uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
+      # u+v, guaranteed positive
+      if u < 0:
+        uv = q / (v - u)
+      else:
+        uv = u + v
       w = (uv - q) / (2 * v)               # positive?
       # Rearrange expression for k to avoid loss of accuracy due to
       # subtraction.  Division by 0 not possible because uv > 0, w >= 0.
@@ -263,18 +274,26 @@ def __init__(self, a, f):
     """
 
     self._a = float(a)
-    self._f = float(f) if f <= 1 else 1.0/f
+    if f <= 1:
+      self._f = float(f)
+    else:
+      self._f = 1.0/f
     self._f1 = 1 - self._f
     self._e2 = self._f * (2 - self._f)
     self._ep2 = self._e2 / Math.sq(self._f1) # e2 / (1 - e2)
     self._n = self._f / ( 2 - self._f)
     self._b = self._a * self._f1
     # authalic radius squared
-    self._c2 = (Math.sq(self._a) + Math.sq(self._b) *
-                (1 if self._e2 == 0 else
-                 (Math.atanh(math.sqrt(self._e2)) if self._e2 > 0 else
-                  math.atan(math.sqrt(-self._e2))) /
-                 math.sqrt(abs(self._e2))))/2
+    if self._e2 == 0:
+      self._c2 = (Math.sq(self._a) + Math.sq(self._b))/2
+    elif self._e2 > 0:
+      self._c2 = (Math.sq(self._a) +
+                  Math.sq(self._b) * Math.atanh(math.sqrt(self._e2)) /
+                  math.sqrt(abs(self._e2)))/2
+    else: # self._e2 < 0:
+      self._c2 = (Math.sq(self._a) +
+                  Math.sq(self._b) * math.atan(math.sqrt(-self._e2)) /
+                  math.sqrt(abs(self._e2)))/2
     self._etol2 = Geodesic.tol2_ / max(0.1, math.sqrt(abs(self._e2)))
     if not(Math.isfinite(self._a) and self._a > 0):
       raise ValueError("Major radius is not positive")
@@ -443,14 +462,17 @@ def InverseStart(self, sbet1, cbet1, sbet2, cbet2, lam12,
     sbet12a += cbet2 * sbet1
 
     shortline = cbet12 >= 0 and sbet12 < 0.5 and lam12 <= math.pi / 6
-    omg12 =  (lam12 / math.sqrt(1 - self._e2 * Math.sq(cbet1)) if shortline
-              else lam12)
+    if shortline:
+      omg12 = lam12 / math.sqrt(1 - self._e2 * Math.sq(cbet1))
+    else:
+      omg12 = lam12
     somg12 = math.sin(omg12); comg12 = math.cos(omg12)
 
     salp1 = cbet2 * somg12
-    calp1 = (
-      sbet12 + cbet2 * sbet1 * Math.sq(somg12) / (1 + comg12)  if comg12 >= 0
-      else sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12))
+    if comg12 >= 0:
+      calp1 = sbet12 + cbet2 * sbet1 * Math.sq(somg12) / (1 + comg12)
+    else:
+      calp1 = sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12)
 
     ssig12 = math.hypot(salp1, calp1)
     csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
@@ -492,8 +514,10 @@ def InverseStart(self, sbet1, cbet1, sbet2, cbet2, lam12,
           self._n, math.pi + bet12a, sbet1, -cbet1, sbet2, cbet2,
           cbet1, cbet2, dummy, False, C1a, C2a)
         x = -1 + m12a/(self._f1 * cbet1 * cbet2 * m0 * math.pi)
-        betscale = (sbet12a / x if x < -real(0.01)
-                    else -self._f * Math.sq(cbet1) * math.pi)
+        if x < -real(0.01):
+          betscale = sbet12a / x
+        else:
+          betscale = -self._f * Math.sq(cbet1) * math.pi
         lamscale = betscale / cbet1
         y = (lam12 - math.pi) / lamscale
 
@@ -502,7 +526,10 @@ def InverseStart(self, sbet1, cbet1, sbet2, cbet2, lam12,
         if self._f >= 0:
           salp1 = min(1.0, -x); calp1 = - math.sqrt(1 - Math.sq(salp1))
         else:
-          calp1 = max((0.0 if x > -Geodesic.tol1_ else -1.0),  x)
+          if x > -Geodesic.tol1_:
+            calp1 = max(0.0, x)
+          else:
+            calp1 = max(-1.0, x)
           salp1 = math.sqrt(1 - Math.sq(calp1))
       else:
         # Estimate alp1, by solving the astroid problem.
@@ -540,8 +567,10 @@ def InverseStart(self, sbet1, cbet1, sbet2, cbet2, lam12,
         #
         # Because omg12 is near pi, estimate work with omg12a = pi - omg12
         k = Geodesic.Astroid(x, y)
-        omg12a = lamscale * ( -x * k/(1 + k) if self._f >= 0
-                               else -y * (1 + k)/k )
+        if self._f >= 0:
+          omg12a = lamscale * -x * k/(1 + k)
+        else:
+          omg12a = lamscale * -y * (1 + k)/k
         somg12 = math.sin(omg12a); comg12 = -math.cos(omg12a)
         # Update spherical estimate of alp1 using omg12 instead of lam12
         salp1 = cbet2 * somg12
@@ -576,15 +605,23 @@ def Lambda12(self, sbet1, cbet1, sbet2, cbet2, salp1, calp1, diffp,
     # about this case, since this can yield singularities in the Newton
     # iteration.
     # sin(alp2) * cos(bet2) = sin(alp0)
-    salp2 = salp0 / cbet2 if cbet2 != cbet1 else salp1
+    if cbet2 != cbet1:
+      salp2 = salp0 / cbet2
+    else:
+      salp2 = salp1
     # calp2 = sqrt(1 - sq(salp2))
     #       = sqrt(sq(calp0) - sq(sbet2)) / cbet2
     # and subst for calp0 and rearrange to give (choose positive sqrt
     # to give alp2 in [0, pi/2]).
-    calp2 = (math.sqrt(Math.sq(calp1 * cbet1) +
-                       ((cbet2 - cbet1) * (cbet1 + cbet2) if cbet1 < -sbet1
-                        else (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2
-             if cbet2 != cbet1 or abs(sbet2) != -sbet1 else abs(calp1))
+    if cbet2 != cbet1 or abs(sbet2) != -sbet1:
+      if cbet1 < -sbet1:
+        calp2 = math.sqrt(Math.sq(calp1 * cbet1) +
+                          (cbet2 - cbet1) * (cbet1 + cbet2)) / cbet2
+      else:
+        calp2 = math.sqrt(Math.sq(calp1 * cbet1) +
+                          (sbet1 - sbet2) * (sbet1 + sbet2)) / cbet2
+    else:
+      calp2 = abs(calp1)
     # tan(bet2) = tan(sig2) * cos(alp2)
     # tan(omg2) = sin(alp0) * tan(sig2).
     ssig2 = sbet2; somg2 = salp0 * sbet2
@@ -633,7 +670,10 @@ def GenInverse(self, lat1, lon1, lat2, lon2, outmask):
     # Not sure this is necessary...
     lon12 = Geodesic.AngRound(lon12)
     # Make longitude difference positive.
-    lonsign = 1 if lon12 >= 0 else -1
+    if lon12 >= 0:
+      lonsign = 1
+    else:
+      lonsign = -1
     lon12 *= lonsign
     if lon12 == 180:
       lonsign = 1
@@ -641,12 +681,17 @@ def GenInverse(self, lat1, lon1, lat2, lon2, outmask):
     lat1 = Geodesic.AngRound(lat1)
     lat2 = Geodesic.AngRound(lat2)
     # Swap points so that point with higher (abs) latitude is point 1
-    swapp = 1 if abs(lat1) >= abs(lat2) else -1
-    if swapp < 0:
+    if abs(lat1) >= abs(lat2):
+      swapp = 1
+    else:
+      swapp = -1
       lonsign *= -1
       lat2, lat1 = lat1, lat2
     # Make lat1 <= 0
-    latsign = 1 if lat1 < 0 else -1
+    if lat1 < 0:
+      latsign = 1
+    else:
+      latsign = -1
     lat1 *= latsign
     lat2 *= latsign
     # Now we have
@@ -666,13 +711,19 @@ def GenInverse(self, lat1, lon1, lat2, lon2, outmask):
     phi = lat1 * Math.degree
     # Ensure cbet1 = +epsilon at poles
     sbet1 = self._f1 * math.sin(phi)
-    cbet1 = Geodesic.tiny_ if lat1 == -90 else math.cos(phi)
+    if lat1 == -90:
+      cbet1 = Geodesic.tiny_
+    else:
+      cbet1 = math.cos(phi)
     sbet1, cbet1 = Geodesic.SinCosNorm(sbet1, cbet1)
 
     phi = lat2 * Math.degree
     # Ensure cbet2 = +epsilon at poles
     sbet2 = self._f1 * math.sin(phi)
-    cbet2 = Geodesic.tiny_ if abs(lat2) == 90 else math.cos(phi)
+    if abs(lat2) == 90:
+      cbet2 = Geodesic.tiny_
+    else:
+      cbet2 = math.cos(phi)
     sbet2, cbet2 = Geodesic.SinCosNorm(sbet2, cbet2)
 
     # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
@@ -685,13 +736,19 @@ def GenInverse(self, lat1, lon1, lat2, lon2, outmask):
 
     if cbet1 < -sbet1:
       if cbet2 == cbet1:
-        sbet2 = sbet1 if sbet2 < 0 else -sbet1
+        if sbet2 < 0:
+          sbet2 = sbet1
+        else:
+          sbet2 = -sbet1
     else:
       if abs(sbet2) == -sbet1:
         cbet2 = cbet1
 
     lam12 = lon12 * Math.degree
-    slam12 = 0 if lon12 == 180 else math.sin(lam12)
+    if lon12 == 180:
+      slam12 = 0.0
+    else:
+      slam12 = math.sin(lam12)
     clam12 = math.cos(lam12)      # lon12 == 90 isn't interesting
 
     # real a12, sig12, calp1, salp1, calp2, salp2
@@ -968,10 +1025,12 @@ def Inverse(self, lat1, lon1, lat2, lon2, outmask = DISTANCE | AZIMUTH):
   # return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
   def GenDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask):
     from geodesicline import GeodesicLine
-    line = GeodesicLine(
-      self, lat1, lon1, azi1,
-      # Automatically supply DISTANCE_IN if necessary
-      outmask | ( Geodesic.NONE if arcmode else Geodesic.DISTANCE_IN))
+    # Automatically supply DISTANCE_IN if necessary
+    if arcmode:
+      dist_mask = Geodesic.NONE
+    else:
+      dist_mask = Geodesic.DISTANCE_IN
+    line = GeodesicLine(self, lat1, lon1, azi1, outmask | dist_mask)
     return line.GenPosition(arcmode, s12_a12, outmask)
 
   def Direct(self, lat1, lon1, azi1, s12,
diff --git a/lib/mpl_toolkits/basemap/geodesicline.py b/lib/mpl_toolkits/basemap/geodesicline.py
index f695a9e3a..e98d6995d 100644
--- a/lib/mpl_toolkits/basemap/geodesicline.py
+++ b/lib/mpl_toolkits/basemap/geodesicline.py
@@ -32,7 +32,7 @@ class GeodesicLine(object):
   """Points on a geodesic path"""
 
   def __init__(self, geod, lat1, lon1, azi1, caps = GeodesicCapability.ALL):
-    from mpl_toolkits.basemap.geodesic import Geodesic
+    from pyproj.geodesic import Geodesic
     self._a = geod._a
     self._f = geod._f
     self._b = geod._b
@@ -51,13 +51,22 @@ def __init__(self, geod, lat1, lon1, azi1, caps = GeodesicCapability.ALL):
     alp1 = azi1 * Math.degree
     # Enforce sin(pi) == 0 and cos(pi/2) == 0.  Better to face the ensuing
     # problems directly than to skirt them.
-    self._salp1 = 0 if     azi1  == -180 else math.sin(alp1)
-    self._calp1 = 0 if abs(azi1) ==   90 else math.cos(alp1)
+    if azi1 == -180:
+      self._salp1 = 0
+    else:
+      self._salp1 = math.sin(alp1)
+    if abs(azi1) == 90:
+      self._calp1 = 0
+    else:
+      self._calp1 = math.cos(alp1)
     # real cbet1, sbet1, phi
     phi = lat1 * Math.degree
     # Ensure cbet1 = +epsilon at poles
     sbet1 = self._f1 * math.sin(phi)
-    cbet1 = Geodesic.tiny_ if abs(lat1) == 90 else math.cos(phi)
+    if abs(lat1) == 90:
+      cbet1 = Geodesic.tiny_
+    else:
+      cbet1 = math.cos(phi)
     sbet1, cbet1 = Geodesic.SinCosNorm(sbet1, cbet1)
 
     # Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
@@ -75,8 +84,10 @@ def __init__(self, geod, lat1, lon1, azi1, caps = GeodesicCapability.ALL):
     # No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
     # With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
     self._ssig1 = sbet1; self._somg1 = self._salp0 * sbet1
-    self._csig1 = self._comg1 = (cbet1 * self._calp1
-                                 if sbet1 != 0 or self._calp1 != 0 else 1)
+    if sbet1 != 0 or self._calp1 != 0:
+      self._csig1 = self._comg1 = cbet1 * self._calp1
+    else:
+      self._csig1 = self._comg1 = 1.0
     # sig1 in (-pi, pi]
     self._ssig1, self._csig1 = Geodesic.SinCosNorm(self._ssig1, self._csig1)
     self._somg1, self._comg1 = Geodesic.SinCosNorm(self._somg1, self._comg1)
@@ -126,7 +137,7 @@ def __init__(self, geod, lat1, lon1, azi1, caps = GeodesicCapability.ALL):
   # return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
   def GenPosition(self, arcmode, s12_a12, outmask):
 
-    from mpl_toolkits.basemap.geodesic import Geodesic
+    from pyproj.geodesic import Geodesic
     a12 = lat2 = lon2 = azi2 = s12 = m12 = M12 = M21 = S12 = Math.nan
     outmask &= self._caps & Geodesic.OUT_ALL
     if not (arcmode or (self._caps & Geodesic.DISTANCE_IN & Geodesic.OUT_ALL)):
@@ -140,8 +151,14 @@ def GenPosition(self, arcmode, s12_a12, outmask):
       sig12 = s12_a12 * Math.degree
       s12a = abs(s12_a12)
       s12a -= 180 * math.floor(s12a / 180)
-      ssig12 = 0 if s12a ==  0 else math.sin(sig12)
-      csig12 = 0 if s12a == 90 else math.cos(sig12)
+      if s12a == 0:
+        ssig12 = 0.0
+      else:
+        ssig12 = math.sin(sig12)
+      if s12a == 90:
+        csig12 = 0.0
+      else:
+        csig12 = math.cos(sig12)
     else:
       # Interpret s12_a12 as distance
       tau12 = s12_a12 / (self._b * (1 + self._A1m1))
@@ -180,7 +197,10 @@ def GenPosition(self, arcmode, s12_a12, outmask):
                   comg2 * self._comg1 + somg2 * self._somg1)
 
     if outmask & Geodesic.DISTANCE:
-      s12 = self._b * ((1 + self._A1m1) * sig12 + AB1) if arcmode else s12_a12
+      if arcmode:
+        s12 = self._b * ((1 + self._A1m1) * sig12 + AB1)
+      else:
+        s12 = s12_a12
 
     if outmask & Geodesic.LONGITUDE:
       lam12 = omg12 + self._A3c * (
@@ -242,14 +262,20 @@ def GenPosition(self, arcmode, s12_a12, outmask):
         # else
         #   csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
         # No need to normalize
-        salp12 = self._calp0 * self._salp0 * (
-          self._csig1 * (1 - csig12) + ssig12 * self._ssig1 if csig12 <= 0
-          else ssig12 * (self._csig1 * ssig12 / (1 + csig12) + self._ssig1))
+        if csig12 > 0:
+          salp12 = self._calp0 * self._salp0 * (
+            ssig12 * (self._csig1 * ssig12 / (1 + csig12) + self._ssig1))
+        else:
+          salp12 = self._calp0 * self._salp0 * (
+            self._csig1 * (1 - csig12) + ssig12 * self._ssig1)
         calp12 = (Math.sq(self._salp0) +
                   Math.sq(self._calp0) * self._csig1 * csig2)
       S12 = self._c2 * math.atan2(salp12, calp12) + self._A4 * (B42 - self._B41)
 
-    a12 = s12_a12 if arcmode else sig12 / Math.degree
+    if arcmode:
+      a12 = s12_a12
+    else:
+      a12 = sig12 / Math.degree
     return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
 
   def Position(self, s12,
@@ -286,7 +312,7 @@ def Position(self, s12,
       Geodesic.ALL
     """
 
-    from mpl_toolkits.basemap.geodesic import Geodesic
+    from pyproj.geodesic import Geodesic
     Geodesic.CheckDistance(s12)
     result = {'lat1': self._lat1, 'lon1': self._lon1, 'azi1': self._azi1,
               's12': s12}
diff --git a/lib/mpl_toolkits/basemap/polygonarea.py b/lib/mpl_toolkits/basemap/polygonarea.py
index 0d0a68137..4b9b9e3a7 100644
--- a/lib/mpl_toolkits/basemap/polygonarea.py
+++ b/lib/mpl_toolkits/basemap/polygonarea.py
@@ -23,23 +23,29 @@ class PolygonArea(object):
   def transit(lon1, lon2):
     # Return 1 or -1 if crossing prime meridian in east or west direction.
     # Otherwise return zero.
-    from mpl_toolkits.basemap.geodesic import Geodesic
+    from pyproj.geodesic import Geodesic
     lon1 = Geodesic.AngNormalize(lon1)
     lon2 = Geodesic.AngNormalize(lon2)
     # treat lon12 = -180 as an eastward geodesic, so convert to 180.
     lon12 = -Geodesic.AngNormalize(lon1 - lon2) # In (-180, 180]
-    cross = (1 if lon1 < 0 and lon2 >= 0 and lon12 > 0
-             else (-1 if lon2 < 0 and lon1 >= 0 and lon12 < 0 else 0))
+    if lon1 < 0 and lon2 >= 0 and lon12 > 0:
+      cross = 1
+    elif lon2 < 0 and lon1 >= 0 and lon12 < 0:
+      cross = -1
+    else:
+      cross = 0
     return cross
   transit = staticmethod(transit)
 
   def __init__(self, earth, polyline = False):
-    from mpl_toolkits.basemap.geodesic import Geodesic
+    from pyproj.geodesic import Geodesic
     self._earth = earth
     self._area0 = 4 * math.pi * earth._c2
     self._polyline = polyline
-    self._mask = Geodesic.DISTANCE | (0 if self._polyline else Geodesic.AREA)
-    if not self._polyline: self._areasum = Accumulator()
+    self._mask = Geodesic.DISTANCE
+    if not self._polyline:
+      self._mask |= Geodesic.AREA
+      self._areasum = Accumulator()
     self._perimetersum = Accumulator()
     self.Clear()
 
@@ -84,7 +90,10 @@ def Compute(self, reverse, sign):
     tempsum.Add(S12)
     crossings = self._crossings + PolygonArea.transit(self._lon1, self._lon0)
     if crossings & 1:
-      tempsum.Add( (1 if tempsum < 0 else -1) * self._area0/2 )
+      if tempsum < 0:
+        tempsum.Add(self._area0 / 2)
+      else:
+        tempsum.Add(-self._area0 / 2)
     # area is with the clockwise sense.  If !reverse convert to
     # counter-clockwise convention.
     if not reverse: tempsum.Negate()
@@ -112,24 +121,38 @@ def TestCompute(self, lat, lon, reverse, sign):
       return 1, perimeter, area
 
     perimeter = self._perimetersum.Sum()
-    tempsum = 0 if self._polyline else self._areasum.Sum()
+    if self._polyline:
+      tempsum = 0
+      dl = [0]
+    else:
+      tempsum = self._areasum.Sum()
+      dl = [0, 1]
     crossings = self._crossings; num = self._num + 1
-    for i in ([0] if self._polyline else [0, 1]):
-      t, s12, t, t, t, t, t, S12 = self._earth.GenInverse(
-        self._lat1 if i == 0 else lat, self._lon1 if i == 0 else lon,
-        self._lat0 if i != 0 else lat, self._lon0 if i != 0 else lon,
-        self._mask)
+    for i in dl:
+      if i == 0:
+        ilon1 = self._lon1
+        ilon0 = lon
+        t, s12, t, t, t, t, t, S12 = self._earth.GenInverse(
+          self._lat1, self._lon1, lat, lon, self._mask)
+      else:
+        ilon1 = lon
+        ilon0 = self._lon0
+        t, s12, t, t, t, t, t, S12 = self._earth.GenInverse(
+          lat, lon, self._lat0, self._lon0, self._mask)
       perimeter += s12
       if not self._polyline:
         tempsum += S12
-        crossings += PolygonArea.transit(self._lon1 if i == 0 else lon,
-                                         self._lon0 if i != 0 else lon)
+        
+        crossings += PolygonArea.transit(ilon1, ilon0)
 
     if self._polyline:
       return num, perimeter, area
 
     if crossings & 1:
-      tempsum += (1 if tempsum < 0 else -1) * self._area0/2
+      if tempsum < 0:
+        tempsum += self._area0 / 2
+      else:
+        tempsum -= self._area0 / 2
     # area is with the clockwise sense.  If !reverse convert to
     # counter-clockwise convention.
     if not reverse: tempsum *= -1