diff --git a/DataProcess/DataProcess.py b/DataProcess/DataProcess.py
index 9693a98c..d83113c3 100644
--- a/DataProcess/DataProcess.py
+++ b/DataProcess/DataProcess.py
@@ -312,10 +312,6 @@ def processData(self, restrictToWindfieldDomain=False):
         if self.ncflag:
             self.data['index'] = indicator
 
-        # ieast : parameter used in latLon2Azi
-        # FIXME: should be a config setting describing the input data.
-        ieast = 1
-
         # Determine the index of initial cyclone observations, excluding
         # those cyclones that have only one observation. This is used
         # for calculating initial bearing and speed
@@ -325,7 +321,7 @@ def processData(self, restrictToWindfieldDomain=False):
 
         # Calculate the bearing and distance (km) of every two
         # consecutive records using ll2azi
-        bear_, dist_ = maputils.latLon2Azi(lat, lon, ieast, azimuth=0)
+        bear_, dist_ = maputils.latLon2Azi(lat, lon)
         assert bear_.size == indicator.size - 1
         assert dist_.size == indicator.size - 1
         bear = np.empty(indicator.size, 'f')
diff --git a/Evaluate/interpolateTracks.py b/Evaluate/interpolateTracks.py
index eb3123db..63d25268 100644
--- a/Evaluate/interpolateTracks.py
+++ b/Evaluate/interpolateTracks.py
@@ -201,7 +201,7 @@ def interpolate(track, delta, interpolation_type=None):
         nrMax = interp1d(timestep, track.rMax, kind='linear')(newtime)
 
     if len(nLat) >= 2:
-        bear_, dist_ = latLon2Azi(nLat, nLon, 1, azimuth=0)
+        bear_, dist_ = latLon2Azi(nLat, nLon)
         nthetaFm = np.zeros(newtime.size, dtype=float)
         nthetaFm[:-1] = bear_
         nthetaFm[-1] = bear_[-1]
diff --git a/Utilities/interpTrack.py b/Utilities/interpTrack.py
index 3da0054a..f781c9a5 100644
--- a/Utilities/interpTrack.py
+++ b/Utilities/interpTrack.py
@@ -186,7 +186,7 @@ def interpolateTrack(configFile, trackFile, source, delta=0.1,
         npEnv = scint.interp1d(timestep, penv, kind='linear')(newtime)
         nrMax = scint.interp1d(timestep, rmax, kind='linear')(newtime)
 
-    bear_, dist_ = maputils.latLon2Azi(nLat, nLon, 1, azimuth=0)
+    bear_, dist_ = maputils.latLon2Azi(nLat, nLon)
     nthetaFm = numpy.zeros(newtime.size, 'f')
     nthetaFm[:-1] = bear_
     nthetaFm[-1] = bear_[-1]
diff --git a/Utilities/loadData.py b/Utilities/loadData.py
index ea3d5521..6061fc63 100644
--- a/Utilities/loadData.py
+++ b/Utilities/loadData.py
@@ -152,7 +152,7 @@ def __getattr__(self, key):
 
 """
 
-def getSpeedBearing(index, lon, lat, deltatime, ieast=1,
+def getSpeedBearing(index, lon, lat, deltatime,
                     missingValue=sys.maxsize):
     """
     Calculate the speed and bearing of a TC.
@@ -182,7 +182,7 @@ def getSpeedBearing(index, lon, lat, deltatime, ieast=1,
 
     """
 
-    bear_, dist_ = maputils.latLon2Azi(lat, lon, ieast, azimuth=0)
+    bear_, dist_ = maputils.latLon2Azi(lat, lon)
     assert bear_.size == index.size - 1
     assert dist_.size == index.size - 1
     bearing = np.zeros(index.size, 'f')
@@ -649,7 +649,7 @@ def ltmPressure(jdays, time, lon, lat, ncfile, ncvar='slp'):
     ncobj = nctools.ncLoadFile(ncfile)
     if ncvar not in ncobj.variables:
         raise KeyError(f"{ncfile} does not contain a variable called '{ncvar}'")
-    
+
     slpunits = getattr(ncobj.variables[ncvar], 'units')
 
     data = nctools.ncGetData(ncobj, ncvar)
@@ -679,9 +679,9 @@ def getPoci(penv, pcentre, lat, jdays, eps,
     :param lat: Latitude of storm (degrees).
     :param jdays: Julian day (day of year).
     :param eps: random variate. Retained as a constant for a single storm.
-    :param list coeffs: Coefficients of the model. Defaults based on 
-                        Southern Hemisphere data 
-                        (IBTrACS v03r06, 1981-2014). 
+    :param list coeffs: Coefficients of the model. Defaults based on
+                        Southern Hemisphere data
+                        (IBTrACS v03r06, 1981-2014).
 
     :returns: Revised estimate for the pressure of outermost closed isobar.
     """
@@ -700,7 +700,7 @@ def getPoci(penv, pcentre, lat, jdays, eps,
         assert len(penv) == len(pcentre)
         assert len(penv) == len(lat)
         assert len(penv) == len(jdays)
-      
+
     poci_model = coeffs[0] + coeffs[1]*penv + coeffs[2]*pcentre \
       + coeffs[3]*pcentre*pcentre + coeffs[4]*lat*lat + \
         coeffs[5]*np.cos(np.pi*2*jdays/365) + eps
@@ -716,9 +716,9 @@ def getPoci(penv, pcentre, lat, jdays, eps,
         poci_model = np.nan
     elif pcentre == missingValue:
         poci_model = np.nan
-    
+
     return poci_model
-    
+
 
 def filterPressure(pressure, inputPressureUnits='hPa',
                    missingValue=sys.maxsize):
diff --git a/Utilities/maputils.py b/Utilities/maputils.py
index 61a60e97..a3f57d40 100644
--- a/Utilities/maputils.py
+++ b/Utilities/maputils.py
@@ -18,20 +18,22 @@
 import math
 from . import metutils
 import warnings
+import pyproj
+
 try:
     from . import fmaputils
 except ImportError:
-    warnings.warn("Compiled maputils not found - defaulting to slower python wind models")
-
+    warnings.warn(
+        "Compiled maputils not found - defaulting to slower python wind models"
+    )
 
-# C weave code disabled for now.  The code speeds up the windfield interface module by ~6% but
-# does not appear to work on some systems.
-#from scipy import weave
-#from scipy.weave import converters
 
 logger = logging.getLogger(__name__)
 logger.addHandler(logging.NullHandler())
 
+geodesic = pyproj.Geod(ellps="WGS84")
+
+
 def xy2r(x, y):
     """
     Given x and y arrays, returns the distance between consecutive
@@ -46,19 +48,18 @@ def xy2r(x, y):
     :rtype: :class:`numpy.ndarray`
     """
 
-    #if len(x) != len(y):
+    # if len(x) != len(y):
     #   raise ArrayMismatch, "Input array sizes do not match"
     return np.sqrt(x**2 + y**2)
 
-def latLon2Azi(lat, lon, ieast=1, azimuth=0, wantdeg=True):
+
+def latLon2Azi(lat, lon, wantdeg=True):
     """
-    Returns the bearing and distance (in km) between consecutive
-    members of the array pair (lat,lon).
+    Returns the forward bearing and distance (in km) between consecutive
+    members of the array pair (lat, lon)
 
     :param lat: Latitudes of positions.
     :param lon: Longitudes of positions.
-    :param int ieast: 1 for longitudes increasing towards East, -1 for
-                      longitudes increasing towards West (default 1).
     :param float azimuth: Local coordinate system constructed with origin at
                           latr,lonr, X axis ('North') in direction of azimuth,
                           and Y axis such that X x Y = Z(down)
@@ -66,35 +67,13 @@ def latLon2Azi(lat, lon, ieast=1, azimuth=0, wantdeg=True):
     :param boolean wantdeg: If ``True`` return bearings as degrees, not radians.
 
     :returns: azimuth (+ve clockwise from north) and distance (in km).
-
     """
-    #if len(lat) != len(lon):
-    #   raise ArrayMismatch, "Input array sizes do not match"
 
-    xr = 0
-    yr = 0
-
-    yn, xe = latLon2XY(xr, yr, lat, lon, ieast, azimuth)
-    length = xy2r(yn, xe)
-    ### for azimuth calculation use atan2 which returns
-    ### angle from -pi to pi. Rules for getting azimuth are:
-    ### 1st quadrant (yn > 0 xe > 0):   0    <= angle <=  pi/2
-    ###      Rule -1* + pi/2 maps to:  pi/2               0
-    ### 2nd quadrant (yn < 0 xe > 0): -pi/2  <= angle <=  0
-    ###      Rule -1* + pi/2 maps to:  pi                 pi/2
-    ### 3rd quadrant (yn < 0 xe < 0): -pi    <= angle <= -pi/2
-    ###      Rule -1* + pi/2 maps to:  3pi/2              pi
-    ### 4th quadrant (yn > 0 xe < 0):  pi/2  <= angle <=  pi
-    ###      Rule -1* + 5pi/2 maps to: 2pi               3pi/2
-    ####################################################################
-    angle = np.arctan2(yn, xe) # yes, in that order
-    bearing = [theta2bearing(i) for i in angle]
-
-    # If bearing in degrees isexpected on return:
-    if wantdeg:
-        bearing = np.array([math.degrees(i) for i in bearing], 'f')
-
-    return bearing, length
+    fwdazimuth, _, distances = geodesic.inv(
+        lon[:-1], lat[:-1], lon[1:], lat[1:]
+        )
+    return fwdazimuth, distances / 1000.0
+
 
 def bear2LatLon(bearing, distance, oLon, oLat):
     """
@@ -109,19 +88,9 @@ def bear2LatLon(bearing, distance, oLon, oLat):
 
     :returns: new longitude and latitude (in degrees)
     """
-    radius = 6367.0 # Earth radius (km)
-    oLon = math.radians(oLon)
-    oLat = math.radians(oLat)
-    bear = math.radians(bearing)
-
-    nLat = math.asin(np.sin(oLat) * np.cos(distance / radius) + \
-            np.cos(oLat) * np.sin(distance / radius) * np.cos(bear))
-    aa = np.sin(bear) * np.sin(distance / radius) * np.cos(oLat)
-    bb = np.cos(distance / radius) - np.sin(oLat) * np.sin(nLat)
-
-    nLon = oLon + np.arctan2(aa, bb)
+    lon, lat, _ = geodesic.fwd(oLon, oLat, bearing, distance * 1000)
+    return lon, lat
 
-    return math.degrees(nLon), math.degrees(nLat)
 
 def latLon2XY(xr, yr, lat, lon, ieast=1, azimuth=0):
     """
@@ -148,17 +117,14 @@ def latLon2XY(xr, yr, lat, lon, ieast=1, azimuth=0):
               distance between consecutive points.
     """
 
-    #if len(lat) != len(lon):
-    #   raise ArrayMismatch, "Input array sizes do not match"
-
-    radius = 6367.0 # Earth radius (km)
+    radius = 6367.0  # Earth radius (km)
 
     lat = np.radians(lat)
     lon = np.radians(lon)
 
     # Is azimuth fixed or variable?
     if np.size(azimuth) == 1:
-        angle = np.radians(azimuth)*np.ones(lat.size - 1)
+        angle = np.radians(azimuth) * np.ones(lat.size - 1)
     else:
         angle = np.radians(azimuth)
 
@@ -172,47 +138,26 @@ def latLon2XY(xr, yr, lat, lon, ieast=1, azimuth=0):
 
     return xn, ye
 
+
 def distGC(lat, lon):
     """
-    Distance based on the great circle navigation between pairs of points.
-
-    :param lat: A pair of latitude values for the two points.
-    :param lon: A pair of longitude values for the two points.
-
-    :returns: Distance (in kilometres) between the two points, based on
-              great circle navigation.
-
-    Example::
-
-        >>> dist = distGC([-20, -40],[120,190])
-        6914.42
+    Distance based on great circle navigation between pairs of points
 
+    :param lat: _description_
+    :type lat: _type_
+    :param lon: _description_
+    :type lon: _type_
     """
-    radius = 6367.0 # Earth radius (km)
-
-    lat = np.radians(lat)
-    lon = np.radians(lon)
-
-    angular_distance = math.acos(math.sin(lat[0]) * math.sin(lat[1]) + \
-                       math.cos(lat[0]) * math.cos(lat[1]) * \
-                       math.cos(lon[0] - lon[1]))
-
-    return radius*angular_distance
-
+    return geodesic.line_lengths(lon, lat, radians=False)
 
 
-def gridLatLonDist(cLon, cLat, lonArray, latArray, units=None):
+def gridLatLonDistBear(cLon, cLat, lonArray, latArray):
     """
-    Generate a grid containing the spherical earth distance
+    Generate a grid containing the spherical earth distance and bearing
     of the points defined by (lonarray, latarray) from the
     point defined by (clon, clat).
     (lonarray,latarray) and (clon,clat) are in degrees.
-    Returns distances in km by default, other units specified by the
-    'units' kwarg.
-
-    Based on m_lldist.m by Rich Pawlowicz (rich@ocgy.ubc.ca)
-    Modified by Craig Arthur 2006-11-13
-
+    Returns distances in km, bearing in radians
 
     :param float cLon: Longitude of the point to measure the distance from.
     :param float cLat: Latitude of the point to measure the distance from.
@@ -220,217 +165,28 @@ def gridLatLonDist(cLon, cLat, lonArray, latArray, units=None):
                      grid over which distances will be calculated.
     :param latArray: 1-d array of latitude values that will define the
                      grid over which distances will be calculated.
-    :param str units: Units of distance to be returned (default is kilometre)
 
-    :returns: 2-d array containing the distance of the points defined in
-             ``lonArray`` and ``latArray`` from the point
+    :returns: 2 2-d arrays containing the bearing and distance of the points
+              defined in ``lonArray`` and ``latArray`` from the point
              (``cLon``, ``cLat``).
 
     Example::
 
         >>> lonArray = np.arange(90.,100.,0.1)
         >>> latArray = np.arange(-20.,-10.,0.1)
-        >>> dist = gridLatLonDist( 105., -15., lonArray, latArray, 'km')
-
-    """
-
-    # #CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-    # cLat_cos = 0.0
-    # cLat_sin = 0.0
-    # lat = empty(len(latArray), 'd')
-    # lon = empty(len(lonArray), 'd')
-
-    # dLon_sin = empty(len(lonArray), 'd')
-    # dLat_sin = empty(len(latArray), 'd')
-    # lat_cos = empty(len(latArray), 'd')
-
-    # dist = empty([len(latArray), len(lonArray)], 'd')
-
-    # code = """
-        # #include <math.h>
-
-        # double radius = 6367.0;
-        # double toRads = 0.017453292519943295;
-
-        # double cLon_ = cLon;
-        # double cLat_ = cLat;
-
-        # cLon_ = cLon_*toRads;
-        # cLat_ = cLat_*toRads;
-
-        # cLat_cos = cos(cLat_);
+        >>> bearing, distance = gridLatLonDistBear( 105., -15., lonArray, latArray)
 
-        # for (int i = 0; i < NlonArray[0]; ++i)
-        # {
-            # lon(i) = lonArray(i)*toRads;
-            # double dLon = (lon(i) - cLon_)/2.0;
-            # dLon_sin(i) = sin(dLon);
-        # }
-
-        # for (int i = 0; i < NlatArray[0]; ++i)
-        # {
-            # lat(i) = latArray(i)*toRads;
-            # lat_cos(i) = cos(lat(i));
-
-            # double dLat = (lat(i) - cLat_)/2.0;
-            # dLat_sin(i) = sin(dLat);
-        # }
-
-        # for (int j = 0; j < NlatArray[0]; ++j)
-        # {
-            # for (int i = 0; i < NlonArray[0]; ++i)
-            # {
-                 # double a = pow(dLat_sin(j), 2) + \
-                            # cLat_cos*lat_cos(j)*pow(dLon_sin(i), 2);
-                 # double c = 2.0*atan2(sqrt(fabs(a)), sqrt(1 - a));
-
-                 # dist(j, i) = radius*c;
-            # }
-        # }
-    # """
-    # err = weave.inline(code,
-                       # ['cLon', 'cLat', 'lonArray', 'latArray', 'lat', 'lon',
-                        # 'dLon_sin', 'dLat_sin', 'lat_cos', 'dist', 'cLat_cos'],
-                       # type_converters=converters.blitz,
-                       # compiler = 'gcc')
-    # #CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-
-    radius = 6367.0
-
-    lat = np.radians(latArray)
-    lon = np.radians(lonArray)
-
-    cLon = math.radians(cLon)
-    cLat = math.radians(cLat)
-    lon_, lat_ = np.meshgrid(lon, lat)
-
-    dLon = lon_ - cLon
-    dLat = lat_ - cLat
-
-    a = np.square(np.sin(dLat / 2.0)) + \
-        np.cos(cLat) * np.cos(lat_) * np.square(np.sin(dLon / 2.0))
-    c = 2.0 * np.arctan2(np.sqrt(np.absolute(a)), np.sqrt(1 - a))
-    dist = radius * c
-
-    dist = metutils.convert(dist, "km", units)
-
-    return dist
-
-def gridLatLonBear(cLon, cLat, lonArray, latArray):
     """
-    Generate a grid containing the bearing of the points defined by
-    (lonArray,latArray) from the point defined by (cLon,cLat).
-    (lonArray,latArray) and (cLon,cLat) are in degrees.
-    Returns bearing in radians.
-
-    :param float cLon: Longitude of the point to measure the distance from.
-    :param float cLat: Latitude of the point to measure the distance from.
-    :param lonArray: 1-d array of longitude values that will define the
-                     grid over which distances will be calculated.
-    :param latArray: 1-d array of latitude values that will define the
-                     grid over which distances will be calculated.
-    :returns: 2-d array containing the bearing (direction) of the points
-              defined in ``lonArray`` and ``latArray`` from the point
-              (``cLon``, ``cLat``)
-
-    Example::
-
-        >>> from maputils import gridLatLonBear
-        >>> import numpy as np
-        >>> lonArray = np.arange(90.,100.,0.1)
-        >>> latArray = np.arange(-20.,-10.,0.1)
-        >>> gridLatLonBear( 105., -15., lonArray, latArray)
-        array([[-1.94475949, -1.94659552, -1.94845671, ..., -2.36416927,
-                -2.37344337, -2.38290081],
-               [-1.93835542, -1.94015859, -1.94198663, ..., -2.35390045,
-                -2.36317282, -2.37263233],
-               [-1.93192776, -1.93369762, -1.93549204, ..., -2.34343069,
-                -2.35269718, -2.36215458],
-                ...,
-               [-1.29066433, -1.28850464, -1.28632113, ..., -0.84374983,
-                -0.83405688, -0.82416555],
-               [-1.28446304, -1.28227062, -1.28005406, ..., -0.83332654,
-                -0.82361918, -0.813717  ],
-               [-1.27828819, -1.27606348, -1.27381433, ..., -0.82310335,
-                -0.81338586, -0.80347714]])
-
-    """
-
-    # #CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-    # lat = empty(len(latArray), 'd')
-    # lon = empty(len(lonArray), 'd')
-
-    # dLon_sin = empty(len(lonArray), 'd')
-    # dLon_cos = empty(len(lonArray), 'd')
-    # lat_sin = empty(len(latArray), 'd')
-    # lat_cos = empty(len(latArray), 'd')
-
-    # bearing = empty([len(latArray), len(lonArray)], 'd')
-
-    # code = """
-        # #include <math.h>
-
-        # double toRads = 0.017453292519943295;
-
-        # double cLon_ = cLon;
-        # double cLat_ = cLat;
-
-        # cLon_ = cLon_*toRads;
-        # cLat_ = cLat_*toRads;
-
-        # double cLat_cos = cos(cLat_);
-        # double cLat_sin = sin(cLat_);
-
-        # for (int i = 0; i < NlonArray[0]; ++i)
-        # {
-            # lon(i) = lonArray(i)*toRads;
-            # double dLon = lon(i) - cLon_;
-            # dLon_sin(i) = sin(dLon);
-            # dLon_cos(i) = cos(dLon);
-        # }
-
-        # for (int i = 0; i < NlatArray[0]; ++i)
-        # {
-            # lat(i) = latArray(i)*toRads;
-            # lat_sin(i) = sin(lat(i));
-            # lat_cos(i) = cos(lat(i));
-        # }
-
-        # for (int j = 0; j < NlatArray[0]; ++j)
-        # {
-            # for (int i = 0; i < NlonArray[0]; ++i)
-            # {
-                # double alpha = dLon_sin(i)*lat_cos(j);
-                # double beta = (cLat_cos*lat_sin(j)) - (cLat_sin*lat_cos(j)*dLon_cos(i));
-
-                # bearing(j, i) = atan2(alpha, beta);
-            # }
-        # }
-    # """
-    # err = weave.inline(code,
-                       # ['cLon', 'cLat', 'lonArray', 'latArray', 'lat', 'lon',
-                        # 'dLon_sin', 'dLon_cos', 'lat_sin', 'lat_cos', 'bearing'],
-                       # type_converters=converters.blitz,
-                       # compiler = 'gcc')
-    # #CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
-
-    lat = np.radians(latArray)
-    lon = np.radians(lonArray)
-
-    cLon = math.radians(cLon)
-    cLat = math.radians(cLat)
-    lon_, lat_ = np.meshgrid(lon, lat)
-
-    dLon = lon_ - cLon
-    #dLat= lat_ - cLat
-
-    alpha = np.sin(dLon) * np.cos(lat_)
-    beta = np.cos(cLat) * np.sin(lat_) - \
-           np.sin(cLat) * np.cos(lat_) * np.cos(dLon)
-
-    bearing = np.arctan2(alpha, beta)
-
-    return bearing
+    rows, cols = len(latArray), len(lonArray)
+    longrid, latgrid = np.meshgrid(lonArray, latArray)
+    lonflat = longrid.flatten()
+    latflat = latgrid.flatten()
+    bflat, _, dflat = geodesic.inv(np.tile(cLon, len(lonflat)),
+                               np.tile(cLat, len(lonflat)),
+                               lonflat, latflat)
+    bearings = bflat.reshape(rows, cols)
+    distances = dflat.reshape(rows, cols) / 1000.
+    return np.radians(bearings), distances
 
 def bearing2theta(bearing):
     """
@@ -444,11 +200,12 @@ def bearing2theta(bearing):
               from east).
 
     """
-    theta = np.pi / 2. - bearing
-    theta = np.mod(theta, 2.*np.pi)
+    theta = np.pi / 2.0 - bearing
+    theta = np.mod(theta, 2.0 * np.pi)
 
     return theta
 
+
 def theta2bearing(theta):
     """
     Converts a cartesian angle (in radians) to an azimuthal bearing
@@ -461,13 +218,22 @@ def theta2bearing(theta):
                     from north).
 
     """
-    bearing = 2. * np.pi - (theta - np.pi / 2.)
-    bearing = np.mod(bearing, 2. * np.pi)
+    bearing = 2.0 * np.pi - (theta - np.pi / 2.0)
+    bearing = np.mod(bearing, 2.0 * np.pi)
 
     return bearing
 
-def makeGrid(cLon, cLat, margin=2, resolution=0.01, minLon=None, maxLon=None,
-             minLat=None, maxLat=None):
+
+def makeGrid(
+    cLon,
+    cLat,
+    margin=2,
+    resolution=0.01,
+    minLon=None,
+    maxLon=None,
+    minLat=None,
+    maxLat=None,
+):
     """
     Generate a grid of the distance and angle of a grid of points
     surrounding a storm centre given the location of the storm. The
@@ -487,8 +253,12 @@ def makeGrid(cLon, cLat, margin=2, resolution=0.01, minLon=None, maxLon=None,
     :returns: 2 2-d arrays containing the distance (km) and bearing (azimuthal)
               of all points in a grid from the ``cLon``, ``cLat``.
     """
-    if (type(cLon)==list or type(cLat)==list or
-        type(cLon)==np.ndarray or type(cLat)==np.ndarray):
+    if (
+        type(cLon) == list
+        or type(cLat) == list
+        or type(cLon) == np.ndarray
+        or type(cLat) == np.ndarray
+    ):
         raise TypeError("Input values must be scalar values")
 
     gridSize = int(resolution * 1000)
@@ -515,23 +285,26 @@ def makeGrid(cLon, cLat, margin=2, resolution=0.01, minLon=None, maxLon=None,
 
     try:
         from ._maputils import beardist
-        lonArray = xGrid / 1000.
-        latArray = yGrid / 1000.
-        R = np.zeros((len(latArray), len(lonArray)), order='F')
-        theta = np.zeros((len(latArray), len(lonArray)), order='F')
+
+        lonArray = xGrid / 1000.0
+        latArray = yGrid / 1000.0
+        R = np.zeros((len(latArray), len(lonArray)), order="F")
+        theta = np.zeros((len(latArray), len(lonArray)), order="F")
 
         beardist(cLon, cLat, lonArray, latArray, theta, R)
         R = np.ascontiguousarray(R)
         theta = np.ascontiguousarray(theta)
     except ImportError:
-        R = gridLatLonDist(cLon, cLat, xGrid / 1000., yGrid / 1000.)
-        theta = np.pi/2. - gridLatLonBear(cLon, cLat, xGrid / 1000., yGrid / 1000.)
+        bearing, R = gridLatLonDistBear(cLon, cLat, xGrid / 1000., yGrid / 1000.)
+        theta = np.pi / 2.0 - bearing
         np.putmask(R, R == 0, 1e-30)
 
     return R, theta
 
-def makeGridDomain(cLon, cLat, minLon, maxLon, minLat, maxLat,
-                   margin=2, resolution=0.01):
+
+def makeGridDomain(
+    cLon, cLat, minLon, maxLon, minLat, maxLat, margin=2, resolution=0.01
+):
     """
     Generate a grid of the distance and angle of a grid of points
     across a complete model domain, given the location of the storm.
@@ -551,8 +324,12 @@ def makeGridDomain(cLon, cLat, minLon, maxLon, minLat, maxLat,
               complete region.
 
     """
-    if (type(cLon)==list or type(cLat)==list or
-        type(cLon)==np.ndarray or type(cLat)==np.ndarray):
+    if (
+        type(cLon) == list
+        or type(cLat) == list
+        or type(cLon) == np.ndarray
+        or type(cLat) == np.ndarray
+    ):
         raise TypeError("Input values must be scalar values")
     gridSize = int(resolution * 1000)
     minLon_ = int(1000 * (minLon)) - int(1000 * margin)
@@ -563,12 +340,12 @@ def makeGridDomain(cLon, cLat, minLon, maxLon, minLat, maxLat,
     xGrid = np.array(np.arange(minLon_, maxLon_, gridSize), dtype=int)
     yGrid = np.array(np.arange(minLat_, maxLat_, gridSize), dtype=int)
 
-    R = gridLatLonDist(cLon, cLat, xGrid / 1000., yGrid / 1000.)
-    np.putmask(R, R==0, 1e-30)
-    theta = np.pi / 2. - gridLatLonBear(cLon, cLat,
-                                        xGrid / 1000., yGrid / 1000.)
+    bearing, R = gridLatLonDistBear(cLon, cLat, xGrid / 1000., yGrid / 1000.)
+    theta = np.pi / 2.0 - bearing
+    np.putmask(R, R == 0, 1e-30)
     return R, theta
 
+
 def meshLatLon(cLon, cLat, margin=2, resolution=0.01):
     """
     Create a meshgrid of the longitudes and latitudes of a grid.
@@ -584,8 +361,12 @@ def meshLatLon(cLon, cLat, margin=2, resolution=0.01):
               degrees of (``cLon``, ``cLat``).
 
     """
-    if (type(cLon)==list or type(cLat)==list or
-        type(cLon)==np.ndarray or type(cLat)==np.ndarray):
+    if (
+        type(cLon) == list
+        or type(cLat) == list
+        or type(cLon) == np.ndarray
+        or type(cLat) == np.ndarray
+    ):
         raise TypeError("Input values must be scalar values")
     gridSize = int(1000 * resolution)
 
@@ -598,10 +379,10 @@ def meshLatLon(cLon, cLat, margin=2, resolution=0.01):
     yy = np.array(np.arange(minLat, maxLat, gridSize))
 
     xGrid, yGrid = np.meshgrid(xx, yy)
-    return xGrid / 1000., yGrid / 1000.
+    return xGrid / 1000.0, yGrid / 1000.0
+
 
-def meshLatLonDomain(minLon, maxLon, minLat, maxLat,
-                     margin=2, resolution=0.01):
+def meshLatLonDomain(minLon, maxLon, minLat, maxLat, margin=2, resolution=0.01):
     """
     Create a meshgrid of the lon/lat grid across th full model domain.
 
@@ -628,7 +409,8 @@ def meshLatLonDomain(minLon, maxLon, minLat, maxLat,
     yy = np.array(np.arange(minLat_, maxLat_, gridSize))
 
     xGrid, yGrid = np.meshgrid(xx, yy)
-    return xGrid / 1000., yGrid / 1000.
+    return xGrid / 1000.0, yGrid / 1000.0
+
 
 def dist2GC(cLon1, cLat1, cLon2, cLat2, lonArray, latArray, units="km"):
     """
@@ -657,10 +439,10 @@ def dist2GC(cLon1, cLat1, cLon2, cLat2, lonArray, latArray, units="km"):
     """
 
     # Calculate distance and bearing from first point to array of points:
-    dist_ = gridLatLonDist(cLon1, cLat1, lonArray, latArray, units="rad")
-    bear_ = gridLatLonBear(cLon1, cLat1, lonArray, latArray)
+    bear_, dist_ = gridLatLonDistBear(cLon1, cLat1, lonArray, latArray)
+    dist_ *= 0.0001570783 # Convert to radians
 
-    #bearing of the cyclone:
+    # bearing of the cyclone:
     cyc_bear_ = latLon2Azi([cLon1, cLon2], [cLat1, cLat2])
 
     dist2GC_ = np.arcsin(np.sin(dist_) * np.sin(bear_ - cyc_bear_))
@@ -668,6 +450,7 @@ def dist2GC(cLon1, cLat1, cLon2, cLat2, lonArray, latArray, units="km"):
     distance = metutils.convert(dist2GC_, "rad", units)
     return distance
 
+
 def coriolis(lat):
     """
     Calculate the Coriolis factor (f) for a given latitude (degrees).
@@ -679,11 +462,12 @@ def coriolis(lat):
     :returns: Coriolis factor.
 
     """
-    omega = 2 * np.pi / 24. / 3600.
+    omega = 2 * np.pi / 24.0 / 3600.0
     f = 2 * omega * np.sin(np.radians(lat))
 
     return f
 
+
 def find_index(array, value):
     """
     Find the index of 'array' with a value closest to 'value'
@@ -707,10 +491,10 @@ def find_index(array, value):
     if type(value) == np.ndarray or type(value) == list:
         raise ValueError("Value cannot be an array")
 
-    if (value > array.max()):
+    if value > array.max():
         # Value is above the largest value in the array - return the last index:
         return len(array) - 1
-    elif (value < array.min()):
+    elif value < array.min():
         # Value is below minimum value in the array - return the first index:
         return 0
     else:
@@ -718,6 +502,7 @@ def find_index(array, value):
         idx = (abs(array - value)).argmin()
         return idx
 
+
 def find_nearest(array, value):
     """
     Find the closest value in 'array' to 'value'
diff --git a/tests/test_maputils.py b/tests/test_maputils.py
index 5092aa0f..f2ed2f46 100644
--- a/tests/test_maputils.py
+++ b/tests/test_maputils.py
@@ -66,13 +66,17 @@ class TestMapUtils(NumpyTestCase.NumpyTestCase):
                    110.2968, 110.0437, 109.7570, 109.4369, 109.0834, 108.6968, 108.2770,
                    107.8242, 107.3386, 106.8203, 106.2695, 105.6863, 105.0709, 104.4234])
 
-    AZI = np.array([135.0044, 135.0175, 135.0393, 135.0699, 135.1092, 135.1574,
-                    135.2143, 135.2802, 135.3549, 135.4385, 135.5312, 135.6329, 135.7437,
-                    135.8637, 135.9930, 136.1315, 136.2795, 136.4370, 136.6041, 136.7808])
-
-    Dist = np.array([157.1427, 157.1068, 157.0470, 156.9633, 156.8559, 156.7248,
-                     156.5700, 156.3918, 156.1902, 155.9654, 155.7176, 155.4470, 155.1538,
-                     154.8382, 154.5004, 154.1407, 153.7595, 153.3570, 152.9336, 152.4895])
+    AZI = np.array([134.81195977, 134.82952992, 134.85594569, 134.89121489,
+                    134.93535056, 134.98837098, 135.05029967, 135.12116539,
+                    135.20100215, 135.2898492 , 135.38775107, 135.4947575 ,
+                    135.61092353, 135.73630944, 135.87098074, 136.01500822,
+                    136.16846789, 136.33144099, 136.50401397, 136.68627847])
+
+    Dist = np.array([156.89956829, 156.8761494 , 156.82932912, 156.75914243,
+                     156.66564185, 156.54889743, 156.40899678, 156.24604514,
+                     156.0601654 , 155.85149814, 155.62020174, 155.36645239,
+                     155.09044421, 154.79238931, 154.4725179 , 154.13107836,
+                     153.76833737, 153.38457998, 152.98010978, 152.55524897])
 
     def test_ValuesXY2r(self):
         """Test xy2R function"""
@@ -102,13 +106,14 @@ def  test_GridLatLonDist(self):
         cLat = 1.0
         lonArray = np.array([0.59297447, 0.20873497, 0.44271653, 0.36579662, 0.06680392])
         latArray = np.array([0.5019297, 0.42174226, 0.23712093, 0.02745615, 0.13316245])
-        expected = ma.array([[56.30400807, 64.11584285, 55.7129098, 57.32175097, 73.35094702],
-                             [65.08411729, 71.94898913, 64.57343322, 65.9665517 , 80.28821237],
-                             [85.4022051, 90.74293694, 85.0136491, 86.07661618, 97.4877433],
-                             [108.56672733, 112.8162381 , 108.2613338, 109.09805046, 118.30941319],
-                             [96.87985127, 101.61920664, 96.537498, 97.47489149, 107.6850083]])
-
-        dist = maputils.gridLatLonDist(cLon, cLat, lonArray, latArray)
+        expected = np.array(
+            [[ 56.03776462,  63.90794904,  55.44173653,  57.06380581, 73.19967679],
+             [ 64.77263753,  71.6903351 ,  64.25767366,  65.66234269, 80.08427248],
+             [ 84.98739328,  90.3708789 ,  84.59556145,  85.66743417, 97.16440531],
+             [108.03529882, 112.31958766, 107.72732251, 108.57108803, 117.85478246],
+             [ 96.40715151, 101.18493545,  96.06190664,  97.00717907, 107.29601751]])
+
+        _, dist = maputils.gridLatLonDistBear(cLon, cLat, lonArray, latArray)
         self.numpyAssertAlmostEqual(dist, expected)
 
     def test_GridLatLonBear(self):
@@ -117,13 +122,14 @@ def test_GridLatLonBear(self):
         cLat = 1.0
         lonArray = np.array([0.59297447, 0.20873497, 0.44271653, 0.36579662, 0.06680392])
         latArray = np.array([0.5019297, 0.42174226, 0.23712093, 0.02745615, 0.13316245])
-        expected = np.array([[2.95705149, -2.61243611, -3.0270878 , -2.87840196, -2.42573357],
-                             [2.98217456, -2.67499093, -3.04285356, -2.91354889, -2.4986278 ],
-                             [3.02031491, -2.77686502, -3.06664317, -2.96745336, -2.62513214],
-                             [3.04627842, -2.85059019, -3.08275714, -3.0044598 , -2.72252399],
-                             [3.03474017, -2.81742223, -3.07560296, -2.9879869 , -2.67812704]])
-
-        bear = maputils.gridLatLonBear(cLon, cLat, lonArray, latArray)
+        expected = np.array(
+            [[ 2.95583639, -2.60950496, -3.02632306, -2.87671014, -2.42240667],
+             [ 2.9811187 , -2.6722875 , -3.04219261, -2.91206557, -2.4954023 ],
+             [ 3.01950581, -2.7746216 , -3.06614006, -2.96630381, -2.62224351],
+             [ 3.0456401 , -2.84873923, -3.08236162, -3.00354743, -2.72002218],
+             [ 3.03402572, -2.81538916, -3.07515962, -2.98696825, -2.67543611]])
+
+        bear, _ = maputils.gridLatLonDistBear(cLon, cLat, lonArray, latArray)
         self.numpyAssertAlmostEqual(bear, expected)
 
     def test_Bearing(self):