1

I was recently asked to create a very odd "buffer" around several line features. I'm not sure if this is innately possible with an ArcGIS Advanced license or not and therefore am willing to explore other solutions such as QGIS. I am proficient in ArcPy, so please include any ArcGIS scripting solutions as well.

I have been given a set of lines that have a field called "Acreage". The ultimate goal would be to generate polygons around each line that have the area of the "Acreage" field. So, when considering the available data... I have the length of each line and the area that each final polygon should be. This alone doesn't seem like it should be to difficult to come up with a formula to calculate the buffer distance, especially if we were to assume each line is straight (though they are not straight, this might be the best I could do). However, to make matters even more complicated those who asked for this would like to see the final result as "bent" rectangles. The image below shows a mock-up of what they would like to see.

enter image description here

In the Image the red line is the starting data and the green "buffer" is what needs to be generated from each line.

The key elements are that the "buffer" needs to have the designated area(one of the line attributes) and also needs to have 90 degree corners.

If getting the corners is much more difficult or not possible. I would also be happy with a solution that just creates a "normal" rounded buffer that contains the right area.

  • Do the corners at the bend need to be angles as well? Because most buffer algorithms work on the nearest-distance principle and so produce curves at bends in straight lines.... – Spacedman Aug 31 '16 at 16:50
  • You can create buffers like this in R using rgeos::gBuffer if you specify join="mitre", capStyle="square" and a large mitreLimit. If those parameters are available for QGIS or Arc's buffering then it might help you get the right shape, if not the right size (which you can get in R with a loop and a search over the width...). – Spacedman Aug 31 '16 at 17:01
  • @Spacedman It would be fine to have the inner angles be curves. I really just drew it like that kind of hastily. – GeoJohn Aug 31 '16 at 19:40
  • @Spacedman Interesting idea using R. I might have to look into that, but unfortunately it is not a language I'm all that familiar with. – GeoJohn Aug 31 '16 at 19:41
1

Input:

enter image description here

Output:

enter image description here

Algorithm:

  • Extend line on both ends by certain amount and do flat ends buffer
  • Repeat by changing extension length until precision required met

Script:

import arcpy
from arcpy import env
env.overwriteOutput = True
singleLine=r"in_memory\line"
singleBuffer=r"in_memory\pgon"
mxd = arcpy.mapping.MapDocument("CURRENT")
lines = arcpy.mapping.ListLayers(mxd,"LINES")[0]
buffers = arcpy.mapping.ListLayers(mxd,"BUFFERS")[0]
curT=arcpy.da.InsertCursor(buffers,"SHAPE@")
with arcpy.da.SearchCursor(lines,["SHAPE@","AREA_M2"]) as cursor:
    for shp,target in cursor:
        part=shp.getPart(0);origList=list(part)
##      get last vertices vector
        p1,p2=origList[-2],origList[-1]
        dX=p2.X-p1.X; dY=p2.Y-p1.Y; d1=pow(dX*dX+dY*dY,0.5)
        dXend=dX/d1;dYend=dY/d1
##      get first vertices vector
        p1,p2=origList[1],origList[0]
        dX=p2.X-p1.X; dY=p2.Y-p1.Y; d2=pow(dX*dX+dY*dY,0.5)
        dXstart=dX/d2;dYstart=dY/d2
##      define bounds
        L=target/shp.length/2
        low=L/2;high=L*2
        while True:
            middle=0.5*(low+high)
##      extend line        
            p1=arcpy.Point(origList[-1].X+middle*dXend,origList[-1].Y+middle*dYend)
            p2=arcpy.Point(origList[0].X+middle*dXstart,origList[0].Y+middle*dYstart)
            extList=[p2]+origList+[p1]
            extLine=arcpy.Polyline(arcpy.Array(extList))
            arcpy.CopyFeatures_management(extLine, singleLine)
##      get buffer ind repeat        
            arcpy.Buffer_analysis(singleLine, singleBuffer, middle,"FULL","FLAT","NONE","","PLANAR")            
            with arcpy.da.SearchCursor(singleBuffer,"SHAPE@") as inMem:
                for r in inMem:
                    pgon=r[0]
#       tolerance                    
            if (high-low)<0.01: break
            curArea=pgon.area
            if curArea<target:low=middle
            else:high=middle
        arcpy.AddMessage("Difference {:8.2f}%".format((curArea-target)/target*100))
        curT.insertRow((pgon,))
arcpy.Delete_management("in_memory")

It is good you are proficient with arcpy, you'll figure out parameters

UPDATE:

Your approach gives the impression of working simply because your lines are not that bendy and most importantly buffers are skinny. Try with the shape below and buffering distance comparable to line length and you’ll notice significant increase in mismatch with target area.

enter image description here

The class of problems you are dealing with has no analytical solution. However an accurate estimate can be found numerically through iterations, tries and errors but much quicker by using one of well known root-finding techniques. I used the method of bisection that is not very efficient. The golden section or others might converge faster. Nevertheless in took only 18 iterations to find buffer distance, accurate to 0.01 m for 12 km long shape with 1000 vertices shown above. As for it’s accuracy, gauge by yourself, I am struggling to express it using percents.

It took 5 seconds only on my out of date machine at home. Unfortunately buffer() method of arcpy geometry does not support flat ends ( this is why I had to use tools working in_memory), otherwise it could be done in no time.

  • This is excellent work! Definitely better for accuracy! Unfortunately I was dealing with a situation that desired time sensitivity over accuracy. However, in the future if this ever comes up again your solution would definitely reign superior as long as time isn't a factor! – GeoJohn Sep 2 '16 at 18:22
2

I was looking into this exact same issue last night and I am exploring using the: toolboxes\system toolboxes\analysis tools.tbx\proximity\multiple ring buffer To try and generate multiple buffers and hopefully giving me a start point for a buffer that is flat on the end. I know its a long shot but it may get me in the ball park allowing just some tweaks to get the final product. This is the Script I found on Arc:

# Name: MultipleRingBuffer_Example2.py
# Description: Create multiple buffers for the input features

# Import system modules
import arcpy

# Set environment settings
arcpy.env.workspace = "C:/data/airport.gdb"

# Set local variables
inFeatures = "schools"
outFeatureClass = "c:/output/output.gdb/multibuffer1"
distances = [10, 20, 30]
bufferUnit = "meters"

# Execute MultipleRingBuffer
arcpy.MultipleRingBuffer_analysis(inFeatures, outFeatureClass, distances, bufferUnit, "", "ALL")

Found at the link below:

http://pro.arcgis.com/en/pro-app/tool-reference/analysis/multiple-ring-buffer.htm

2

You could try one of those solutions (in QGIS):

Buffer with a surface

https://github.com/jdugge/BufferByPercentage

Buffer by Percentage

Buffer polygon features so the buffered area is a specified percentage of the original area. Instead of buffering a polygon using a specified buffer distance, this plugin lets the user specify the area the buffered polygon should cover, as a percentage of the original polygon's area.

Author: Juernjakob Dugge

Available version: 0.2.4.2

Buffer with square ends

https://grass.osgeo.org/grass73/manuals/v.buffer.html

GRASS - v.buffer

v.buffer creates a buffer around features of given type, which have a category in the given layer. The tolerance controls the number of vector segments being generated (the smaller the value, the more vector segments are generated). Straight corners with caps are created by -s flag (red color on the figure below), while -c flag doesn't make caps at the ends of polylines (green color on the figure below):

  • I don't think BufferByPercentage accepts linestring features – LeasMaps Jan 2 at 8:28
0

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:

enter image description here

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.

enter image description here

  • I don't see it working, result are different from target at 3rd signigicant number – FelixIP Sep 2 '16 at 2:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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