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) """ 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 = 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""" 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.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.

Jeff nice implementation...now for other shapes... see my blog post

Hexagons, Rectangles... Sampling Frameworks