Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Given: A single part polygon, total number of points to place inside that polygon, and table specifying the size of each subset (number of subsets n >= 1) of that total number of points and each subset's minimum spacing.

Find: A (semi-)random distribution of points inside the polygon such that no point is closer than its minimum spacing to any other point in the polygon.

For example, you are given a polygon and 100 points. Of the set of 100 points there are three disjoint subsets--A, B and C. Subset A contains 50 points having a minimum spacing of 7 m, subset B contains 25 points having a minimum spacing of 4 m, and subset C contains 25 points having a minimum spacing of 2 m. Distribute the points such that no point in subset A is within 7 m of any of the remaining 99 points, no point in subset B is within 4 m, no point in subset C is within 2 m.

One possible method to solve this would be to start with subset A, distribute the points randomly using Create Random Points, create a 7 m buffer around the points, erase those buffers from the original polygon and use the erased polygon to randomly distribute the points in subset B. . .repeat through all subsets.

I'm looking for a way to implement this (or a superior) solution using scripting or model builder. Ideally, the user would input the number of points in each subset, the minimum spacing for each subset, and the polygon feature class into which the points would be placed.

share|improve this question
1  
I do not understand the objection to your proposed method. When you exclude a 7 m buffer around the points in A, the points in B and C can both come arbitrarily close to 7 m to those in A--but no closer. The sum of buffer radii does not appear relevant. –  whuber Aug 16 '13 at 20:16
    
As long as the buffer sizes in the subsets are decreasing (or remaining constant), it should work. –  Paul Aug 16 '13 at 20:25
    
@whuber I see your point and you're correct. I suppose I was incorrectly visualizing the buffers of subsets B, C, etc being placed inside the ever-shrinking polygon as opposed to points. Should I take this down? –  Ben Aug 16 '13 at 20:26
    
Although I think it's an interesting problem, unless you broaden the question to probe for alternative implementations you probably ought to delete it because at present there may be little more one can say. –  whuber Aug 16 '13 at 20:28
    
If you are looking for a script implementation, you could edit your question accordingly. –  Paul Aug 16 '13 at 20:30

2 Answers 2

up vote 4 down vote accepted

A script implementation of this would be relatively straightforward:

  1. Create random points for subset A, creating points α.
  2. Buffer α by distance A, creating polygon buff1.
  3. Use Erase to remove buff1 from erase0 (input polygon), creating polygon erase1.
  4. Create random points for subset B, based on erase1, creating points β.
  5. Buffer β by distance B, creating polygon buff2.
  6. Use Erase to remove buff2 from erase1, creating erase2.
  7. Create random points for subset C, based on erase2, creating points γ.
  8. Use Append to bring α, β, and γ together.

Repeat steps 2-4 as needed..

The following script should get you started:

from os.path import join
import arcpy

#name input polygon erase0
inputpoly = "path to polygon FC"
#final output will be named point0
outpath = "path to output directory"  
arcpy.env.workspace = outpath

subsets = [[50, "7 meters"], [25, "4 meters"], [25, "2 meters"]]

arcpy.CreateRandomPoints_management(outpath, "point0", "erase0", "",
                                    subsets[0][0], subsets[0][1])

pointnames = []
for i in xrange(len(subsets)-1):        
    arcpy.Buffer_analysis("point{}".format(i), "buff{}".format(i),
                          subsets[i][1], dissolve_option="ALL")
    arcpy.Erase_analysis("erase{}".format(i), "buff{}".format(i), "erase{}".format(i+1))
    arcpy.CreateRandomPoints_management(outpath, "point{}".format(i+1), "erase{}".format(i+1),
                                        "", subsets[i+1][0], subsets[i+1][1])

    pointnames.append("point{}".format(i+1))


points = ";".join(pointnames)
arcpy.Append_management(points, "point0")
share|improve this answer
    
This is helpful. Is there a way to may it more general? Specifically, I'd like to incorporate it into a tool where the user can specify how many subsets there are, how many points in each subset, and the spacing in each subset. I'm new to python so I'm not sure how to make subsets sized dynamically. I don't expect the finished code, but any specific resources would be helpful. –  Ben Aug 20 '13 at 16:10
    
@Ben, it would be very simple to do it like that. One easy way would be to input a comma delimited text file for your three parameters. Since I'm looping over subsets it can be as small or as big as needed. Look into reading text files with open(). –  Paul Aug 20 '13 at 18:13

Here's the code for a python toolbox I came up with to randomly place a finite number of sets of circles (whose centers are, naturally, points) with different radii into a finite number of arbitrarily shaped polygons. The code is based on the seed Paul provided in his answer above. The parameter param1 is a GPValueTable data type, which allows you to specify the number of columns and each column's own data type. It's actually just a list of lists, but each list can have multiple data types in it.

The idea for the problem was based on figuring out spacing for placing plants which require a minimum spacing from other plants, but each species could have a different spacing. . .thus the language centered around "Species".

The way the tool is currently set up produces a feature class of only circles, but point feature classes are created, and it is easy enough to remark out the lines of code that delete the point feature classes.

There is some validation included, which keeps you from executing the tool unless the radii/spacing is in descending order in the value table.

It works fine as a brute force tool, but only if the workspace is a geodatabase (not just a folder), and I can't figure out how to add the resulting feature class to a dataframe in the current .mxd using this python toolbox. It's probably a bit clunky, but it's my first try.

import arcpy
import arcpy.mapping
import os.path


class Toolbox(object):
def __init__(self):
    """Define the toolbox (the name of the toolbox is the name of the
    .pyt file)."""
    self.label = "Restoration Planting Toolbox"
    self.alias = ""

    # List of tool classes associated with this toolbox
    self.tools = [SpeciesSpacing]


class SpeciesSpacing(object):
def __init__(self):
    """Define the tool (tool name is the name of the class)."""
    self.label = "Species Spacing Tool"
    self.description = "This tool takes a list of species, number of individuals in each species, their minimum spacing, and a bounding polygon and distributes circles within
    the polygon such that no species' minimum spacing is violated."
    self.canRunInBackground = False

def getParameterInfo(self):
    """Define parameter definitions"""

    # First parameter
    param0 = arcpy.Parameter(
       displayName="Bounding polygon feature class:",
            name="in_features",
            datatype="GPLayer",
            parameterType="Required",
            direction="Input")

    # Second parameter
    param1 = arcpy.Parameter(
            displayName="Add each species by name ***LIST SPECIES IN ORDER OF DESCENDING MINIMUM SPACING***:",
            name="species_table",
            datatype="GPValueTable",
            parameterType="Required",
            direction="Input")

    param1.columns = [['String','Species Name'],['Long','Number of Individuals'],['Long','Minimum Spacing (feet)']]

    param2 = arcpy.Parameter(
        displayName="Output workspace:",
        name="output_ws",
        datatype="DEWorkspace",
        parameterType="Required",
        direction="Input")

    param3 = arcpy.Parameter(
        displayName="Name of output feature class:",
        name="output_fc",
        datatype="GPString",
        parameterType="Required",
        direction="Input")

    parameters = [param0,param1,param2,param3]
    return parameters

def isLicensed(self):
    """Set whether tool is licensed to execute."""
    return True

def updateParameters(self, parameters):
    """Modify the values and properties of parameters before internal
    validation is performed.  This method is called whenever a parameter
    has been changed."""

    if parameters[1].values is not None:

        vtCheck = parameters[1].values

        for i in range(len(vtCheck)-1):
            if vtCheck[i][2] < vtCheck[i+1][2]:
                    parameters[2].enabled = False
                    parameters[3].enabled = False


            elif vtCheck[i][2] >= vtCheck[i+1][2]:
                    parameters[2].enabled = True
                    parameters[3].enabled = True

    return


def updateMessages(self, parameters):
    """Modify the messages created by internal validation for each tool
    parameter.  This method is called after internal validation."""
    return

def execute(self, parameters, messages):
    """The source code of the tool."""
    # Specify the current map
    mxd = arcpy.mapping.MapDocument("CURRENT")

    # get the data frame 
    df = arcpy.mapping.ListDataFrames(mxd,"*")[0]

    # Specify input polygon
    inputPoly = parameters[0].valueAsText

    # Set output path
    outpath = parameters[2].valueAsText
    arcpy.env.workspace = outpath

    # Specify final feature class name
    outname = parameters[3].valueAsText

    fullpath = os.path.join(outpath,outname)

    messages.addMessage("\nOutput is located at {0}".format(fullpath))


    # Get values from value table parameter
    vt=parameters[1].values

    # Make a copy of the area
    arcpy.CopyFeatures_management(inputPoly,"area0")

    # Get the spatial reference of the input polygon
    sr = arcpy.Describe(inputPoly).spatialReference

    # Make polygon feature class to place circles in
    arcpy.CreateFeatureclass_management(outpath,outname,"POLYGON",spatial_reference=sr)
    # Add field for species name
    arcpy.AddField_management(outname,"Species","text")

    # Create the initial set of random points
    arcpy.CreateRandomPoints_management(outpath, "{0}_points".format(vt[0][0]), "area0", "",
                                vt[0][1], "{0} feet".format(vt[0][2]))

    for i in range(len(vt)):
        # Create buffers
        arcpy.Buffer_analysis(outpath+"/{0}_points".format(vt[i][0]), "{0}_buffers".format(vt[i][0]),
                      "{0} feet".format(vt[i][2]), dissolve_option="NONE")
        messages.addMessage("\nPlacing {0}. . .".format(vt[i][0]))
        # Add 'Species' field
        arcpy.AddField_management("{0}_buffers".format(vt[i][0]),"Species","text")
        # Calculate 'Species' field based on value table input
        arcpy.CalculateField_management("{0}_buffers".format(vt[i][0]),"Species",'"{0}"'.format(vt[i][0]))
        # Erase buffers from 'area{i}' to make 'area{i+1}'
        arcpy.Erase_analysis("area{0}".format(i), "{0}_buffers".format(vt[i][0]), "area{0}".format(i+1))
        # Append '{species}_buffers' to 'randomized_planting' feature class
        arcpy.Append_management("{0}_buffers".format(vt[i][0]),outname,"NO_TEST")
        # Delete 'area{i}' feature class
        arcpy.Delete_management(outpath+"/area{0}".format(i))
        messages.addMessage("\nCleaning up. . .")
        # Delete '{species}_buffers' feature class
        arcpy.Delete_management(outpath+"/{0}_buffers".format(vt[i][0]))
        # Delete '{species}_points' feature class
        arcpy.Delete_management(outpath+"/{0}_points".format(vt[i][0]))

        # Generate random points in the next iteration of the area (except for the final iteration)
        if i < len(vt)-1:
            arcpy.CreateRandomPoints_management(outpath, "{0}_points".format(vt[i+1][0]), "area{0}".format(i+1), "",
                                vt[i+1][1], "{0} feet".format(vt[i+1][2]))

    # Delete final iteration of the area
    arcpy.Delete_management(outpath+"/area{0}".format(len(vt)))
    messages.addMessage("\nWait for it. . .")

    newLayer = arcpy.MakeFeatureLayer_management(fullpath, "Random")
    add_layer = arcpy.mapping.Layer("Random")
    arcpy.mapping.AddLayer(df,add_layer)

    # Refresh things
    arcpy.RefreshActiveView()
    arcpy.RefreshTOC()
    #del mxd, df

    return
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.