Creating wedge-shaped territories

[JuanetteWillis]

Good morning, I am trying to create new territories for our inspectors based on the number of locations to be inspected. Our office is centrally located in the county. We are wanting to have wedge-shaped territories radiating from our office, like pieces of a pie. Any suggestions on how to have the GIS assist with the division?

Thanks,

Juanette

[XanderBakker]

Cool! I should have used numpy, but didn't, since I'm not that experienced in it as you are ...

So what I did still has a lot of room for enhancement and fails to create the last wedge. Not sure why, and I don't have the time to look into it. I will post what I have and if anyone wants to play around with it, please do and post back if you get a decent result.

I drew a polygon with the area and created some random points (4400). Additionally I create another point featureclass with the point representing the office.

The code works as follows:

• pnt_off is the point geometry of the office
• it creates a list of slices (19*2 and 6*1) representing the 19 large and 6 smaller areas to be created
• It loops through the slices
  • for each slice it starts to create the wedge using small angle steps
  • it creates a polygon of the wedge and counts the points inside
  • it is is close enough or higher than the number of points searched for it breaks and add the polygon to a list
  • at the end it writes the list of wedges (polygons) to a featureclass
• I manually did a union with the area, selected those records where the FID_Area <> -1 and joined the points to the wedges to count the results.

I´m getting an error when I try to activate the advanced editor, so I will post the code later...

[XanderBakker]

Code used:

import arcpy
import math
import os

def main():
    arcpy.env.overwriteOutput = True

    fc_off = r"D:\Xander\GeoNet\Wedges2\data.gdb\Office"
    fc_pnt = r"D:\Xander\GeoNet\Wedges2\data.gdb\RandomPoints"
    fc_tst = r"D:\Xander\GeoNet\Wedges2\data.gdb\test07"

    sr = arcpy.Describe(fc_pnt).spatialReference
    ws, name = os.path.split(fc_pnt)

    pnt_off = arcpy.da.SearchCursor(fc_off, ("SHAPE@")).next()[0]
    cnt_pts = Count(fc_pnt)

    radius = 6000 # value higher than the max distance of point to office
    slices = [2] * 19 + [1] * 6 # create the slices and corresponding sizes
    tol = 3 # a tolerance value (will stop if current value >= serach value - tolerance)
    steps = 1000 # number of steps (each bearing is 360 / steps)

    lst, lst_pols  = [], []
    i = 0
    last_step = 0
    for s in slices:
        i += 1
        print "Slice: {0},   Size: {1}".format(i, s)
        print "Percentage: {0}%,  Points: {1}".format(round(float(s)/sum(slices)*100, 2), round(float(s)/sum(slices) * cnt_pts,0))
        max_pnts = round(float(s)/sum(slices) * cnt_pts,0)

        for step in range(last_step, steps):
            bearing = float(step) / float(steps) * 360
            print " - bearing: {0}".format(round(bearing,2))
            pnt = getPoint(pnt_off.firstPoint, bearing, radius)

            if step == last_step:
                # create list of first points
                pnt0 = arcpy.Point(pnt.X, pnt.Y)
                lst = [pnt_off.firstPoint, arcpy.Point(pnt.X, pnt.Y), pnt_off.firstPoint]
            else:
                # update list, create polygon and count points in polygon
                lst.insert(len(lst)-1, arcpy.Point(pnt.X, pnt.Y))
                polygon = arcpy.Polygon(arcpy.Array(lst), sr)
                cnt = CountInPolygon(fc_pnt, polygon, ws)
                print "   - point count: {0}".format(cnt)
                last_step = step
                if cnt >= max_pnts - tol:
                    print "     - add wedge"
                    lst_pols.append(polygon)
                    break

    # add last polygon
    lst.insert(len(lst)-1, arcpy.Point(pnt0.X, pnt0.Y))
    polygon = arcpy.Polygon(arcpy.Array(lst), sr)
    print "     - add last wedge"
    lst_pols.append(polygon)

    # create polygon and count
    arcpy.CopyFeatures_management(lst_pols, fc_tst)

def getPoint(pnt, bearing, radius):
    angle = math.radians(90 - bearing)
    dist_x, dist_y = (radius * math.cos(angle), radius * math.sin(angle))
    return arcpy.Point(pnt.X + dist_x, pnt.Y + dist_y)

def Count(fc):
    return int(arcpy.GetCount_management(fc).getOutput(0))

def CountInPolygon(fc_pnt, polygon, ws):
    arcpy.env.overwriteOutput = True
    fc_tmp = "IN_MEMORY\polygon"
    arcpy.CopyFeatures_management([polygon], fc_tmp)
    arcpy.MakeFeatureLayer_management(fc_tmp, "pol", None, ws)
    arcpy.MakeFeatureLayer_management(fc_pnt, "pnts", None, ws)
    arcpy.SelectLayerByLocation_management("pnts", "INTERSECT", "pol")
    return Count("pnts")

if __name__ == '__main__':
    main()

[JuanetteWillis]

Oh my...this is all great and has my mind spinning with the various possibilities.  I am still deep in the depths of hand placing the last couple of hundred locations that didn't geocode. But I am really looking forward to trying at least a couple of these ideas. I have not used code before in ArcMap but am excited about the possibilities. Unfortunately GIS is only a tiny portion of my job duties (and one that I specifically sought out myself) so I don't get near enough time to work on much beyond the basics. I do really appreciate all this input and frankly am a bit in awe of the great responses. Thank you so much.  (And feel free to keep additional ideas coming if you are so inspired.) 

[JustinHollenbach1]

You could just make a shapefile with equal area wedges and then spatially join or intersect the point file to the newly created shapefile to get points that are within a certain area.

[JonathanGillespie]

This may or may not be what you're looking for, but I had a similar problem a couple of years ago.  The solution was a custom script written by someone else.  Basically it creates a circular buffer around a central point (or points) split into wedges using user defined angles.

Creating Wedge/Pie Shaped &amp;amp;quot;Buffers&amp;amp;quot; Around Points

As far as I remember it works with the basic license, though I don't recall which extensions are required.

[XanderBakker]

It is an interesting script and you gotta love the message at the end: "** SUSAN IS VERY COOL **". However, it is based on the user defining the angles. This can be complex (and time consuming for the trial and error) if you have unevenly distributed points and a central point which is not in the center. The script I posted (although it has its flaws) will determine the angle while evaluating a growth of the polygon by performing a count on the points inside.

[DavidBollinger]

equal-angle and equal-area approaches are useful first approximations, but there's a "relatively easy" way to approximate what is probably truly desired (equal-distribution) via iterative methods, requiring only some small modifications if given a script like the one you posted earlier.

finding a set of N radial dividing lines that would exactly equally divide a population P of sites would be a non-trivial optimization problem.  however, if we're already admitting straight-line travel distance and other such approximations, then perhaps something like the following approach would suffice, fwiw:  (only intent is slight improvement of distribution versus a rigid set of two_pi/N angles, not "perfection")

we are given N and P, (or, at least, the OP knows those values), so we also know ideal P/N (aka "sites per inspector").  now let's also define a bit of a "fudge factor", say 5%? around that ideal value, and consider the value P/N-F.

now iterate your wedge angles one degree at a time (or some other small amount - there's a speed/accuracy trade-off impleid here) and keep track of the total accumulated angle.  create a "test" wedge at the current total accumulated angle, and count its contents.  if within tolerance of the fudge factor, accept it and reset your total;  if not, continue to next iteration and test again with a slightly larger wedge.

the last wedge processed will accumulate the fudge, and might not fall within tolerance.  but it should help indicate if your tolerance is too loose, or if you need smaller angular increments to be able to "divvy up" more finely, or any other such "tweaks".  an additional run after thoughtful tweaking ought to get you really close.

[DavidBollinger]

aside, fwiw (and replying to self):  we did exactly as described above to generate requested "what if radial" scenarios for 2010 Supervisorial Redistricting, using Census centroids as the sites and their summed population as the value (ie, rather than number of included sites, the value of all included sites).  we did it a couple of times using various "starting angles" (which occasionally had dramatic influence on final result), but they were all just "curiosities" - they all failed other conditions required of those particular boundaries, though it was pretty easy to satisfy the equal-within-tolerance distributions.

[JuanetteWillis]

David Bollinger, I am updating this distribution which was discussed back in 2016.  As you suggested, I am looking at using values instead of count. You mentioned you had done so for a 2010 Supervisorial Redistricting. If by some chance you happen to still have your script, I would love to have a look at it. I haven't touched python since 2016 and am just now working on my script and appreciate any assistance. I am interested in both the using values and adjusting the starting angles. 

The gang on this site were incredibly helpful in 2016! Xander Bakker‌, Your script worked beautifully with some personal tweaks. You have no idea how grateful I am to you. 

[XanderBakker]

Hi ojwillis , nice to hear that you have been able to tweak it and get some nice results... Thanks for posting this!