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:


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)  
  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) """  
  angle = np.deg2rad(angle)               # convert to radians  
  s = np.sin(angle);  c = np.cos(angle)   # rotation terms  
  aff_matrix = np.array([[c, -s],[s, c]]) # rotation matrix  
  XY_r =, aff_matrix)         # numpy magic to rotate pnts  
  return XY_r  
def output_polygons(output_shp,SR,pnts):  
  """produce the output polygon shapefile"""  
  msg = '\nRead the script header... A projected coordinate system required'  
  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.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.