**Numpy Snippets**

**Updated: **2016-09-09

The purpose of this post is to show how numpy can play nicely with arcpy to produce geometry and perform tasks like shape translation, rotation and scaling which are not readily available in arcpy's available functions.

To pay homage to ArcMap's ...** Fish Net** ... I propose a very scaled down version aptly named **Phish_Nyet**. A fuller version will appear when ArcScript 2.0 opens. This is a demo script. All you do is vary the parameters in the **demo()** function within the script.

Only the rectangular option will be presented here, however, it is possible to make all kinds of sampling grids, including hexagonal ones as shown above.

The following points apply:

- Output is to a shapefile.
- The output file must have a projected coordinate system. To create grids in Geographic coordinates, use Fishnet in ArcMap. Why? because invariably people create grids in Geographic coordinates, then project them without densifying the grid blocks. This results in shapes which are in error because the curvature associated with lines of latitude are not accounted for.
- A corner point is specified as the origin of the net. One can determine which corner to use by specifying a dX and dY with positive and/or negative values. In this fashion, you can get your corner to be the top or bottom, left or right of the output. If you want the middle to be the origin...do the math and get a corner.
- The output is controlled by the cell widths (dX and dY), the number of columns and rows an (i.e. the X and Y directions) and a rotation angle, which is positive for clockwise rotation.

Notes:

- **grid_array** - does the numpy magic of generating the grid shapes. I have tried to be as verbose as possible. In short, I generate a

seed shape and propagate it. I have intentionally kept **rotate** and **output_polygons** as visible function so you can see how they work.

A fuller version will surface as I stated when ArcScripts 2.0 appears. Change the parameters in the demo() function and run it.

`"""`

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

Of course other regular geometric shapes can be generated in a similar fashion, but not all may pack like rectangles and others do.

I've had fun messing around with this. I have successfully created a grid of hexagons. Thanks for indirectly teaching me about list comprehensions.