I've been thinking a lot about hexbin maps lately, after an earlier exploration of employment centers. In the previous map, I generated the hexbin with a QGIS plug-in and manually fixed up coastal hexes in a very non-reproducible way.

I since discovered that the CartoDB project have an open repository of postgresql functions which include functions for generating hexagonal grids.

This makes it easy to create hexagonal grids in postgresql/postgis and automate the process of creating land overlays for specific areas.

First of all I grabbed the LINZ coastlines and islands polygon layer, simplified it considerably, removed some of the smaller off-shore islands and loaded this into postgres as the table land_cover.

I then created a hexagonal grid with a side length of 10km that covers the entire extent of the coastline layer:

SELECT
    ST_SetSRID(
        CDB_HexagonGrid(
            (SELECT ST_Extent("geom") as extent from "land_cover"), 10000),
            2193
    ) as geom
INTO "hexgrid";

The next step is to remove the hexagons that don't cover land. The follow query finds hexagons that intersect with the land cover layer and contain at least 30% land. This factor is a somewhat arbitrary - I wanted it to be slightly optimistic about coverage:

SELECT a.*
INTO hexgrid_cover
FROM hexgrid a, land_cover b
WHERE
    ST_INTERSECTS(a.geom, b.geom) AND
    ST_AREA(ST_INTERSECTION(a.geom, b.geom)) / ST_AREA(a.geom) >= 0.3;

(Note that since the area of each hexagon is constant, the ST_AREA(a.geom) can be optimized out, but it's provided here for simplicity)

The result is a hexagonal grid coverage for the country. I've uploaded a few nationwide hexbin sets I've generated for the whole country as geojson on GitHub.

For this particular map, I've taken the NZ meshblock data with 'usually resident population' counts summed the totals for the hexagonal intersections weighted by area. I won't provide the detail here but:

1. Pre-calculate the area for each meshblock as a new column in the meshblock table.
2. Create a new layer of intersections between the meshblock layer and hexagons, including the hexagon ID, original meshblock area and population count.
3. Calculate the new (sliced) meshblock area, and divide this by the original area to get the proportion.
4. Multiply the proportional area by the original resident population count to get the weighted population.
5. Sum the weighted population over each hexagon id and assign that to the relevant hexagon.

Datasources

• NZ Coastlines and Islands Polygons (LINZ)
• Statistics NZ Pattern Files (Statistics NZ)
• Statistics NZ Tables (Statistics NZ)