Saturday, 11 June 2016

qgis - Creating buffers of a specific size and shape


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.



Answer



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 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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...