There may be slight variations of a single approach that can be used. Let me elaborate that single approach step-wise.
1. Define or calculate intersection of roads
To split the road intersections into quadrants, intersections must be 4-way intersections only. Better if you already have only 4-way intersections point shapefile. Otherwise, you may need to calculate the roads intersections and delete each intersection that is not 4-way. This'll require additional scripting.
2. Generate 10 m buffers around intersections
arcpy.Buffer_analysis can be used to do the task.
3. Split buffers by "intersecting roads" into quadrants
This can be done by using a sequence of tools. Let me list down the sequence.
- arcpy.FeatureToLine_management (Convert bounding geometry to lines)
- arcpy.Merge_management (Merge lines of bounding geometry and intersecting roads)
- arcpy.ExtendLine_edit (if required)
- arcpy.FeatureToPolygon_management (Convert merged lines into polygons)
- arcpy.Identity_analysis (To finally divide buffers into quadrants)
4. Convert quadrants & original buffers into point for angle calculation
This step can be performed using arcpy.FeatureToPoint_management. Converting original buffers into points isn't actually needed. Original intersection may be used. This was included initially to ease clean-up at the end of whole process.
5. Calculate angles between 2 point shapefiles
6. Join the point shapefile (angles calculated) with the polygon shapefile (quadrants)
7. Delete all extra fields, add new field to hold intercardinal directions and calculate directions
8. Perform final clean-up
arcpy.Delete_management (Delete all intermediate layers)
Arcpy script of the whole process is given below.
Note:- Please note that there must be at least 4 intersecting lines. If there are more, for example, a fifth line bisecting NE quadrant into two, both quadrants will still be marked as NE
# import system modules
import arcpy, os
rdshp = <path to shapefile> # Input Road Shapefile
ptshp = <path to shapefile> # Input Road Intersection Shapefile (Optional)
qdshp = <path to shapefile> # Output Quadrant Shapefile
path = os.path.split(rdshp)
# Check if intersection shapefile is defined
# Create one if not defined
inter = ptshp
inter = os.path.join(path,"ip.shp")
clusterTolerance = 0
arcpy.Intersect_analysis(rdshp, inter, "", clusterTolerance, "point")
buffs = os.path.join(path,"bf.shp") # Define: 10 m Buffer shapefile
qushp = os.path.join(path,"iq.shp") # Define: Shapefile of buffers splitted into quadrants
qptss = os.path.join(path,"qp.shp") # Define: Convert each quadrant to point for angle calculation
bufpt = os.path.join(path,"bp.shp") # Define: 10 m Buffer shapefile back into points
# Process: Generate 10 meters buffer around intersections
distanceField = "10 Meter"
sideType = "FULL"
endType = "ROUND"
dissolveType = "NONE"
arcpy.Buffer_analysis(inter, buffs, distanceField, sideType, endType, dissolveType)
# Process: Split buffers by intersecting roads into quadrants
Geometry_Type = "ENVELOPE"
if Geometry_Type == '#' or not Geometry_Type:
Geometry_Type = "ENVELOPE" # provide a default value if unspecified
Extend_Lines_Limit = 0
Extend_to_Extensions = "true"
if Extend_to_Extensions == '#' or not Extend_to_Extensions:
Extend_to_Extensions = "true" # provide a default value if unspecified
# Local variables:
tmpMinimumBoundingGeometry = "in_memory\\tmpMinimumBoundingGeometry"
tmpFeatureToLine = "in_memory\\tmpFeatureToLine"
tmpMergedLines = "in_memory\\tmpMergedLines"
ExtendedLines = tmpMergedLines
tmpNewPolygon = "in_memory\\tmpNewPolygon"
inmerge = [tmpFeatureToLine, rdshp]
# Process: Minimum Bounding Geometry
arcpy.MinimumBoundingGeometry_management(buffs, tmpMinimumBoundingGeometry, Geometry_Type, "NONE", "", "NO_MBG_FIELDS")
# Process: Feature To Line
arcpy.FeatureToLine_management("in_memory\\tmpMinimumBoundingGeometry", tmpFeatureToLine, "", "NO_ATTRIBUTES")
# Process: Merge
arcpy.Merge_management(inmerge, tmpMergedLines, "")
# Process: Extend Line
arcpy.ExtendLine_edit(tmpMergedLines, Extend_Lines_Limit, Extend_to_Extensions)
# Process: Feature To Polygon
arcpy.FeatureToPolygon_management("in_memory\\tmpMergedLines", tmpNewPolygon, "", "ATTRIBUTES", "")
# Process: Identity
arcpy.Identity_analysis(buffs, tmpNewPolygon, qushp, "ONLY_FID", "", "NO_RELATIONSHIPS")
# Convert splitted polygons to points
arcpy.FeatureToPoint_management(qushp, qptss, "INSIDE")
# Convert 10 m buffer to point
arcpy.FeatureToPoint_management(buffs, bufpt, "INSIDE")
# Calculate angles in qptss
# find features only within search radius
search_radius = "15 Meters"
# Avoid finding locations of nearest features
location = "NO_LOCATION"
# Get angles of nearest features
angle = "ANGLE"
arcpy.Near_analysis(qptss, bufpt, search_radius, location, angle)
# Create a feature layer to perform join
# Join polygons FID to its converted points
arcpy.AddJoin_management("qulyr", "FID", qptss, "ORIG_FID")
# Perform final shapefile table clean-up by deleting unnecessary fields
fieldList = arcpy.ListFields(qdshp)
afield = "qp_NEAR_AN"
bfield = "qp_NEAR_FI"
nfield = "iDir"
for field in fieldList:
if not field.required and field.name != afield and field.name != bfield:
arcpy.AddMessage(str(field.name) + " failed")
# Create & calculate new field to hold intercardinal directions of quadrants
exp = "getClass(float(!qp_NEAR_AN!))"
codeblock = """
if afield <= 0 and afield >=-90:
if afield < -90 and afield >-180:
if afield >= 0 and afield <=90:
if afield > 90 and afield <180:
arcpy.AddField_management(qdshp, nfield, "TEXT", 2, "", "", "Directions", "NULLABLE", "REQUIRED")
arcpy.CalculateField_management(qdshp, nfield, exp, "PYTHON_9.3", codeblock)
# Perform clean-up by deleting all intermediate shapefiles