[GIS] Create Hexbin (honeycomb) grid using command line or Python

Question

Qgis has this amazing Create Grid function. Is there an alternative for python / unix command line?

I've found Gdal_grid but it is not the same. The goal is to eventually create a hexbin heatmap.

Answer 1

Use GDAL with the SQLite/Spatialite dialect. Then you can utilize the ST_HexagonalGrid function that is documented in http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html.

return a grid of hexagonal cells (having the edge length of size) precisely covering the input Geometry. The returned Geometry will usually be of the MultiPolygon type (a collection of Squares), but will be a MultiLinestring if the optional edges_only argument is set to TRUE If the optional origin argument (expected to be a Point) is not specified then the (0,0) grid origin will be assumed by default.

The result is naturally collection of hexagons, not squares. Obviously manual text is copied from ST_SquareGrid function.

Usage example with ogrinfo follows. For ogrinfo there must be some bogus datasource, I used a shapefile "temp.shp".

ogrinfo -dialect SQLite -sql "SELECT ST_HexagonalGrid(ST_GeomFromText('POLYGON (( 0 47, 0 50, 3 50, 3 47, 0 47 ))'),1) from test LIMIT 1" test.shp

The result as WKT:

  MULTIPOLYGON (((0.0 46.7653718043597,0.5 45.8993464005753,1.5 45.8993464005753,2.0 46.7653718043597,1.5 47.6313972081442,0.5 47.6313972081442,0.0 46.7653718043597)),((-1.5 47.6313972081442,-1.0 46.7653718043597,0.0 46.7653718043597,0.5 47.6313972081442,0.0 48.4974226119286,-1.0 48.4974226119286,-1.5 47.6313972081442)),((1.5 47.6313972081442,2.0 46.7653718043597,3.0 46.7653718043597,3.5 47.6313972081442,3.0 48.4974226119286,2.0 48.4974226119286,1.5 47.6313972081442)),((0.0 48.4974226119286,0.5 47.6313972081442,1.5 47.6313972081442,2.0 48.4974226119286,1.5 49.363448015713,0.5 49.363448015713,0.0 48.4974226119286)),((-1.5 49.363448015713,-1.0 48.4974226119286,0.0 48.4974226119286,0.5 49.363448015713,0.0 50.2294734194975,-1.0 50.2294734194975,-1.5 49.363448015713)),((1.5 49.363448015713,2.0 48.4974226119286,3.0 48.4974226119286,3.5 49.363448015713,3.0 50.2294734194975,2.0 50.2294734194975,1.5 49.363448015713)),((0.0 50.2294734194975,0.5 49.363448015713,1.5 49.363448015713,2.0 50.2294734194975,1.5 51.0954988232819,0.5 51.0954988232819,0.0 50.2294734194975)))

And this is how it looks. The input geometry that was to be covered by the grid is shown as a yellow rectangle.