After many different approaches. I have found a series of processes that seem to work rather well to generate the buffers I need. The solution I came up with uses ArcGIS, heavily relying on ArcPy and some basic pre-calculus level algebra.
For starters I assume that any given line receiving a buffer is a straight line. This is not always the case with the data in question but overall my margin of error was only .35% at it's highest. the more a line curves however, the more this margin of error will go up as a buffer starts "folding back" onto itself. So... Considering the following image:
We know D (the length of the line feature) and we also know the total area of the buffer (again we assume straight lines).
We can then turn this into a quadratic formula to solve for C:
A = D + 2C
B = 2C
A * B = AreaOfBuffer(AoB)
AoB = 2C(D + 2C)
Quadratic:
4C^2 + 2CD - AoB = 0
Using the quadratic formula we can solve for C however, as with all quadratic functions we get two answers. In this case we will always have a negative answer(we can ignore this one) and a positive answer(this is C).
After acquiring C buffers can be generated for each line that with the necessary dimensions. However, a standard buffer is rounded on each end and if we were to change the buffer ends parameter to "flat", we end up with a buffer that contains less area than desired because the buffer ends at the line beginning/end. To account for this I was able to find a script on SE that extends lines by an amount. The author of the script knows their python. This is one of the more pythonic spatial scripts I've ever seen. I was able to adapt it to take a field as the extension distance as apposed to a constant value for all features. This was really the only change I made to the line extending script. Again, I did not author this one, credit goes to @Paul
from math import hypot
import collections
from operator import add
import arcpy
layer = arcpy.GetParameterAsText(0)
distanceField = arcpy.GetParameterAsText(1)
#Computes new coordinates x3,y3 at a specified distance
#along the prolongation of the line from x1,y1 to x2,y2
def newcoord(coords, dist):
(x1,y1),(x2,y2) = coords
dx = x2 - x1
dy = y2 - y1
linelen = hypot(dx, dy)
x3 = x2 + dx/linelen * dist
y3 = y2 + dy/linelen * dist
return x3, y3
#accumulate([1,2,3,4,5]) --> 1 3 6 10 15
#Equivalent to itertools.accumulate() which isn't present in Python 2.7
def accumulate(iterable):
it = iter(iterable)
total = next(it)
yield total
for element in it:
total = add(total, element)
yield total
#OID is needed to determine how to break up flat list of data by feature.
coordinates = [[row[0], row[1]] for row in
arcpy.da.SearchCursor(layer, ["OID@", "SHAPE@XY"], explode_to_points=True)]
distances = []
with arcpy.da.SearchCursor(layer, distanceField) as cursor:
for row in cursor:
distances.append(row[0])
oid,vert = zip(*coordinates)
#Construct list of numbers that mark the start of a new feature class.
#This is created by counting OIDS and then accumulating the values.
vertcounts = list(accumulate(collections.Counter(oid).values()))
#distCol = list(accumulate(collections.Counter(dist).values()))
#Grab the last two vertices of each feature
lastpoint = [point for x,point in enumerate(vert) if x+1 in vertcounts or x+2 in vertcounts]
for pt in lastpoint:
arcpy.AddMessage(str(pt))
for d in distances:
arcpy.AddMessage(str(d))
#Convert flat list of tuples to list of lists of tuples.
#Obtain list of tuples of new end coordinates.
newvert = [newcoord(y, distances[lastpoint.index(y[1]) / 2]) for y in zip(*[iter(lastpoint)]*2)]
j = 0
with arcpy.da.UpdateCursor(layer, "SHAPE@XY", explode_to_points=True) as rows:
for i,row in enumerate(rows):
if i+1 in vertcounts:
row[0] = newvert[j]
j+=1
rows.updateRow(row)
I ran this script on my line data and then flipped each line and ran it again. This created a line that ran the extent of a normal buffer. I could then run the buffer tool again with "flat ends" selected. The final buffers are 99.99% accurate when run on straight lines. They do vary slightly with many curves, but even then the highest margin of error seen was .35%
Here is the quadratic formula solver script:
import arcpy, math
#Length of line and area need to be the same base unit I.E. (ft and sqft)
inFC = arcpy.GetParameterAsText(0) #Input line features
retSQFT = arcpy.GetParameterAsText(1) #Desired Buffer Area field
lenFT = arcpy.GetParameterAsText(2) # Length Of Line field
bufDist = arcpy.GetParameterAsText(3) #Blank field to populate with buffer distance
#List of all fields needed
fields = [bufDist, retSQFT, lenFT]
#feature counter
index = 0
#Iterate through features and solve quadratic formula based on inputs.
with arcpy.da.UpdateCursor(inFC, fields) as cursor:
for row in cursor:
quadSol1 = "NA"
quadSol2 = "NA"
#Assign quadratic variables
a = 4
b = row[2] * 2
c = (-1)*(row[1])
arcpy.AddMessage(" ")
#Positive version of the quadratic
try:
quadSol1 = (-(b) + math.sqrt((b * b) - (4 * (a) * (c)))) / ((2) * (a))
except:
#In case we try to divide by zero, not sure we would ever need this
arcpy.AddMessage("Index: " + str(index) + " -----> real number exception in addition quadratic")
#Negative version of the quadratic
try:
quadSol2 = (-(b) - math.sqrt((b * b) - (4 * (a) * (c)))) / ((2) * (a))
except:
#In case we try to divide by zero, not sure we would ever need this
arcpy.AddMessage("Index: " + str(index) + " -----> real number exception in subtraction quadratic")
arcpy.AddMessage("Quadratic Solution 1: " + str(quadSol1))
arcpy.AddMessage("Quadratic Solution 2: " + str(quadSol2))
arcpy.AddMessage(" ")
#Get whichever solution comes out as a positive real number
if(quadSol1 != "NA"):
if(quadSol1 > 0):
buffDist = quadSol1
else:
if(quadSol2 != "NA" and quadSol2 > 0):
buffDist = quadSol2
#populate the buffer distance field
row[0] = buffDist
#feature counter up
index += 1
#Update the row
cursor.updateRow(row)
I plan on combining these two scripts together and adding in the line flipping tools, etc. in order to create one seamless tool to run for the scenario.
Here is an example of the results.
