I'm putting this here for posterity. This is how pretty accurate centroids can be determined from a set of points. This was originally written using the shapefile library, and I no longer need it in my code because I'm switching everything to the geopandas package (because of another step that needs geoprocessing that the shapefile library doesn…

get_centroid.py

    def getCentroidFromPoints(self, s):
        points = s.points
        p_x = 0
        p_y = 0
        pp = np.array(s.points)
        pp = np.append(pp, [pp[0]], axis = 0)
        area = self.shoelace(pp)
        for p in range(0, len(pp) - 1):
            p_x += (pp[p, 0] + pp[p + 1, 0])*(pp[p, 0] * pp[p + 1, 1] - pp[p + 1, 0] * pp[p, 1]) 
            p_y += (pp[p, 1] + pp[p + 1, 1])*(pp[p, 0] * pp[p + 1, 1] - pp[p + 1, 0] * pp[p, 1]) 
        return((-1.0 / (6 * area) * p_x, -1.0 / (6 * area) * p_y))

    def shoelace(self, x_y):
        # Thanks to https://stackoverflow.com/a/58515054/877387
        x_y = x_y.reshape(-1,2)
        x = x_y[:,0]
        y = x_y[:,1]
        S1 = np.sum(x * np.roll(y,-1))
        S2 = np.sum(y * np.roll(x,-1))
        area = .5 * np.absolute(S1 - S2)
        return area