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]
        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

grid of regular hexagons

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:enter image description here

Here is the original polygon geojson:


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?


The is no solution to your problem, unfortunately.

On one hand, as soon as you project coordinates on a 2 dimensional map, you need to choose between an equal-area projection (preserving the area) or a conformal projection (preserving local angles), but no projection has both. Of course, some projections make a compromise so that neither are totally preserved in this case, buth distortions are limited. This is the case of the Miller projection, for instance. Using many local projection, if you can tolerate some discontinuities, will minimize the errors because the surface of the Earth is more and more similar to a plane when you zoom in (most can be optimized on your study area by changing the meridian and parallels of reference). EDIT: a practical solution could be to used the UTM projections, which are conformal but with relatively small area distortion. There are 60 zones of 6 degrees, so you will need to manage the discontinuities every 6 degrees. The good thing with UTM is that the distortion will be very little affected by the latitude.

enter image description here

On the other hand, it is not possible to cover a spherical object with only hexagonal faces (this is why a soccer ball is made of hexagons AND pentagons). If you look at this demonstration, no solid can be composed solely on hexagonal faces (so, even if the Earth is not a true sphere, it cannot be covered by hexagons only.

  • I'm aware that my problem isn't solvable if I want my grid to cover the whole surface of the sphere. But as you said, Earth is more and more similar to a plane when you zoom in. At the scale of a city, there should be okayish solutions. For example, the Albers Equal Area projection renders a pretty okay result over french cities. But an awful one over Argentinian cities. Considering the earth is a sphere, is there a way to "recenter" the Albers Equal Area projection over the center of my polygon for instance? – Borbag Nov 8 '19 at 13:42
  • 2
    @Borbag You are misunderstanding the problem. There are an infinite number of Albers Equal Area projections (multiple latitude and longitude variables). There is no one Albers that's applicable everywhere on the globe. And the Earth is a spheroid, not a sphere. – Vince Nov 8 '19 at 14:18
  • 1
    I wonder how they cut the hexagons and pentagons when they sew footballs and if it is possible to use just hexagons quora.com/…. – user30184 Nov 8 '19 at 14:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.