diff --git a/libpysal/cg/alpha_shapes.py b/libpysal/cg/alpha_shapes.py
index 9f59993ff..e033ef187 100644
--- a/libpysal/cg/alpha_shapes.py
+++ b/libpysal/cg/alpha_shapes.py
@@ -15,21 +15,33 @@
 
 if not HAS_JIT:
     from warnings import warn
-    NUMBA_WARN = "Numba not imported, so alpha shape construction may be slower than expected."
+
+    NUMBA_WARN = (
+        "Numba not imported, so alpha shape construction may be slower than expected."
+    )
+
+try:
+    import pygeos
+
+    HAS_PYGEOS = True
+except ModuleNotFoundError:
+    HAS_PYGEOS = False
+
 
 EPS = np.finfo(float).eps
 
-__all__ = ['alpha_shape', 'alpha_shape_auto']
+__all__ = ["alpha_shape", "alpha_shape_auto"]
+
 
 @jit
 def nb_dist(x, y):
-    '''
+    """
     numba implementation of distance between points `x` and `y`
     ...
 
     Parameters
     ----------
-    
+
     x       : ndarray
               Coordinates of point `x`
     y       : ndarray
@@ -37,7 +49,7 @@ def nb_dist(x, y):
 
     Returns
     -------
-    
+
     dist    : float
               Distance between `x` and `y`
 
@@ -49,17 +61,18 @@ def nb_dist(x, y):
     >>> dist = nb_dist(x, y)
     >>> dist
     1.4142135623730951
-    
-    '''
+
+    """
     sum = 0
     for x_i, y_i in zip(x, y):
-        sum += (x_i - y_i)**2
+        sum += (x_i - y_i) ** 2
     dist = np.sqrt(sum)
     return dist
 
+
 @jit(nopython=True)
 def r_circumcircle_triangle_single(a, b, c):
-    '''
+    """
     Computation of the circumcircle of a single triangle
     ...
 
@@ -71,7 +84,7 @@ def r_circumcircle_triangle_single(a, b, c):
 
     Parameters
     ----------
-    
+
     a       : ndarray
               (2,) Array with coordinates of vertex `a` of the triangle
     b       : ndarray
@@ -81,7 +94,7 @@ def r_circumcircle_triangle_single(a, b, c):
 
     Returns
     -------
-    
+
     r       : float
               Circumcircle of the triangle
 
@@ -94,31 +107,29 @@ def r_circumcircle_triangle_single(a, b, c):
     >>> r = r_circumcircle_triangle_single(a, b, c)
     >>> r
     0.2500000000000001
-    
-    '''
+
+    """
     ab = nb_dist(a, b)
     bc = nb_dist(b, c)
     ca = nb_dist(c, a)
 
     num = ab * bc * ca
-    den = np.sqrt( (ab + bc + ca) * \
-                   (bc + ca - ab) * \
-                   (ca + ab - bc) * \
-                   (ab + bc - ca) )
+    den = np.sqrt((ab + bc + ca) * (bc + ca - ab) * (ca + ab - bc) * (ab + bc - ca))
     if den == 0:
         return np.array([ab, bc, ca]).max() / 2.0
     else:
         return num / den
 
+
 @jit(nopython=True)
 def r_circumcircle_triangle(a_s, b_s, c_s):
-    '''
+    """
     Computation of circumcircles for a series of triangles
     ...
 
     Parameters
     ----------
-    
+
     a_s     : ndarray
               (N, 2) array with coordinates of vertices `a` of the triangles
     b_s     : ndarray
@@ -128,7 +139,7 @@ def r_circumcircle_triangle(a_s, b_s, c_s):
 
     Returns
     -------
-    
+
     radii   : ndarray
               (N,) array with circumcircles for every triangle
 
@@ -141,37 +152,36 @@ def r_circumcircle_triangle(a_s, b_s, c_s):
     >>> rs = r_circumcircle_triangle(a_s, b_s, c_s)
     >>> rs
     array([3.53553391, 2.5       , 1.58113883])
-    '''
+    """
     len_a = len(a_s)
-    r2 = np.zeros( (len_a,) )
+    r2 = np.zeros((len_a,))
     for i in range(len_a):
-        r2[i] = r_circumcircle_triangle_single(a_s[i],
-                                               b_s[i],
-                                               c_s[i])
+        r2[i] = r_circumcircle_triangle_single(a_s[i], b_s[i], c_s[i])
     return r2
 
+
 @jit
 def get_faces(triangle):
-    '''
+    """
     Extract faces from a single triangle
     ...
 
     Parameters
     ----------
-    
+
     triangles       : ndarray
                       (3,) array with the vertex indices for a triangle
 
     Returns
     -------
-    
+
     faces           : ndarray
                       (3, 2) array with a row for each face containing the
                       indices of the two points that make up the face
 
     Examples
     --------
-    
+
     >>> triangle = np.array([3, 1, 4], dtype=np.int32)
     >>> faces = get_faces(triangle)
     >>> faces
@@ -179,23 +189,23 @@ def get_faces(triangle):
            [1., 4.],
            [4., 3.]])
 
-    '''
+    """
     faces = np.zeros((3, 2))
     for i, (i0, i1) in enumerate([(0, 1), (1, 2), (2, 0)]):
         faces[i] = triangle[i0], triangle[i1]
     return faces
 
+
 @jit
-def build_faces(faces, triangles_is,
-        num_triangles, num_faces_single):
-    '''
+def build_faces(faces, triangles_is, num_triangles, num_faces_single):
+    """
     Build facing triangles
 
     ...
 
     Parameters
     ----------
-    
+
     faces               : ndarray
                           (num_triangles * num_faces_single, 2) array of
                           zeroes in int form
@@ -210,14 +220,14 @@ def build_faces(faces, triangles_is,
 
     Returns
     -------
-    
+
     faces               : ndarray
                           Two dimensional array with a row for every facing
                           segment containing the indices of the coordinate points
 
     Examples
     --------
-    
+
     >>> import scipy.spatial as spat
     >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
     >>> triangulation = spat.Delaunay(pts)
@@ -242,23 +252,24 @@ def build_faces(faces, triangles_is,
            [1, 0],
            [0, 2]])
 
-    '''
+    """
     for i in range(num_triangles):
         from_i = num_faces_single * i
-        to_i = num_faces_single * (i+1)
-        faces[from_i: to_i] = get_faces(triangles_is[i])
+        to_i = num_faces_single * (i + 1)
+        faces[from_i:to_i] = get_faces(triangles_is[i])
     return faces
 
+
 @jit
 def nb_mask_faces(mask, faces):
-    '''
+    """
     Run over each row in `faces`, if the face in the following row is the
     same, then mark both as False on `mask`
     ...
 
     Parameters
     ----------
-    
+
     mask    : ndarray
               One-dimensional boolean array set to True with as many
               observations as rows in `faces`
@@ -268,13 +279,13 @@ def nb_mask_faces(mask, faces):
 
     Returns
     -------
-    
+
     masked  : ndarray
               Sequence of outward-facing faces
 
     Examples
     --------
-    
+
     >>> import numpy as np
     >>> faces = np.array([[0, 1], [0, 2], [1, 2], [1, 2], [1, 3], [1, 4], [1, 4], [2, 4], [3, 4]])
     >>> mask = np.ones((faces.shape[0], ), dtype=np.bool_)
@@ -285,35 +296,36 @@ def nb_mask_faces(mask, faces):
            [1, 3],
            [2, 4],
            [3, 4]])
-    
-    '''
-    for k in range(faces.shape[0]-1):
+
+    """
+    for k in range(faces.shape[0] - 1):
         if mask[k]:
-            if np.all(faces[k] == faces[k+1]):
+            if np.all(faces[k] == faces[k + 1]):
                 mask[k] = False
-                mask[k+1] = False
+                mask[k + 1] = False
     return faces[mask]
 
+
 def get_single_faces(triangles_is):
-    '''
+    """
     Extract outward facing edges from collection of triangles
     ...
 
     Parameters
     ----------
-    
+
     triangles_is    : ndarray
                       (D, 3) array, where D is the number of Delaunay triangles,
                       with the vertex indices for each triangle
 
     Returns
     -------
-    
+
     single_faces    : ndarray
 
     Example
     -------
-    
+
     >>> import scipy.spatial as spat
     >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
     >>> alpha = 0.33
@@ -329,15 +341,14 @@ def get_single_faces(triangles_is):
            [2, 4],
            [3, 4]])
 
-    '''
+    """
     num_faces_single = 3
     num_triangles = triangles_is.shape[0]
     num_faces = num_triangles * num_faces_single
     faces = np.zeros((num_faces, 2), dtype=np.int_)
     mask = np.ones((num_faces,), dtype=np.bool_)
 
-    faces = build_faces(faces, triangles_is,
-                        num_triangles, num_faces_single)
+    faces = build_faces(faces, triangles_is, num_triangles, num_faces_single)
 
     orderlist = ["x{}".format(i) for i in range(faces.shape[1])]
     dtype_list = [(el, faces.dtype.str) for el in orderlist]
@@ -349,16 +360,17 @@ def get_single_faces(triangles_is):
     single_faces = nb_mask_faces(mask, faces)
     return single_faces
 
-@requires('geopandas', 'shapely')
+
+@requires("geopandas", "shapely")
 def alpha_geoms(alpha, triangles, radii, xys):
-    '''
+    """
     Generate alpha-shape polygon(s) from `alpha` value, vertices of `triangles`,
     the `radii` for all points, and the points themselves
     ...
 
     Parameters
     ----------
-    
+
     alpha       : float
                   Alpha value to delineate the alpha-shape
     triangles   : ndarray
@@ -372,7 +384,7 @@ def alpha_geoms(alpha, triangles, radii, xys):
 
     Returns
     -------
-    
+
     geoms       : GeoSeries
                   Polygon(s) resulting from the alpha shape algorithm. The
                   GeoSeries object remains so even if only a single polygon is
@@ -380,7 +392,7 @@ def alpha_geoms(alpha, triangles, radii, xys):
 
     Examples
     --------
-    
+
     >>> import scipy.spatial as spat
     >>> pts = np.array([[0, 1], [3, 5], [4, 1], [6, 7], [9, 3]])
     >>> alpha = 0.33
@@ -406,29 +418,29 @@ def alpha_geoms(alpha, triangles, radii, xys):
     >>> geoms
     0    POLYGON ((0.00000 1.00000, 3.00000 5.00000, 4....
     dtype: geometry
-    
-    '''
+
+    """
     from shapely.geometry import LineString
     from shapely.ops import polygonize
     from geopandas import GeoSeries
 
-    triangles_reduced = triangles[radii < 1/alpha]
+    triangles_reduced = triangles[radii < 1 / alpha]
     outer_triangulation = get_single_faces(triangles_reduced)
     face_pts = xys[outer_triangulation]
-    geoms = GeoSeries(list(polygonize(list(map(LineString,
-                                               face_pts)))))
+    geoms = GeoSeries(list(polygonize(list(map(LineString, face_pts)))))
     return geoms
 
-@requires('geopandas', 'shapely')
+
+@requires("geopandas", "shapely")
 def alpha_shape(xys, alpha):
-    '''
+    """
     Alpha-shape delineation (Edelsbrunner, Kirkpatrick &
     Seidel, 1983) from a collection of points
     ...
 
     Parameters
     ----------
-    
+
     xys     : ndarray
               (N, 2) array with one point per row and coordinates structured as X
               and Y
@@ -437,7 +449,7 @@ def alpha_shape(xys, alpha):
 
     Returns
     -------
-    
+
     shapes  : GeoSeries
               Polygon(s) resulting from the alpha shape algorithm. The
               GeoSeries object remains so even if only a single polygon is
@@ -463,15 +475,14 @@ def alpha_shape(xys, alpha):
     Edelsbrunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of
         a set of points in the plane. IEEE Transactions on information theory,
         29(4), 551-559.
-    
-    '''
+
+    """
     if not HAS_JIT:
         warn(NUMBA_WARN)
     if xys.shape[0] < 4:
         from shapely import ops, geometry as geom
-        return ops.cascaded_union([geom.Point(xy)
-                                   for xy in xys])\
-                  .convex_hull.buffer(0)
+
+        return ops.cascaded_union([geom.Point(xy) for xy in xys]).convex_hull.buffer(0)
     triangulation = spat.Delaunay(xys)
     triangles = xys[triangulation.simplices]
     a_pts = triangles[:, 0, :]
@@ -482,39 +493,46 @@ def alpha_shape(xys, alpha):
     geoms = alpha_geoms(alpha, triangulation.simplices, radii, xys)
     return geoms
 
+
 def _valid_hull(geoms, points):
-    '''
+    """
     Sanity check within ``alpha_shape_auto()`` to verify the generated
     alpha shape actually contains the original set of points (xys).
-    
+
     Parameters
     ----------
-    
+
     geoms   : GeoSeries
               see alpha_geoms()
     points  : list
               xys parameter cast as shapely.geometry.Point objects
-    
+
     Returns
     -------
-    
+
     flag    : bool
               Valid hull for alpha shape [True] or not [False]
-    
-    '''
+
+    """
     flag = True
     # if there is not exactly one polygon
     if geoms.shape[0] != 1:
-        flag = False
+        return False
     # if any (xys) points do not intersect the polygon
-    for point in points:
-        if not point.intersects(geoms[0]):
-            flag = False
-    return flag
-
-@requires('geopandas', 'shapely')
-def alpha_shape_auto(xys, step=1, verbose=False):
-    '''
+    if HAS_PYGEOS:
+        return pygeos.intersects(pygeos.from_shapely(geoms[0]), points).all()
+    else:
+        for point in points:
+            if not point.intersects(geoms[0]):
+                return False
+        return True
+
+
+@requires("geopandas", "shapely")
+def alpha_shape_auto(
+    xys, step=1, verbose=False, return_radius=False, return_circles=False
+):
+    """
     Computation of alpha-shape delineation with automated selection of alpha.
     ...
 
@@ -529,7 +547,7 @@ def alpha_shape_auto(xys, step=1, verbose=False):
 
     Parameters
     ----------
-    
+
     xys     : ndarray
               Nx2 array with one point per row and coordinates structured as X
               and Y
@@ -564,60 +582,155 @@ def alpha_shape_auto(xys, step=1, verbose=False):
     Edelsbrunner, H., Kirkpatrick, D., & Seidel, R. (1983). On the shape of
         a set of points in the plane. IEEE Transactions on information theory,
         29(4), 551-559.
-    
-    '''
+
+    """
     if not HAS_JIT:
         warn(NUMBA_WARN)
     from shapely import geometry as geom
+
+    if return_circles:
+        return_radius = True
     if xys.shape[0] < 4:
         from shapely import ops
-        return ops.cascaded_union([geom.Point(xy)
-                                   for xy in xys])\
-                  .convex_hull.buffer(0)
+
+        if xys.shape[0] == 3:
+            multipoint = ops.cascaded_union([geom.Point(xy) for xy in xys])
+            alpha_shape = multipoint.convex_hull.buffer(0)
+        else:
+            alpha_shape = geom.Polygon([])
+        if xys.shape[0] == 1:
+            if return_radius:
+                if return_circles:
+                    out = [alpha_shape, 0, alpha_shape]
+                return alpha_shape, 0
+            return alpha_shape
+        elif xys.shape[0] == 2:
+            if return_radius:
+                r = spat.distance.euclidean(xys[0], xys[1]) / 2
+                if return_circles:
+                    circle = _construct_centers(xys[0], xys[1], r)
+                    return [alpha_shape, r, circle]
+                return [alpha_shape, r]
+            return alpha_shape
+        elif return_radius:  # this handles xys.shape[0] == 3
+            radius = r_circumcircle_triangle_single(xys[0], xys[1], xys[2])
+            if return_circles:
+                circles = construct_bounding_circles(alpha_shape, radius)
+                return [alpha_shape, radius, circles]
+            return [alpha_shape, radius]
+        return alpha_shape
     triangulation = spat.Delaunay(xys)
     triangles = xys[triangulation.simplices]
     a_pts = triangles[:, 0, :]
     b_pts = triangles[:, 1, :]
     c_pts = triangles[:, 2, :]
     radii = r_circumcircle_triangle(a_pts, b_pts, c_pts)
-    radii[np.isnan(radii)] = 0 # "Line" triangles to be kept for sure
+    radii[np.isnan(radii)] = 0  # "Line" triangles to be kept for sure
     del triangles, a_pts, b_pts, c_pts
     radii_sorted_i = radii.argsort()
     triangles = triangulation.simplices[radii_sorted_i][::-1]
     radii = radii[radii_sorted_i][::-1]
-    geoms_prev = alpha_geoms((1/radii.max())-EPS, triangles, radii, xys)
-    xys_bb = np.array([*xys.min(axis=0), *xys.max(axis=0)])
-    points = [geom.Point(pnt) for pnt in xys]
+    geoms_prev = alpha_geoms((1 / radii.max()) - EPS, triangles, radii, xys)
+    if HAS_PYGEOS:
+        points = pygeos.points(xys)
+    else:
+        points = [geom.Point(pnt) for pnt in xys]
     if verbose:
-        print('Step set to %i'%step)
+        print("Step set to %i" % step)
     for i in range(0, len(radii), step):
         radi = radii[i]
         alpha = (1 / radi) - EPS
         if verbose:
-            print('%.2f%% | Trying a = %f'\
-            %((i+1)/radii.shape[0], alpha))
+            print("%.2f%% | Trying a = %f" % ((i + 1) / radii.shape[0], alpha))
         geoms = alpha_geoms(alpha, triangles, radii, xys)
         if _valid_hull(geoms, points):
             geoms_prev = geoms
+            radi_prev = radi
         else:
             break
     if verbose:
         print(geoms_prev.shape)
-    return geoms_prev[0] # Return a shapely polygon
+    if return_radius:
+        out = [geoms_prev[0], radi_prev]
+        if return_circles:
+            out.append(construct_bounding_circles(out[0], radi_prev))
+        return out
+    # Return a shapely polygon
+    return geoms_prev[0]
+
+
+def construct_bounding_circles(alpha_shape, radius):
+    """
+    Construct the bounding circles for an alpha shape, given the radius 
+    computed from the `alpha_shape_auto` method. 
+
+    Arguments
+    ---------
+    alpha_shape : shapely.Polygon
+                  an alpha-hull with the input radius. 
+    radius      : float
+                  the radius of the input alpha_shape.
+
+    Returns
+    -------
+    numpy.ndarray of shape (n,2) containing the centers of the circles
+    defining the alpha_shape. 
+
+    """
+    coordinates = list(alpha_shape.boundary.coords)
+    n_coordinates = len(coordinates)
+    centers = []
+    for i in range(n_coordinates - 1):
+        a, b = coordinates[i], coordinates[i + 1]
+        centers.append(_construct_centers(a, b, radius))
+    return centers
 
-if __name__ == '__main__':
+
+@jit(nopython=True)
+def _construct_centers(a, b, radius):
+    midpoint_x = (a[0] + b[0]) * 0.5
+    midpoint_y = (a[1] + b[1]) * 0.5
+    d = ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** 0.5
+    if b[0] - a[0] == 0:
+        m = np.inf
+        axis_rotation = np.pi / 2
+    else:
+        m = (b[1] - a[1]) / (b[0] - a[0])
+        axis_rotation = np.arctan(m)
+    # altitude is perpendicular bisector of AB
+    interior_angle = np.arccos(0.5 * d / radius)
+    chord = np.sin(interior_angle) * radius
+
+    dx = chord * np.sin(axis_rotation)
+    dy = chord * np.cos(axis_rotation)
+
+    up_x = midpoint_x - dx
+    up_y = midpoint_y + dy
+    down_x = midpoint_x + dx
+    down_y = midpoint_y - dy
+
+    # sign gives us direction of point, since
+    # shapely shapes are clockwise-defined
+    sign = np.sign((b[0] - a[0]) * (up_y - a[1]) - (b[1] - a[1]) * (up_x - a[0]))
+    if sign == 1:
+        return up_x, up_y
+    else:
+        return down_x, down_y
+
+
+if __name__ == "__main__":
 
     import matplotlib.pyplot as plt
     import time
     import geopandas as gpd
-    plt.close('all')
+
+    plt.close("all")
     xys = np.random.random((1000, 2))
     t0 = time.time()
     geoms = alpha_shape_auto(xys, 1)
     t1 = time.time()
-    print('%.2f Seconds to run algorithm'%(t1-t0))
+    print("%.2f Seconds to run algorithm" % (t1 - t0))
     f, ax = plt.subplots(1)
-    gpd.GeoDataFrame({'geometry':[geoms]}).plot(ax=ax, color='orange', alpha=0.5)
+    gpd.GeoDataFrame({"geometry": [geoms]}).plot(ax=ax, color="orange", alpha=0.5)
     ax.scatter(xys[:, 0], xys[:, 1], s=0.1)
     plt.show()
-
diff --git a/libpysal/cg/tests/data/alpha_05.cpg b/libpysal/cg/tests/data/alpha_05.cpg
deleted file mode 100644
index cd89cb975..000000000
--- a/libpysal/cg/tests/data/alpha_05.cpg
+++ /dev/null
@@ -1 +0,0 @@
-ISO-8859-1
\ No newline at end of file
diff --git a/libpysal/cg/tests/data/alpha_05.dbf b/libpysal/cg/tests/data/alpha_05.dbf
deleted file mode 100644
index 38725f2db..000000000
Binary files a/libpysal/cg/tests/data/alpha_05.dbf and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_05.gpkg b/libpysal/cg/tests/data/alpha_05.gpkg
new file mode 100644
index 000000000..dc4c53609
Binary files /dev/null and b/libpysal/cg/tests/data/alpha_05.gpkg differ
diff --git a/libpysal/cg/tests/data/alpha_05.shp b/libpysal/cg/tests/data/alpha_05.shp
deleted file mode 100644
index 72d619ba5..000000000
Binary files a/libpysal/cg/tests/data/alpha_05.shp and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_05.shx b/libpysal/cg/tests/data/alpha_05.shx
deleted file mode 100644
index ebe0ea695..000000000
Binary files a/libpysal/cg/tests/data/alpha_05.shx and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_auto.cpg b/libpysal/cg/tests/data/alpha_auto.cpg
deleted file mode 100644
index cd89cb975..000000000
--- a/libpysal/cg/tests/data/alpha_auto.cpg
+++ /dev/null
@@ -1 +0,0 @@
-ISO-8859-1
\ No newline at end of file
diff --git a/libpysal/cg/tests/data/alpha_auto.dbf b/libpysal/cg/tests/data/alpha_auto.dbf
deleted file mode 100644
index 088c3d1b5..000000000
Binary files a/libpysal/cg/tests/data/alpha_auto.dbf and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_auto.gpkg b/libpysal/cg/tests/data/alpha_auto.gpkg
new file mode 100644
index 000000000..c39c1dcd5
Binary files /dev/null and b/libpysal/cg/tests/data/alpha_auto.gpkg differ
diff --git a/libpysal/cg/tests/data/alpha_auto.shp b/libpysal/cg/tests/data/alpha_auto.shp
deleted file mode 100644
index 929f81edc..000000000
Binary files a/libpysal/cg/tests/data/alpha_auto.shp and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_auto.shx b/libpysal/cg/tests/data/alpha_auto.shx
deleted file mode 100644
index ba90ee2b8..000000000
Binary files a/libpysal/cg/tests/data/alpha_auto.shx and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_fifth.cpg b/libpysal/cg/tests/data/alpha_fifth.cpg
deleted file mode 100644
index cd89cb975..000000000
--- a/libpysal/cg/tests/data/alpha_fifth.cpg
+++ /dev/null
@@ -1 +0,0 @@
-ISO-8859-1
\ No newline at end of file
diff --git a/libpysal/cg/tests/data/alpha_fifth.dbf b/libpysal/cg/tests/data/alpha_fifth.dbf
deleted file mode 100644
index 38725f2db..000000000
Binary files a/libpysal/cg/tests/data/alpha_fifth.dbf and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_fifth.gpkg b/libpysal/cg/tests/data/alpha_fifth.gpkg
new file mode 100644
index 000000000..b045b9229
Binary files /dev/null and b/libpysal/cg/tests/data/alpha_fifth.gpkg differ
diff --git a/libpysal/cg/tests/data/alpha_fifth.shp b/libpysal/cg/tests/data/alpha_fifth.shp
deleted file mode 100644
index 5467d404e..000000000
Binary files a/libpysal/cg/tests/data/alpha_fifth.shp and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_fifth.shx b/libpysal/cg/tests/data/alpha_fifth.shx
deleted file mode 100644
index 1fe56681c..000000000
Binary files a/libpysal/cg/tests/data/alpha_fifth.shx and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_fourth.cpg b/libpysal/cg/tests/data/alpha_fourth.cpg
deleted file mode 100644
index cd89cb975..000000000
--- a/libpysal/cg/tests/data/alpha_fourth.cpg
+++ /dev/null
@@ -1 +0,0 @@
-ISO-8859-1
\ No newline at end of file
diff --git a/libpysal/cg/tests/data/alpha_fourth.dbf b/libpysal/cg/tests/data/alpha_fourth.dbf
deleted file mode 100644
index 38725f2db..000000000
Binary files a/libpysal/cg/tests/data/alpha_fourth.dbf and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_fourth.gpkg b/libpysal/cg/tests/data/alpha_fourth.gpkg
new file mode 100644
index 000000000..a9edb5b86
Binary files /dev/null and b/libpysal/cg/tests/data/alpha_fourth.gpkg differ
diff --git a/libpysal/cg/tests/data/alpha_fourth.shp b/libpysal/cg/tests/data/alpha_fourth.shp
deleted file mode 100644
index c25b9ad5a..000000000
Binary files a/libpysal/cg/tests/data/alpha_fourth.shp and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_fourth.shx b/libpysal/cg/tests/data/alpha_fourth.shx
deleted file mode 100644
index 2f27b1dc2..000000000
Binary files a/libpysal/cg/tests/data/alpha_fourth.shx and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_tenth.cpg b/libpysal/cg/tests/data/alpha_tenth.cpg
deleted file mode 100644
index cd89cb975..000000000
--- a/libpysal/cg/tests/data/alpha_tenth.cpg
+++ /dev/null
@@ -1 +0,0 @@
-ISO-8859-1
\ No newline at end of file
diff --git a/libpysal/cg/tests/data/alpha_tenth.dbf b/libpysal/cg/tests/data/alpha_tenth.dbf
deleted file mode 100644
index 38725f2db..000000000
Binary files a/libpysal/cg/tests/data/alpha_tenth.dbf and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_tenth.gpkg b/libpysal/cg/tests/data/alpha_tenth.gpkg
new file mode 100644
index 000000000..1f5b51b9e
Binary files /dev/null and b/libpysal/cg/tests/data/alpha_tenth.gpkg differ
diff --git a/libpysal/cg/tests/data/alpha_tenth.shp b/libpysal/cg/tests/data/alpha_tenth.shp
deleted file mode 100644
index f212eabee..000000000
Binary files a/libpysal/cg/tests/data/alpha_tenth.shp and /dev/null differ
diff --git a/libpysal/cg/tests/data/alpha_tenth.shx b/libpysal/cg/tests/data/alpha_tenth.shx
deleted file mode 100644
index 5ef429313..000000000
Binary files a/libpysal/cg/tests/data/alpha_tenth.shx and /dev/null differ
diff --git a/libpysal/cg/tests/data/eberly_bounding_circles.gpkg b/libpysal/cg/tests/data/eberly_bounding_circles.gpkg
new file mode 100644
index 000000000..df2bacec2
Binary files /dev/null and b/libpysal/cg/tests/data/eberly_bounding_circles.gpkg differ
diff --git a/libpysal/cg/tests/test_ashapes.py b/libpysal/cg/tests/test_ashapes.py
index b4f36f45d..c366268dd 100644
--- a/libpysal/cg/tests/test_ashapes.py
+++ b/libpysal/cg/tests/test_ashapes.py
@@ -7,6 +7,7 @@
 try:
     import geopandas
     from shapely import geometry
+
     GEOPANDAS_EXTINCT = False
 except ImportError:
     GEOPANDAS_EXTINCT = True
@@ -14,26 +15,58 @@
 this_directory = os.path.dirname(__file__)
 
 
-@skipIf(GEOPANDAS_EXTINCT, 'Geopandas is missing, so test will not run')
+@skipIf(GEOPANDAS_EXTINCT, "Geopandas is missing, so test will not run")
 class Test_Alpha_Shapes(TestCase):
     def setUp(self):
-        eberly = geopandas.read_file(get_path('eberly_net.shp'))
-        eberly_vertices = eberly.geometry.apply(lambda x: np.hstack(x.xy).reshape(2, 2).T).values
+        eberly = geopandas.read_file(get_path("eberly_net.shp"))
+        eberly_vertices = eberly.geometry.apply(
+            lambda x: np.hstack(x.xy).reshape(2, 2).T
+        ).values
         eberly_vertices = np.vstack(eberly_vertices)
         self.vertices = eberly_vertices
 
-        self.a05 = geopandas.read_file(os.path.join(this_directory, 'data/alpha_05.shp')).geometry.to_numpy().item()
-        self.a10 = geopandas.read_file(os.path.join(this_directory, 'data/alpha_tenth.shp')).geometry.to_numpy().item()
-        self.a2  = geopandas.read_file(os.path.join(this_directory, 'data/alpha_fifth.shp')).geometry.to_numpy().item()
-        self.a25 = geopandas.read_file(os.path.join(this_directory, 'data/alpha_fourth.shp')).geometry.to_numpy().item()
-        self.a25 = geopandas.read_file(os.path.join(this_directory, 'data/alpha_fourth.shp')).geometry.to_numpy().item()
+        self.a05 = (
+            geopandas.read_file(os.path.join(this_directory, "data/alpha_05.gpkg"))
+            .geometry.to_numpy()
+            .item()
+        )
+        self.a10 = (
+            geopandas.read_file(os.path.join(this_directory, "data/alpha_tenth.gpkg"))
+            .geometry.to_numpy()
+            .item()
+        )
+        self.a2 = (
+            geopandas.read_file(os.path.join(this_directory, "data/alpha_fifth.gpkg"))
+            .geometry.to_numpy()
+            .item()
+        )
+        self.a25 = (
+            geopandas.read_file(os.path.join(this_directory, "data/alpha_fourth.gpkg"))
+            .geometry.to_numpy()
+            .item()
+        )
+        self.a25 = (
+            geopandas.read_file(os.path.join(this_directory, "data/alpha_fourth.gpkg"))
+            .geometry.to_numpy()
+            .item()
+        )
+        circles = geopandas.read_file(
+            os.path.join(this_directory, "data/eberly_bounding_circles.gpkg")
+        )
+        self.circle_radii = circles.radius.iloc[0]
+        self.circle_verts = np.column_stack(
+            (circles.geometry.x.values, circles.geometry.y.values)
+        )
+
+        self.autoalpha = geopandas.read_file(
+            os.path.join(this_directory, "data/alpha_auto.gpkg")
+        ).geometry[0]
 
-        self.autoalpha = geopandas.read_file(os.path.join(this_directory, 'data/alpha_auto.shp')).geometry[0]
     def test_alpha_shapes(self):
-        new_a05 = alpha_shape(self.vertices, .05).to_numpy().item()
-        new_a10 = alpha_shape(self.vertices, .10).to_numpy().item()
-        new_a2  = alpha_shape(self.vertices, .2).to_numpy().item()
-        new_a25 = alpha_shape(self.vertices, .25).to_numpy().item()
+        new_a05 = alpha_shape(self.vertices, 0.05).to_numpy().item()
+        new_a10 = alpha_shape(self.vertices, 0.10).to_numpy().item()
+        new_a2 = alpha_shape(self.vertices, 0.2).to_numpy().item()
+        new_a25 = alpha_shape(self.vertices, 0.25).to_numpy().item()
 
         assert new_a05.equals(self.a05)
         assert new_a10.equals(self.a10)
@@ -46,9 +79,22 @@ def test_auto(self):
         assert self.autoalpha.equals(auto_alpha)
 
     def test_small_n(self):
-        new_singleton = alpha_shape(self.vertices[0].reshape(1,-1), .5)
+        new_singleton = alpha_shape(self.vertices[0].reshape(1, -1), 0.5)
         assert isinstance(new_singleton, geometry.Polygon)
-        new_duo = alpha_shape(self.vertices[:1], .5)
+        new_duo = alpha_shape(self.vertices[:1], 0.5)
         assert isinstance(new_duo, geometry.Polygon)
-        new_triple = alpha_shape(self.vertices[:2], .5)
+        new_triple = alpha_shape(self.vertices[:2], 0.5)
         assert isinstance(new_triple, geometry.Polygon)
+        new_triple = alpha_shape_auto(
+            self.vertices[0].reshape(1, -1), return_circles=True
+        )
+        assert isinstance(new_triple[0], geometry.Polygon)
+        new_triple = alpha_shape_auto(self.vertices[:1], return_circles=True)
+        assert isinstance(new_triple[0], geometry.Polygon)
+        new_triple = alpha_shape_auto(self.vertices[:2], return_circles=True)
+        assert isinstance(new_triple[0], geometry.Polygon)
+
+    def test_circles(self):
+        ashape, radius, centers = alpha_shape_auto(self.vertices, return_circles=True)
+        np.testing.assert_allclose(radius, self.circle_radii)
+        np.testing.assert_allclose(centers, self.circle_verts)