A new tool for creating sampling hexagons

by Tim Whiteaker, Research Associate, Center for Research in Water Resources, The University of Texas at Austin, twhit@mail.utexas.edu.

I recently needed to create a mesh of hexagons covering a study area in the Arctic for a research proposal.  These hexagons define sampling locations, helping to ensure that all regions within the study area are represented by the sampling results.

Such sampling schemes are often used by wildlife biologists, aquatic chemists, and other researchers.  However, a hexagon creation tool is not included with ArcGIS, so I decided to create my own using ModelBuilder.

Sampling locations defined by hexagons

This tool creates a mesh of hexagons overlapping a study area. I only had to create one script tool for the model, relying heavily on standard ArcGIS tools in the model to keep code (and future maintenance!) to a minimum.

Create hexagons model

The model tool is fairly simple.  The model uses a script tool to create a mesh of points.  The points are spaced such that if you run the Create Thiessen Polygons tool included with ArcGIS, regular hexagons (those of equal side length with internal angles of 120 degrees) are generated.  The model runs these tools and then intersects the result with the study area to produce the final hexagon mesh.

Study area split up into hexagons

If you don’t like the result, you can choose a different hexagon size and run the tool again, or you can use ArcMap’s editor to fine tune the hexagon locations.  You can also dig into the relatively small bit of script code or change tools within the model to get the result you’re looking for.

I’ve shared this tool as a geoprocessing package in ArcGIS Online.  I welcome feedback and in fact, I would love to hear ideas on how to tackle some challenges I encountered while developing the tool.  You can read more about those challenges on the tool’s ArcGIS Online page.  Happy hexagon-ing!

This entry was posted in Analysis & Geoprocessing, Hydro, Local Government, Location Analytics, National Government, Petroleum, Public Safety, Spatial Statistics, State Government, Transportation, Uncategorized, Water Utilities and tagged , , , , . Bookmark the permalink.

Leave a Reply

10 Comments

  1. graber says:

    That is a good idea using Thiessen Polygons. I’ll have to try it and see if it is any faster. I attached a script below that I wrote about a year ago that generates the hexs using math, which I found was not that much after the center points had already been created.

    # Import system modules
    import sys, string, os, time, string, math, arcpy
    from arcpy import env
    import datetime

    starttime = datetime.datetime.now()

    env.overwriteOutput = True

    # Change parameters below to suit your needs
    extentlayer = r’C:\TestData\tydixton_watershed.shp’ #Input to get Extent from
    output_tot = env.scratchGDB + r”/outputhex” #Ouput Location
    hexareastr = 20 #area of hexegons (in hectare if projection is in meters)
    hexarea = float(hexareastr) * float(10000) #area conversion betoween ha and sq meteres.

    triarea = hexarea / 6
    oneside = math.sqrt((triarea * 4) / math.sqrt(3))
    halfside = oneside / 2
    onehi = math.sqrt((oneside * oneside) – ((halfside) * (halfside)))
    longway = oneside * 2
    shortway = onehi * 2

    desc = arcpy.Describe(extentlayer)
    sr = desc.spatialreference

    minx = desc.extent.XMin
    miny = desc.extent.YMin
    maxx = desc.extent.XMax
    maxy = desc.extent.YMax
    xrange = maxx – minx
    yrange = maxy – miny

    distancex = oneside + halfside
    distancey = shortway

    th = int((xrange + longway + oneside + 0.02) / distancex)* int((yrange + shortway + shortway) / distancey)
    oum = str(th) + ” Will Cover the Extent Area”
    print oum
    arcpy.AddMessage(oum)
    arcpy.AddMessage(“Generating Hexagons…”)

    featureList = []

    for xc in range(int((xrange + longway + oneside + 0.02) / distancex)):
    centerx = (minx + (xc * distancex)) – 0.01
    if (xc%10) == 0:
    arcpy.AddMessage(” On row ” + str(xc) + ” of ” + str(int((xrange + shortway) / distancex)))
    print ” On row ” + str(xc) + ” of ” + str(int((xrange + shortway) / distancex))
    if (xc%2) == 0:
    ystart = miny
    else:
    ystart = miny – onehi
    for yc in range(int((yrange + shortway + shortway) / distancey)):
    centery = ystart + (yc * distancey)
    #Hex Generation
    gonarray = arcpy.Array()

    pt1 = arcpy.Point()
    pt1.X = centerx + halfside
    pt1.Y = centery + onehi
    gonarray.add(pt1)

    pt2 = arcpy.Point()
    pt2.X = centerx + oneside
    pt2.Y = centery
    gonarray.add(pt2)

    pt3 = arcpy.Point()
    pt3.X = centerx + halfside
    pt3.Y = centery – onehi
    gonarray.add(pt3)

    pt4 = arcpy.Point()
    pt4.X = centerx – halfside
    pt4.Y = centery – onehi
    gonarray.add(pt4)

    pt5 = arcpy.Point()
    pt5.X = centerx – oneside
    pt5.Y = centery
    gonarray.add(pt5)

    pt6 = arcpy.Point()
    pt6.X = centerx – halfside
    pt6.Y = centery + onehi
    gonarray.add(pt6)

    pt7 = arcpy.Point()
    pt7.X = centerx + halfside
    pt7.Y = centery + onehi
    gonarray.add(pt7)

    op = arcpy.Polygon (gonarray, sr)
    gonarray.removeAll()

    featureList.append(op)

    tmplay = “in_memory/all_hex” #env.scratchGDB + r”/all_hex”
    arcpy.CopyFeatures_management(featureList, tmplay)

    arcpy.MakeFeatureLayer_management(tmplay, ‘hexes’)
    arcpy.MakeFeatureLayer_management(extentlayer, ‘extent’)
    arcpy.SelectLayerByLocation_management(‘hexes’, ‘intersect’, ‘extent’)

    arcpy.CopyFeatures_management(‘hexes’, output_tot)

    #gp.calculatefield_management(output_tot, “UNIT_ID”, “[" + oidnam + "] + 1″)
    arcpy.AddMessage(“HEXGEN Complete”)
    print “Complete”

    del desc
    del arcpy

    endtime = datetime.datetime.now()

    td = endtime – starttime

    print “Total Processing Time: “, td

    • twhiteaker says:

      @graber
      Thanks for sharing your hexagon code. I was curious about which would be faster so I turned your code into a script tool and removed some of the AddMessages and deleted some other lines (like import statements that weren’t critical such as datetime) to make it run as fast as possible, and I also let it grab parameters from a script tool instead of hard-coding them, so that your tool would run in a manner similar to mine. I then created a model and added your tool to it so that both examples would be encapsulated by their own model. I ran each model three times to produce hexagons with side length 50 meters, area 0.64951905283 hectares, in the rough outline of Torrey Pines that I included in the package. This resulted in 1089 hexagons for my tool and 1090 for yours, reflecting that our “hexagon origin” is different but that results are comparable.

      In all three trials, my tool ran for 9 seconds whereas yours ran for 15 seconds as reported by the geoprocessing window. The difference may be due to some of your hexagons being written twice with the two CopyFeatures calls, or maybe the Esri Thiessen algorithm is really fast. The Esri Thiessen tool might also be getting a slight boost if the code it’s written in runs faster than Python. The difference may also be due to how I set your code up to run as a script tool, so I’ve included your modified code below for reference.

      On another note, try setting the map projection to WGS84 (or any geographic system) while your data are in a projected coordinate system. When you run the tool, you may get unexpected results since arcpy.Describe returns results in map units, not feature units. Any ideas to work around this?

      I’ve included your modified code below. If anyone tries to execute this, you may need to replace single quotes, double quotes, and minus signs, as these might be converted to other characters when posted in the blog. I had to replace those characters from graber’s post to make it work for me. I’ll see if using the HTML code tag helps.


      # Import system modules
      import sys, string, os, string, math, arcpy

      # Change parameters below to suit your needs
      #extentlayer = r'C:\TestData\tydixton_watershed.shp' #Input to get Extent from
      extentlayer = arcpy.GetParameterAsText(0)
      #output_tot = env.scratchGDB + r"/outputhex" #Ouput Location
      output_tot = arcpy.GetParameterAsText(2)
      #hexareastr = 20 #area of hexegons (in hectare if projection is in meters)
      hexareastr = arcpy.GetParameterAsText(1)
      hexarea = float(hexareastr) * float(10000) #area conversion between ha and sq meteres.

      triarea = hexarea / 6
      oneside = math.sqrt((triarea * 4) / math.sqrt(3))
      halfside = oneside / 2
      onehi = math.sqrt((oneside * oneside) - ((halfside) * (halfside)))
      longway = oneside * 2
      shortway = onehi * 2

      desc = arcpy.Describe(extentlayer)
      sr = desc.spatialreference

      minx = desc.extent.XMin
      miny = desc.extent.YMin
      maxx = desc.extent.XMax
      maxy = desc.extent.YMax
      xrange1 = maxx - minx
      yrange = maxy - miny

      distancex = oneside + halfside
      distancey = shortway

      th = int((xrange1 + longway + oneside + 0.02) / distancex)* int((yrange + shortway + shortway) / distancey)
      oum = str(th) + " Will Cover the Extent Area"
      arcpy.AddMessage(oum)
      arcpy.AddMessage("Generating Hexagons...")

      featureList = []

      for xc in range(int((xrange1 + longway + oneside + 0.02) / distancex)):
      centerx = (minx + (xc * distancex)) - 0.01
      if (xc%2) == 0:
      ystart = miny
      else:
      ystart = miny - onehi

      for yc in range(int((yrange + shortway + shortway) / distancey)):
      centery = ystart + (yc * distancey)
      #Hex Generation
      gonarray = arcpy.Array()

      pt1 = arcpy.Point()
      pt1.X = centerx + halfside
      pt1.Y = centery + onehi
      gonarray.add(pt1)

      pt2 = arcpy.Point()
      pt2.X = centerx + oneside
      pt2.Y = centery
      gonarray.add(pt2)

      pt3 = arcpy.Point()
      pt3.X = centerx + halfside
      pt3.Y = centery - onehi
      gonarray.add(pt3)

      pt4 = arcpy.Point()
      pt4.X = centerx - halfside
      pt4.Y = centery - onehi
      gonarray.add(pt4)

      pt5 = arcpy.Point()
      pt5.X = centerx - oneside
      pt5.Y = centery
      gonarray.add(pt5)

      pt6 = arcpy.Point()
      pt6.X = centerx - halfside
      pt6.Y = centery + onehi
      gonarray.add(pt6)

      pt7 = arcpy.Point()
      pt7.X = centerx + halfside
      pt7.Y = centery + onehi
      gonarray.add(pt7)

      op = arcpy.Polygon (gonarray, sr)
      gonarray.removeAll()

      featureList.append(op)

      tmplay = "in_memory/all_hex" #env.scratchGDB + r"/all_hex"
      arcpy.CopyFeatures_management(featureList, tmplay)

      arcpy.MakeFeatureLayer_management(tmplay, 'hexes')
      arcpy.MakeFeatureLayer_management(extentlayer, 'extent')
      arcpy.SelectLayerByLocation_management('hexes', 'intersect', 'extent')

      arcpy.CopyFeatures_management('hexes', output_tot)

      #gp.calculatefield_management(output_tot, "UNIT_ID", "[" + oidnam + "] + 1")
      arcpy.AddMessage("HEXGEN Complete")

      • Kenneth Field says:

        Hex binning tools do exist for ArcGIS. Click the following links to blogs that describe and use the tool:

        Create hexagons tool
        Using a binning technique for point based multiscale web maps

        • twhiteaker says:

          Thanks for the reply, Kenneth. Actually, for a while I had been using the very same beehive tool that you linked to, except I got it from here. I downloaded it many months ago and found that I had to tweak the code to get it to run without errors, and even then it ran slowly for me which was part of the motivation for writing my own tool. I did a fresh download and am glad to see that it runs without errors now. Using their sample shapefile with their default hexagon size of 500, it took 1 hour 49 minutes 52 seconds to construct 38110 hexagons. I ran my tool with the same parameters and it took 45 seconds. I also discovered a couple of bugs with my tool (model didn’t work with shapefiles; right side doesn’t extend far enough) and have updated the package with fixes for those bugs.

          There are some other hexagon tools out there. Today I found Hexagon Maker, an Add-In that runs much faster than mine, although I found it a bit quirky in how you set up the tool to run, e.g., it adds hexagons to your original study area feature class; it only allows one study area polygon to be selected. But it sure is fast! There’s also one from Jenness Enterprises although I couldn’t install it (.exe, .bat) due to security restrictions.

          Thanks for sharing those links. I found your binning post quite interesting.

  2. tjksmit says:

    I’m interested in the way you created the sample points for the hexagons. It is quite a issue to create a proper grid of points in ArcGIS, and the script you are using looks very promising and is quick. Is there possibility to look into it?

    • Tim Whiteaker says:

      The script is available within the geoprocessing package. When you download and double-click the package, ArcMap should open, and a “Shared” category should appear in the Geoprocessing Results window. Expand Shared and Create_Hexagon_Tessellation to see the Create Hexagon Tessellation tool. Right-click the tool and click Locate. The Catalog window will navigate to the tool location, and you can expand the create_hexagon_tessellations subfolder to find the original model and script.

      Or, you could just navigate directly to where packages are installed. For me, the path to the script looks like:
      [user]\Documents\ArcGIS\Packages\Create_Hexagon_Tessellation\v101\create_hexagon_tessellations

  3. slevental says:

    Thanks for posting!

  4. bsneish says:

    Using 10.2 and a feature class from a file geodatabase,
    the line
    if desc.dataType == “FeatureLayer”:
    desc = arcpy.Describe(desc.featureClass.catalogPath)
    breaks for me, resulting in an error message. If i remove that check it works. (i made sure my data frame was the same coordinate as the layer)

    • ruaricolquhoun says:

      Hi Tim,
      How would you go about selecting an area say in Australia for making a grid of hexagons? I have tried navigating to a polygon outline using your tool an ArGIS Desktop 10.1, however I keep getting error messages (however it works fine for your example of Torrey Pines State Reserve). Would I have to edit the script to navigate to my polygon shape file?
      Thanks,

      Ruari

      • Tim Whiteaker says:

        In the off chance that you’re still working on this, I’m curious as to why the tool fails and would like to help. I just created hexagons for Australia using a shapefile from Natural Earth. What data are you using for your Australia boundary? Feel free to email me directly.