See this Gist for the full source.
A python script that generates hexagons from a Gosper hex tiling. It returns 7^ORDER hexagons, where ORDER is the Gosper fractal order (hardcoded to 7, yielding 823,543 hexagons). The script can be modified to generate more hexagons. If you want to accommodate for N hexagons, you can find the order you need by using this equation: ORDER = ceil(log(N, 7)).
Unlike the client-side L-system implementations (see for example, the Gosper curve or the Hilbert curve), this program does not load at once into memory a huge string representing the whole fractal, but instead computes the L-system rules while writing the results to the output file. The memory occupation is always O(ORDER*L), where L is the length of strings in the system rules.
In order to write to the file line by line, Python’s open() line buffering has been set up.
The result is saved in a CSV-compliant file (without header) named hexes.wkt.csv, in which each line represents a hexagon as a WKT POLYGON. This simple format offers the possibility for consumer scripts to load only the necessary amount of hexagons into memory, by reading them line by line. Therefore, the generator script can be run once to create a big file that can be used by different scripts, even if they do not need that many hexagons.