Oriented buffers

davesacco · ‎02-07-2015

Hi there,

I am using arcview 9.3.1 and I need to create oriented buffers from points in a shapefile. The euclidean direction method almost works, but because it calculates to the nearest point, buffers associated with points within the maximum will be truncated. I need each point to have a complete buffer in a specific direction, regardless of overlap with other points.

Any ideas would be greatly appreciated.

Thanks,

Dave

XanderBakker · ‎02-09-2015

I haven't tested this, but it should do the trick (I hope...):

#-------------------------------------------------------------------------------
# Name:       createSectors2.py
# Purpose:    create circle sectors for points
#
# Author:      Xander
#
# Created:    08-02-2015
#-------------------------------------------------------------------------------
import arcpy

def main():
    import os
    fc_in = r"C:\Forum\CircleSectors\data.gdb\points2"
    fc_out = r"C:\Forum\CircleSectors\data.gdb\sectors2"

    fld_start = "StartAngle"
    fld_end = "EndAngle"
    fld_dist = "Length"
    fld_oid = "OIDpoints" # output field with oid of points

    sr = arcpy.Describe(fc_in).spatialReference
    arcpy.env.overwriteOutput = True

    # create empty output featureclass
    fc_ws, fc_name = os.path.split(fc_out)
    arcpy.CreateFeatureclass_management(fc_ws, fc_name, "POLYGON", spatial_reference=sr)

    # add field to polygon fc to store oid of points
    arcpy.AddField_management(fc_out, fld_oid, "LONG")

    # start insert cursor on polygon fc
    flds_out = ("SHAPE@", fld_oid)
    with arcpy.da.InsertCursor(fc_out, flds_out) as curs_out:

        flds = ("SHAPE@", fld_start, fld_end, fld_dist, "OID@")
        lst_polygons = []
        with arcpy.da.SearchCursor(fc_in, flds) as curs:
            for row in curs:
                pnt_g = row[0]
                start = row[1]
                end = row[2]
                length = row[3]
                oid = row[4]
                circle = pnt_g.buffer(length)
                pnt = pnt_g.firstPoint
                arrPnts = arcpy.Array()
                arrPnts.add(pnt)
                for bearing in range(int(start), int(end) + 1):
                    arrPnts.add(createPointAtAngleWithBearing(pnt, bearing, length))
                polygon = arcpy.Polygon(arrPnts, sr)
                curs_out.insertRow((polygon, oid, ))

def createPointAtAngleWithBearing(pnt, angle, distance):
    import math
    angle = math.radians(angle)
    dist_x, dist_y = (distance * math.cos(angle), distance * math.sin(angle))
    return arcpy.Point(pnt.X + dist_x, pnt.Y + dist_y)

if __name__ == '__main__':
    main()

davesacco · ‎02-09-2015

I wish I could upgrade, but it's not in the budget. And, thanks for the offer to process the data. I'll decline as I still need to determine the angles to be used for each point and incorporate a surface roughness factor to scale the length. Once I have all the pieces together I should be able to find a machine that I can use for a bit.

As for the data extraction, I converted the bedrock data to a raster and then used the tabulate area tool to create a table that I can join back to the polygons.

Thanks again for your help Xander. I will follow your link on learning python and one day return the favour back to the community!

ZackTucker · ‎12-12-2019

Xander,

I have found that this code does not produce oriented sectors of the buffer when the angles cross the 0/360 degrees axis (due East). Do you have any idea why?

Example: If I have Angles 5 and 335 with a distance of 20 miles, no buffer is created, and I have found this to be true of 830+ of 4,500 records.

I am trying to use 60 degree segments, so anything in between the ranges from 59 - 359 degrees & 301 - 1 degrees does not compute in the above script. Instead the shape length and Area of the polygon is returned as 0.

Is the issue stemming from an exception needes in the "For bearing In range" function and parameters?

I'm running ArcMap 10.4

If you can see this still, Thank you.

-ZT

XanderBakker · ‎12-12-2019

Hi zmotuck ,

This is a really old thread... brings back some memories... 🙂

Would you have any data that you can share to be able to dig into the code and test and adjust it?

ZackTucker · ‎12-12-2019

OBJECTID_1	StationID	Shape_STAr	Shape_STLe	Angle	X_Point	Y_Point	Plus30	Minus30	DistanceMeters
733	159	8.63835E-05	0.048079561	88.08735	1574539.163	2241287.79	118.0873451	58.08734511	32186.9
734	159	8.63835E-05	0.048079561	87.8675	1574539.163	2241287.79	117.8674966	57.86749663	32186.9
735	159	8.63835E-05	0.048079561	3.742908	1574539.163	2241287.79	33.7429079	333.7429079	32186.9
736	159	8.63835E-05	0.048079561	87.95309	1574539.163	2241287.79	117.9530876	57.9530876	32186.9
737	159	8.63835E-05	0.048079561	6.656265	1574539.163	2241287.79	36.65626485	336.6562648	32186.9
738	159	8.63835E-05	0.048079561	7.620967	1574539.163	2241287.79	37.62096675	337.6209667	32186.9
739	159	8.63835E-05	0.048079561	357.6368	1574539.163	2241287.79	27.63680285	327.6368028	32186.9
740	159	8.63835E-05	0.048079561	359.3138	1574539.163	2241287.79	29.31379995	329.3138	32186.9
741	160	1.36592E-05	0.015307084	19.55481	1623294.907	2247645.701	49.55480509	349.5548051	32186.9

Here is a small sample that I have ran with your code.

The spatial reference is:

You can make those points an XY event layer and then run your script as is. (changing the field names of course).

I use plus30 and minus30 columns as my begin and end angles (that should be 30 +/- from the 'Angle' Column). The distance is in meters (20 miles).

Three of the oriented buffers appear when I run that data, as referenced below, While the seven others fall within a range of angles that cross the 0/360 degree line.

ZackTucker · ‎12-12-2019

Sorry for such a small sample, I can send a larger one with other values in the other quadrants of the Arithmetic Planes.

However, I am new to this forum and finding it not too easy to upload that here. How does one provide shapefiles or excel data sheets?

Regards,

ZT

XanderBakker · ‎12-13-2019

Hi Zack Tucker ,

Good catch. The code simply will create a point every degree using the range from Minus30 to Plus30 in your case and if Minus30 is larger than Plus30 no points will be created. I added a small fix at lines 48-49 and it should work now:

#-------------------------------------------------------------------------------
# Name:       createSectors3.py
# Purpose:    create circle sectors for points
#
# Author:      Xander
#
# Created:    13-12-2019
#-------------------------------------------------------------------------------
import arcpy

def main():
    import os
    fc_in = r"C:\GeoNet\OrientedBuffers\data.gdb\points"
    fc_out = r"C:\GeoNet\OrientedBuffers\data.gdb\sectors_v03"

    fld_start = "Minus30" # "StartAngle"
    fld_end = "Plus30" # "EndAngle"
    fld_dist = "DistanceMeters" # "Length"
    fld_oid = "OIDpoints" # output field with oid of points

    sr = arcpy.Describe(fc_in).spatialReference
    arcpy.env.overwriteOutput = True

    # create empty output featureclass
    fc_ws, fc_name = os.path.split(fc_out)
    arcpy.CreateFeatureclass_management(fc_ws, fc_name, "POLYGON", spatial_reference=sr)

    # add field to polygon fc to store oid of points
    arcpy.AddField_management(fc_out, fld_oid, "LONG")

    # start insert cursor on polygon fc
    flds_out = ("SHAPE@", fld_oid)
    with arcpy.da.InsertCursor(fc_out, flds_out) as curs_out:

        flds = ("SHAPE@", fld_start, fld_end, fld_dist, "OID@")
        lst_polygons = []
        with arcpy.da.SearchCursor(fc_in, flds) as curs:
            for row in curs:
                pnt_g = row[0]
                start = row[1]
                end = row[2]
                length = row[3]
                oid = row[4]
                circle = pnt_g.buffer(length)
                pnt = pnt_g.firstPoint
                arrPnts = arcpy.Array()
                arrPnts.add(pnt)
                if end < start:
                    end += 360.0
                for bearing in range(int(start), int(end) + 1):
                    arrPnts.add(createPointAtAngleWithBearing(pnt, bearing, length))
                polygon = arcpy.Polygon(arrPnts, sr)
                curs_out.insertRow((polygon, oid, ))

def createPointAtAngleWithBearing(pnt, angle, distance):
    import math
    angle = math.radians(angle)
    dist_x, dist_y = (distance * math.cos(angle), distance * math.sin(angle))
    return arcpy.Point(pnt.X + dist_x, pnt.Y + dist_y)

if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Result (after changing some of the overlapping point locations):

ZackTucker · ‎12-13-2019

Xander,

The code edit worked!!! This helped indeed. All sectors can now be represented and created into buffer sections, regardless of if they pass due East (0/360 degree axis)

Thank you for helping me understand the For loop and Point creation process to solve the problem. I think this geoprocessing is amazing, especially in regards to linear features and directions of flow or travel. This oriented buffers page seems to be the only resource out there to do something like this. This method should be added in the standard buffer tool when using vector point data.

Thank you so much for your quick response!

XanderBakker · ‎12-13-2019

You're welcome, glad it works and thanks for catching the error!

Hayley_DelMaynard · ‎02-18-2017

Xander Bakker‌ Thanks so much for this! I'm trying to create 4 90degree buffers around all points in my data frame. I tried using your code and I keep getting this below error. I'm pretty new to Python, do you know how to fix this? Thanks again!

>>> import arcpy
...
... def main():
... fc_in = "C:\Users\A0701216\Documents\ArcGIS\Default.gdb\ALEX"
... fc_out = "C:\Users\A0701216\Documents\ArcGIS\Default.gdb\ALEX_"
...
... fld_start = 180
... fld_end = 225
... fld_dist = 15
...
... sr = arcpy.Describe(fc_in).spatialReference
... arcpy.env.overwriteOutput = True
...
... flds = ("SHAPE@", fld_start, fld_end, fld_dist)
... lst_polygons = [1]
... with arcpy.da.SearchCursor(fc_in, flds) as curs:
... for row in curs:
... pnt_g = row[0]
... start = row[1]
... end = row[2]
... length = row[3]
... circle = pnt_g.buffer(length)
... pnt = pnt_g.firstPoint
... arrPnts = arcpy.Array()
... arrPnts.add(pnt)
... for bearing in range(int(start), int(end) + 1):
... arrPnts.add(createPointAtAngleWithBearing(pnt, bearing, length))
... polygon = arcpy.Polygon(arrPnts, sr)
... lst_polygons.append(polygon)
...
... arcpy.CopyFeatures_management(lst_polygons, fc_out)
...
... def createPointAtAngleWithBearing(pnt, angle, distance):
... import math
... angle = math.radians(angle)
... dist_x, dist_y = (distance * math.cos(angle), distance * math.sin(angle))
... return arcpy.Point(pnt.X + dist_x, pnt.Y + dist_y)
...
... if __name__ == '__main__':
... main()
...
Runtime error
Traceback (most recent call last):
File "<string>", line 40, in <module>
File "<string>", line 16, in main
TypeError: 'field_names' must be string or non empty sequence of strings
>>>