2022 / Day 14: Hexagons
Inadvertently I used up my hexagon idea for 2022 / Day 07: Raster. So instead we’ll look into another grid-creation issue I bumped into this summer.
Imagine You’d have to create a hexgrid over a large area in quite a large scale. PostGIS can generate this really nicely with st_hexagongrid. It consumes a requested grid size and the bounding geometry.
The only problem being that the Estonian territorial waters area that I needed it for is for one banana-shaped - from the Bay of Finland in the North-East through the Baltic Sea in the North-West, towards the Gulf of Riga down in the South-West. And if that wasn’t enough there’s islands, and some of them are of some size. Which means that doing a 10-meter grid for the sea would essentially mean generating millions upon millions of hexagons, only to dump at least 3/4 of them afterwards. Because they are on land :)
So, what if we used a tile-based approach instead. Meaning first do a more general grid with st_squaregrid get the grid cells that intersect the area of interest and then generate hexagons only into the those squares (and then dump the bits you don’t need).
The following query does just that - casts a half moon full of hexes without generating them everywhere.
with
minmax as (
select 1 as oid,
st_point(
40500.000000,5993000.000000,3301
) as ll,
st_point(
1064500.000000,7017000.000000,3301
) as ur
),
bounds as (
select
st_envelope(
st_collect(
array[ll, ur]
)
) as geom,
(st_x(ur)-st_x(ll)) as width,
st_x(ll) + (st_x(ur)-st_x(ll)) * (2.0 / 3.0) as r_x,
st_y(ll) + (st_y(ur)-st_y(ll)) * (1.0 / 3.0) as r_y,
st_x(ll) + (st_x(ur)-st_x(ll)) / 3.0 as l_x,
st_y(ll) + (st_y(ur)-st_y(ll)) * (2.0 / 3.0) as l_y
from minmax
),
moon as (
/* Represents the area of interest. I'll make it a half-moon-like
but can be anything. even islands e.g*/
select
m.geom, e.minx, e.maxx, focal_point, env
from (
select
width,
st_point(r_x, r_y, 3301) as focal_point,
st_difference(
st_buffer(st_point(l_x, l_y, 3301), width / 3.0),
st_buffer(st_point(r_x, r_y, 3301), width / 2.0)
) as geom
from
bounds
) m
join lateral
st_envelope(m.geom) env on true
join lateral (
select
min(st_x(d.geom)) as minx,
max(st_x(d.geom)) as maxx
from
st_dumppoints(env) d
) e on true
),
tiles as (
/* create smaller "tiles" we can use for generating a hexes in
reasonable places*/
select
(moon.maxx - moon.minx)/10.0 as tilewidth, g.*
from
moon
join lateral
st_squaregrid(
(moon.maxx - moon.minx)/10.0,
moon.env
) g on true
where
st_intersects(moon.geom, g.geom)
)
/* and create the hexgrid, tile by tile. Adding a something extra (st_scale)
for the image to not look just like a hexgrid. As if it was
measuring something :D*/
select
row_number() over()::int as oid, *
from (
select
distinct on (h.i, h.j)
st_scale(
h.geom,
st_point(
(min(l) over() / l),
(min(l) over() / l)
),
st_centroid(h.geom)
) as geom
from
moon,
tiles
join lateral
/* generate a hexagongrid within a specific tile using
tilewidth/20.0 as grid step*/
st_hexagongrid(
tiles.tilewidth/20.0,
tiles.geom
) h on true
join lateral
/* just fooling around to scale the hexes
afterwards*/
st_distance(
moon.focal_point,
st_centroid(h.geom)
) as l on true
where
st_intersects(moon.geom, h.geom)
) f
;