I have a road polyline feature class in an Oracle 18c/10.7.1 geodatabase.

I would like to generate a polygon FC that represents the quadrants of the road intersections.

enter image description here

  • The polygon FC would have a QUADRANT column.
  • Each polygon would have a value that represents the approximate intercardinal direction of the polygon (NE, SE, SW, NW).
  • The road polylines are angled at -19 degrees.

Is there a way to do this using ArcGIS Desktop 10.7.1 (advanced license)?

  • 1
    Should be pretty easy, if you have the intersection points as points and road lines that intersect at those points. Probably no more than two dozen lines of ArcPy, certainly less than two score for initial functionality. What do you have so far, and what problem have you encountered? – Vince Jan 18 at 19:18
  • You had a tag for ArcPy but have not included a code snippet. I’ve removed that rather than voting to close your question. If you start to code something then feel free to put it back. You can add your QUADRANT column using the Add Field tool. – PolyGeo Jan 18 at 20:54

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.

  1. arcpy.MinimumBoundingGeometry_management
  2. arcpy.FeatureToLine_management (Convert bounding geometry to lines)
  3. arcpy.Merge_management (Merge lines of bounding geometry and intersecting roads)
  4. arcpy.ExtendLine_edit (if required)
  5. arcpy.FeatureToPolygon_management (Convert merged lines into polygons)
  6. 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

arcpy.DeleteField_management arcpy.AddField_management arcpy.CalculateField_management

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)[0]

# Check if intersection shapefile is defined
# Create one if not defined
if ptshp:
    inter = ptshp
    inter = os.path.join(path,"ip.shp")
    clusterTolerance = 0
    arcpy.Intersect_analysis(rdshp, inter, "", clusterTolerance, "point")
    arcpy.DeleteIdentical_management(inter, "Shape")

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
arcpy.MakeFeatureLayer_management(qushp, "qulyr")

# Join polygons FID to its converted points
arcpy.AddJoin_management("qulyr", "FID", qptss, "ORIG_FID")
arcpy.CopyFeatures_management("qulyr", qdshp)

arcpy.DeleteIdentical_management(qdshp, "Shape")

# 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.DeleteField_management(qdshp, field.name)
            arcpy.AddMessage(str(field.name) + " failed")

# Create & calculate new field to hold intercardinal directions of quadrants
exp = "getClass(float(!qp_NEAR_AN!))"

codeblock = """
def getClass(afield):
    if afield <= 0 and afield >=-90:
        return "NW"
    if afield < -90 and afield >-180:
        return "NE"
    if afield >= 0 and afield <=90:
        return "SW"
    if afield > 90 and afield <180:
        return "SE"

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

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.