Friday, August 29, 2014

Speed Up Polygon Grid Generation

There are many tools that generate regular uniform grids with geoprocessing, but I found a novel was of generating grids pretty quickly for large areas.  Rasters are fast and can be created on the fly easily using numpy.  To make the raster a feature class, ArcPy has a Raster to Polygon tool that should be available at all levels of desktop.

import arcpy
import numpy
import math
fc = r"c:\temp\somedata.shp" # assume projection is projected coordinate system not Lat/Long where measurements are in decimal degrees
width = 10000 # meters
height = 10000 # meters
#  Gets the feature classes' extent
#
desc = arcpy.Describe(fc)
extent = desc.extent
#   Set the Mask
#
arcpy.env.mask = fc
#   Calculate the number of columns
#

dist_x = math.sqrt(math.pow((extent.XMax - extent.XMin), 2))
dist_y = math.sqrt(math.pow((extent.YMax - extent.YMin),2))
col_x = int(dist_x / width)
col_y = int(dist_y / height)
if col_x % distanceMeters > 0:
   col_x += 1
if col_y % distanceMeters > 0:
  col_y += 1
#   Construct the structure numpy array with random values
#
numpyArray = (numpy.random.permutation(col_x * col_y) + 1).reshape(col_y, col_x)
#   Convert array to raster 
#
myData = arcpy.NumPyArrayToRaster(numpyArray,
                                  arcpy.Point(extent.XMin, extent.YMin),
                                  width, height)
#   Convert the raster to a polygon feature class
#
grid = r"%s\gridData" % arcpy.env.scratchGDB
grid = arcpy.RasterToPolygon_conversion(in_raster=myData,
                                            out_polygon_features=grid,
                                            simplify="NO_SIMPLIFY",
                                            raster_field="VALUE")[0]
The code takes the extent of the data and creates a grid. The mask environmental parameter essentially clips the raster to the shape of the feature class. This way you only gets cells within the given polygon area. If you wanted cells to fall outside the feature classes' shape, then just comment out that line.

So that is that.


If you haven't done so, check out ArcREST.
Also vote for my idea for a GUI builder for Python Add-Ins.

1 comment:

Rob Choucroun said...

FYI, you haven't defined distanceMeters anywhere. Did you mean to use dist_x and dist_y?