Divide Irregular Polygon into 730 parts

6545
18
Jump to solution
05-06-2014 06:25 AM
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
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!
18 Replies
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.
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]
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.
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!
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!
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

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

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

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!