Creating a hexagonal grid of regular hexagons of definite area anywhere on the globe using shapely

I have a geojson file describing a polygon. I want to create an hexagonal grid on top of this polygon, with regular hexagons of area 90000 square meters.

Right now, I can either guarantee the area, or the regularity, but not both.

Here is my code:

def create_hexagon(l, x, y):
    """
    Create a hexagon centered on (x, y)
    :param l: length of the hexagon's edge
    :param x: x-coordinate of the hexagon's center
    :param y: y-coordinate of the hexagon's center
    :return: The polygon containing the hexagon's coordinates
    """
    c = [[x + math.cos(math.radians(angle)) * l, y + math.sin(math.radians(angle)) * l] for angle in range(0, 360, 60)]
    return Polygon(c)

def create_hexgrid(bbox, side):
    """
    returns an array of Points describing hexagons centers that are inside the given bounding_box
    :param bbox: The containing bounding box. The bbox coordinate should be in Webmercator.
    :param side: The size of the hexagons'
    :return: The hexagon grid
    """
    grid = []

    v_step = math.sqrt(3) * side
    h_step = 1.5 * side

    x_min = min(bbox[0], bbox[2])
    x_max = max(bbox[0], bbox[2])
    y_min = min(bbox[1], bbox[3])
    y_max = max(bbox[1], bbox[3])

    h_skip = math.ceil(x_min / h_step) - 1
    h_start = h_skip * h_step

    v_skip = math.ceil(y_min / v_step) - 1
    v_start = v_skip * v_step

    h_end = x_max + h_step
    v_end = y_max + v_step

    if v_start - (v_step / 2.0) < y_min:
        v_start_array = [v_start + (v_step / 2.0), v_start]
    else:
        v_start_array = [v_start - (v_step / 2.0), v_start]

    v_start_idx = int(abs(h_skip) % 2)

    c_x = h_start
    c_y = v_start_array[v_start_idx]
    v_start_idx = (v_start_idx + 1) % 2
    while c_x < h_end:
        while c_y < v_end:
            grid.append((c_x, c_y))
            c_y += v_step
        c_x += h_step
        c_y = v_start_array[v_start_idx]
        v_start_idx = (v_start_idx + 1) % 2

    return grid

I can either apply it on my polygon reprojected to webmercator (3857):

edge = math.sqrt(RESOLUTION**2/(3/2 * math.sqrt(3)))
hex_centers = create_hexgrid(reprojected.bounds, edge)
hexagons = GeometryCollection([
    shapely.ops.transform(webmercator_to_spherical, create_hexagon(edge, center[0], center[1]))
    for center in hex_centers 
    if any([zone.intersects(
        shapely.ops.transform(webmercator_to_spherical, create_hexagon(edge, center[0], center[1]))
    ) for zone in geometry.geoms])
])

and obtain a grid of regular hexagons, but the area would not be RESOLUTION**2

Or, I can apply it to my polygon reprojected to Albers Equal Area projection :

edge = math.sqrt(RESOLUTION**2/(3/2 * math.sqrt(3)))
hex_centers = create_hexgrid(reprojected_true.bounds, edge)
hexagons = GeometryCollection([
    reproject_from_true_meters(create_hexagon(edge, center[0], center[1]))
    for center in hex_centers 
    if any([zone.intersects(
        reproject_from_true_meters(create_hexagon(edge, center[0], center[1]))
    ) for zone in geometry.geoms])
])

To endup with a grid of hexagons with the right area, but not regular:

Here is the original polygon geojson:

{"type":"FeatureCollection","features":[{"type":"Feature","properties":{},"geometry":{"type":"Polygon","coordinates":[[[-64.30624008178711,-36.600369790618835],[-64.30709838867188,-36.60016307220288],[-64.30684089660645,-36.60250584848266],[-64.30615425109862,-36.60243694431346],[-64.30615425109862,-36.61366751126687],[-64.29619789123535,-36.61311635595829],[-64.29645538330078,-36.6160787694227],[-64.29774284362793,-36.61635433698211],[-64.2989444732666,-36.6169054691463],[-64.30117607116699,-36.61814550211169],[-64.30212020874023,-36.61897217967516],[-64.3033218383789,-36.62117660984155],[-64.30435180664062,-36.6229676629369],[-64.30512428283691,-36.62482755864398],[-64.30392265319824,-36.625171978850105],[-64.30632591247559,-36.62765175890191],[-64.28272247314453,-36.64762485464038],[-64.27800178527832,-36.64438818727331],[-64.27396774291992,-36.64108251425615],[-64.27250862121582,-36.637845571969294],[-64.27250862121582,-36.632542202379454],[-64.27216529846191,-36.62572304797872],[-64.26993370056152,-36.620212179400006],[-64.26976203918457,-36.606777786771396],[-64.27070617675781,-36.604297335278915],[-64.27302360534668,-36.60257475259031],[-64.27757263183594,-36.60147227948239],[-64.28186416625977,-36.60085213143531],[-64.28993225097656,-36.59954291363163],[-64.30624008178711,-36.600369790618835]]]}}]}

reprojected is the projection of this polygon in webmercator
reprojected_true is the projection of this polygon in Albers Equal Area projection

Is there just an other projection that I could use?