# Hexagons with Python

Blog Post created by jmward on May 8, 2015

I was reading Dan Patterson's blog about his Phish_Nyet sample from his Numpy Code Snippets series and thought I would give it a whirl.  My first time through I entered too many rows and columns and it ran out of memory.  So after scaling it back, I got the tool to run.  I wanted to see if I could get the script to create hexagons instead of squares or rectangles.

Here is Dan's original code:

"""
Phish_Nyet.py
Author:  Dan.Patterson@carleton.ca
Purpose: Produce a sampling grid with user defined parameters.
"""
import arcpy
import numpy as np

def demo():
"""Generate the grid using the following parameter"""
output_shp = r'C:\temp\Phish_Nyet.shp'
SR =  arcpy.SpatialReference(2951) # u'NAD_1983_CSRS_MTM_9' YOU NEED ONE!!!
corner = [340000.0, 5022000.0]     # corner of grid
dX = 1000.0;  dY = 1000.0          # X and Y cell widths
cols = 3;  rows = 3                # columns/rows...grids in X and Y direction
angle = 0                          # rotation angle, clockwise +ve
# create the grid
pnts = grid_array(corner,dX,dY,cols,rows,angle)
output_polygons(output_shp,SR,pnts)
print('\nPhish_Nyet has created... {}'.format(output_shp))

def grid_array(corner=[0,0],dX=1,dY=1,cols=1,rows=1,angle=0):
"""create the array of pnts to pass on to arcpy using numpy magic"""
X = [0.0,0.0,1.0,1.0,0.0]                  # X,Y values for a unit square
Y = [0.0,1.0,1.0,0.0,0.0]                  #
seed = np.column_stack((X,Y)) * [dX,dY]    # first array corner values scaled
u = [seed + [j*dX,i*dY] for i in range(0,rows) for j in range(0,cols)]
pnts = np.array(u)                         #
x = [rotate(p,angle) for p in pnts]        # rotate the scaled points
pnts = [ p + corner for p in x]            # translate them
return pnts

def rotate(pnts,angle=0):
"""rotate points about the origin in degrees, (+ve for clockwise) """
s = np.sin(angle);  c = np.cos(angle)   # rotation terms
aff_matrix = np.array([[c, -s],[s, c]]) # rotation matrix
XY_r = np.dot(pnts, aff_matrix)         # numpy magic to rotate pnts
return XY_r

def output_polygons(output_shp,SR,pnts):
"""produce the output polygon shapefile"""
assert (SR != None) and (SR.type=='Projected'), msg
polygons = []
for pnt in pnts:                           # create the polygon geometry
polygons.append(arcpy.Polygon(arcpy.Array([arcpy.Point(*xy) for xy in pnt]),SR))
if arcpy.Exists(output_shp):               # overwrite any existing versions
arcpy.Delete_management(output_shp)
arcpy.CopyFeatures_management(polygons, output_shp)

if __name__ == '__main__':
"""Generate the grid using the listed parameters"""
demo()  # modify the parameters in demo to run


To modify the shape of the feature to be hexagons you need to modify these two lists:

X = [0.0,0.0,1.0,1.0,0.0]
Y = [0.0,1.0,1.0,0.0,0.0]


To be:

X = [0.0, 0.5, 1.5, 2.0, 1.5, 0.5, 0.0]
Y = [0.0, 0.866, 0.866, 0.0, -0.866, -0.866, 0.0]


Line 26 in Dan's original script above creates a two dimensional array.  Essentially it turns each list on its end and makes a column out of it, so in this case X gets upended next to Y and the two together end up forming coordinate pairs.  It is then multiplied by the side lengths.

The next line (27) uses list comprehension to create a three dimensional array.  I had to modify this line and add another like it to get my hexagons to stack.  With no modifications the hexagons overlap in a fish scale type pattern, which is also cool but I wanted them stacked.  My modifications are these:

u = [seed + [j*(dX*3.0),i*(dY*0.866)] for i in range(0,rows) for j in range(0,cols) if i % 2 == 0]
t = [seed + [(j*(dX*3.0))+1.5*dX, i*(dY*0.866)] for i in range(0,rows) for j in range(0,cols) if not (i % 2 == 0)]
u = u + t


I had to create two of these arrays to overcome the overlapping hexagons and get them to stack properly.  The mod division at the end of each line is like a bit switch that determines whether to offset the hexagons on that row.  So I am creating two arrays that contain coordinates for polygons on every other row.  One row looks like this:

Two rows like:

Now I'm thinking of all sorts of things to do with hexagons.  Thanks Dan for helping me hone my python chops and introducing me to list comprehension.