# How to create an 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, bbox)
x_max = max(bbox, bbox)
y_min = min(bbox, bbox)
y_max = max(bbox, bbox)

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, center))
for center in hex_centers
if any([zone.intersects(
shapely.ops.transform(webmercator_to_spherical, create_hexagon(edge, center, center))
) 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, center))
for center in hex_centers
if any([zone.intersects(
reproject_from_true_meters(create_hexagon(edge, center, center))
) 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? 