Oriented buffers

7888
28
Jump to solution
02-07-2015 07:58 PM
davesacco
New Contributor II

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

28 Replies
XanderBakker
Esri Esteemed Contributor

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
New Contributor II

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!

0 Kudos
ZackTucker
New Contributor II

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

0 Kudos
XanderBakker
Esri Esteemed Contributor

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?

0 Kudos
ZackTucker
New Contributor II
OBJECTID_1StationIDShape_STArShape_STLeAngleX_PointY_PointPlus30Minus30DistanceMeters
7331598.63835E-050.04807956188.087351574539.1632241287.79118.087345158.0873451132186.9
7341598.63835E-050.04807956187.86751574539.1632241287.79117.867496657.8674966332186.9
7351598.63835E-050.0480795613.7429081574539.1632241287.7933.7429079333.742907932186.9
7361598.63835E-050.04807956187.953091574539.1632241287.79117.953087657.953087632186.9
7371598.63835E-050.0480795616.6562651574539.1632241287.7936.65626485336.656264832186.9
7381598.63835E-050.0480795617.6209671574539.1632241287.7937.62096675337.620966732186.9
7391598.63835E-050.048079561357.63681574539.1632241287.7927.63680285327.636802832186.9
7401598.63835E-050.048079561359.31381574539.1632241287.7929.31379995329.313832186.9
7411601.36592E-050.01530708419.554811623294.9072247645.70149.55480509349.554805132186.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.

0 Kudos
ZackTucker
New Contributor II

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

0 Kudos
XanderBakker
Esri Esteemed Contributor

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
New Contributor II

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
Esri Esteemed Contributor

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

0 Kudos
Hayley_DelMaynard
New Contributor II

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

0 Kudos