Divide Irregular Polygon into 730 parts

6471
18
Jump to solution
05-06-2014 06:25 AM
TysonBarlow
New Contributor III
Hi guys and gals!

Is there a utility or something that I could use to divide an irregularly shaped polygon into 730 pieces, preferably of equal area? For example, take any given state polygon in the US and divide it into 730 pieces.

Thank you!!
Tags (2)
1 Solution

Accepted Solutions
JoshuaChisholm
Occasional Contributor III
Hello Tyson,

I think this might work (sorry if the code is a little messy):

import arcpy arcpy.env.overwriteOutput = True  #set number of desired cells: x=730  #set feature classes: inFC=r'C:\Path\To\Poly.shp' tempGrid=r'C:\Path\To\tempGrid.shp' #just a temporary file outFC=r'C:\Path\To\grid.shp' #final grid  #set default grid variables: extent = arcpy.Describe(inFC).extent origin_coord=str(extent.XMin)+" "+str(extent.YMin) oppositeCoorner=str(extent.XMax)+" "+str(extent.YMax) y_axis_coord=str(extent.XMin)+" "+str(extent.YMin+10) number_rows=number_columns=0 geometryType = 'POLYGON'  #set starting cellSize and increment: cellSize=max(float(extent.width),float(extent.height)) inc=cellSize  countOutFC=0.0 #needed to start while loop  #Find an appropriate cell size: while countOutFC!=x:     if cellSize==inc:         inc=inc/2.0         continue     cell_width=cell_height=cellSize     arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType)      arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr")     arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC)      countOutFC=float(arcpy.GetCount_management("temp_lyr").getOutput(0))          print "Cell size: "+str(cellSize)+"  Count: "+str(countOutFC)     if countOutFC>x:         print "Overshot..."         cellSize=prevCellSize         inc=inc/2.0     else:         prevCellSize=cellSize         cellSize-=inc  cellSize=prevCellSize print '\nFound an appropriate cell size: '+str(cellSize)  #Find a simplified cell size: print "\nTrying to simplify..." for i in range(len(str(cellSize))):     simpCellSize=round(cellSize,i)     if simpCellSize==0:         continue     cell_width=cell_height=simpCellSize     arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType)     arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr")     arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC)      countOutFC=float(arcpy.GetCount_management("temp_lyr").getOutput(0))     print "Cell size: "+str(simpCellSize)+"  Count: "+str(countOutFC)     if countOutFC==x:         break  if simpCellSize==cellSize:     print "Could not easily find a simplified cell size."     #rerun for last good cellsize (not simplified):     cell_width=cell_height=cellSize     arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType)     arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr")     arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC) else:     print "A simplified cell size: "+str(simpCellSize)  arcpy.Select_analysis("temp_lyr",outFC)  print "Done!!!"


The code is a little slow (it took 35 seconds on my fast machine). Let me know how it goes!

View solution in original post

18 Replies
JoshuaChisholm
Occasional Contributor III
Hello Tyson,

How do you want to split the polygon? By a road network? Do you want the polygon split into little squares, long slivers or something else?

Let me know, and we might be able to figure something out.
0 Kudos
TysonBarlow
New Contributor III
I would like to split it into squares. Being an irregularly shaped polygon, I understand that not all of the 730 squares would be equal area, which is fine. I'd appreciate any help you can give! Something like the screenshot below, but with 730 pieces.
[ATTACH=CONFIG]33624[/ATTACH]
0 Kudos
Zeke
by
Regular Contributor III
Not sure this will work, but look into Create Fishnet. This can create lines that you could use to split the polygon. Part I'm not sure of is how to guarantee 730 cells. You might be able to divide the area of the polygon by 730 to get cell height & width, but if the fishnet is created on the extent of the polygon rather than the feature itself this might not work.
0 Kudos
TysonBarlow
New Contributor III
Yeah, I've tried the fishnet function, and run into the problem you outline here. It creates features on the extent of the polygon and not the feature itself, so I'm having a hard time getting exactly 730 pieces. Thanks for the thought though!
0 Kudos
JoshuaChisholm
Occasional Contributor III
Hello Tyson,

I think this might work (sorry if the code is a little messy):

import arcpy arcpy.env.overwriteOutput = True  #set number of desired cells: x=730  #set feature classes: inFC=r'C:\Path\To\Poly.shp' tempGrid=r'C:\Path\To\tempGrid.shp' #just a temporary file outFC=r'C:\Path\To\grid.shp' #final grid  #set default grid variables: extent = arcpy.Describe(inFC).extent origin_coord=str(extent.XMin)+" "+str(extent.YMin) oppositeCoorner=str(extent.XMax)+" "+str(extent.YMax) y_axis_coord=str(extent.XMin)+" "+str(extent.YMin+10) number_rows=number_columns=0 geometryType = 'POLYGON'  #set starting cellSize and increment: cellSize=max(float(extent.width),float(extent.height)) inc=cellSize  countOutFC=0.0 #needed to start while loop  #Find an appropriate cell size: while countOutFC!=x:     if cellSize==inc:         inc=inc/2.0         continue     cell_width=cell_height=cellSize     arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType)      arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr")     arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC)      countOutFC=float(arcpy.GetCount_management("temp_lyr").getOutput(0))          print "Cell size: "+str(cellSize)+"  Count: "+str(countOutFC)     if countOutFC>x:         print "Overshot..."         cellSize=prevCellSize         inc=inc/2.0     else:         prevCellSize=cellSize         cellSize-=inc  cellSize=prevCellSize print '\nFound an appropriate cell size: '+str(cellSize)  #Find a simplified cell size: print "\nTrying to simplify..." for i in range(len(str(cellSize))):     simpCellSize=round(cellSize,i)     if simpCellSize==0:         continue     cell_width=cell_height=simpCellSize     arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType)     arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr")     arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC)      countOutFC=float(arcpy.GetCount_management("temp_lyr").getOutput(0))     print "Cell size: "+str(simpCellSize)+"  Count: "+str(countOutFC)     if countOutFC==x:         break  if simpCellSize==cellSize:     print "Could not easily find a simplified cell size."     #rerun for last good cellsize (not simplified):     cell_width=cell_height=cellSize     arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType)     arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr")     arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC) else:     print "A simplified cell size: "+str(simpCellSize)  arcpy.Select_analysis("temp_lyr",outFC)  print "Done!!!"


The code is a little slow (it took 35 seconds on my fast machine). Let me know how it goes!
JackDavis
New Contributor

Hi Joshua,

   Do you think it would be possible to modify your script to operate on individual polygons inside of a map instead of the map extents? I'm struggling a bit trying to figure out how to divide a map of many ~ 1 mi sections into 100 ~.1 mi sub-sections with labels in each.

Thanks,

Jack

0 Kudos
JoshuaChisholm
Occasional Contributor III

Hello Jack,

I think this should work (but I didn't have a chance to test).

Let me know how it does. Good luck!

import arcpy 

arcpy.env.overwriteOutput = True 

targetLayers=['farm fields','Municipalities'] #specify all the layers you want (they must have unique names)

mxdpath=r'C:\Path\To\MapDoc.mxd' #point to your mxd

mxd=arcpy.mapping.MapDocument("CURRENT")

lyrs=arcpy.mapping.ListLayers(mxd)

dfName='Layers' #name of the data frame you want to add the new grids to

for lyr in lyrs:

  if not lyr.name in targetLayers:

    continue #Skip all non-target layers

  #set number of desired cells: 

  x=730 

  

  #set feature classes: 

  inFC=lyr

  tempGrid=r'C:\Path\To\tempGrid.shp' #just a temporary file 

  outFC=r'C:\Path\To\grid_'+lyr.name+'.shp' #final grid

 

  #delete tempGrid if exists:

  if arcpy.Exists(tempGrid):

    arcpy.Delete_management(tempGrid)

  

  #set default grid variables: 

  extent = arcpy.Describe(inFC).extent 

  origin_coord=str(extent.XMin)+" "+str(extent.YMin) 

  oppositeCoorner=str(extent.XMax)+" "+str(extent.YMax) 

  y_axis_coord=str(extent.XMin)+" "+str(extent.YMin+10) 

  number_rows=number_columns=0 

  geometryType = 'POLYGON' 

  

  #set starting cellSize and increment: 

  cellSize=max(float(extent.width),float(extent.height)) 

  inc=cellSize 

  

  countOutFC=0.0 #needed to start while loop 

  

  #Find an appropriate cell size: 

  while countOutFC!=x: 

    if cellSize==inc: 

      inc=inc/2.0 

      continue 

    cell_width=cell_height=cellSize 

    arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType) 

  

    arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr") 

    arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC) 

  

    countOutFC=float(arcpy.GetCount_management("temp_lyr").getOutput(0)) 

  

    print "Cell size: "+str(cellSize)+"  Count: "+str(countOutFC) 

    if countOutFC>x: 

      print "Overshot..." 

      cellSize=prevCellSize 

      inc=inc/2.0 

    else: 

      prevCellSize=cellSize 

      cellSize-=inc 

  

  cellSize=prevCellSize 

  print '\nFound an appropriate cell size: '+str(cellSize) 

  

  #Find a simplified cell size: 

  print "\nTrying to simplify..." 

  for i in range(len(str(cellSize))): 

    simpCellSize=round(cellSize,i) 

    if simpCellSize==0: 

      continue 

    cell_width=cell_height=simpCellSize 

    arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType) 

    arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr") 

    arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC) 

  

    countOutFC=float(arcpy.GetCount_management("temp_lyr").getOutput(0)) 

    print "Cell size: "+str(simpCellSize)+"  Count: "+str(countOutFC) 

    if countOutFC==x: 

      break 

  

  if simpCellSize==cellSize: 

    print "Could not easily find a simplified cell size." 

    #rerun for last good cellsize (not simplified): 

    cell_width=cell_height=cellSize 

    arcpy.CreateFishnet_management(tempGrid,origin_coord,y_axis_coord,cell_width,cell_height,number_rows,number_columns,oppositeCoorner,"LABELS","#",geometryType) 

    arcpy.MakeFeatureLayer_management(tempGrid, "temp_lyr") 

    arcpy.SelectLayerByLocation_management("temp_lyr", "INTERSECT", inFC) 

    else: 

      print "A simplified cell size: "+str(simpCellSize) 

  

    arcpy.Select_analysis("temp_lyr",outFC)

    #add new data to mxd

    df = arcpy.mapping.ListDataFrames(mxd, dfName)[0]

    newlayer = arcpy.mapping.Layer(outFC

    arcpy.mapping.AddLayer(df, newlayer, "AUTO_ARRANGE")

 

print "Done!!!" 

Edit (July 14, 2014): Added lines 25-27

JackDavis
New Contributor

Thanks for the reply Joshua! I didn't have much time to try this out Friday before I left the office, but on initial run, I got an error on the AddLayer call on line 94. The second argument 'addLayer' was undefined. I'll read through it again Monday morning and see if I can figure it out...I may have left something out.

Thanks again,

Jack

0 Kudos
JoshuaChisholm
Occasional Contributor III

Hello Jack,

I had made an edit to my code shortly after posting. You must have grabbed the wrong code. Line 94 should read: arcpy.mapping.AddLayer(df, newlayer, "AUTO_ARRANGE"

Let me know if that works for you!

0 Kudos