Script for numbering roads based on Cardinal Direction

Question

I am kind of new to python, and am attempting to write a script that will give me a unique identifier for roads based on the cardinal direction. Ultimatley, I would like to have a system set up such as the interstate system in the united states, where North -South roads are odd numbered, and east west roads are even numbered. I would also like to reflect these numbers based on their position within the area I am working in. For example, the southern most w to e road would take the number 2, while the western most n-s road would take 1. So far, the only code I could find is code to conver the endpoints of the roads into a numerical direction:

180 + math.atan2((!Shape.lastpoint.Y! - !Shape.firstpoint.Y!),(!Shape.lastpoint.X! - !Shape.firstpoint.X!)) * (180 / math.pi)

I am assuming I will have to write some if/else statements, but I am not entirely sure how to set them up.

I also found this block of code, but it didnt seem to work for me for numbering the n-s roads:

rec=0 
def autoIncrement(): 
  global rec 
  pStart = 1 #adjust startvalue, if req'd
  pInterval = 2 #adjust interval value, if req'd
  if (rec == 0): 
    rec = pStart
  else: rec = rec + pInterval
    return rec

Even if I were to get this to work, I still don't know how to get it to number roads based on how far north/south, east/west they are.

Answer 1

Update you roads direction using a formula something like

if (45 <= (360 + math.atan2((!Shape.lastpoint.X! - !Shape.firstpoint.X!),(!Shape.lastpoint.Y! - !Shape.firstpoint.Y!)) * (180 / math.pi)) % 360 % 180) < 135):
    dir = 'W-E' 
else:
    dir = 'N-S'

Then for each direction group number them in order using the first opposite ordinate to sort. For example for the north south roads number in order based on the first X. West east roads would be ordered by the first Y.

I'm not a python scripter, but as a database query it would be along the lines of

SELECT ID, Direction, CalcDirection,
    (ROW_NUMBER() OVER (
        PARTITION BY CalcDirection 
        ORDER BY CASE CalcDirection 
            WHEN 'NS' THEN Geom.STPointN(1).STX 
            ELSE Geom.STPointN(1).STY 
            END
            ) 
        * 2) - 
    CASE CalcDirection WHEN 'WE' THEN 0 ELSE 1 END CalcID
FROM (
    SELECT ID, Direction, Geom,
        CASE WHEN
            CAST(360 + ROUND(ATN2((Geom.STPointN(Geom.STNumPoints()).STX - Geom.STPointN(1).STX),(Geom.STPointN(Geom.STNumPoints()).STY - Geom.STPointN(1).STY)) * (180.0 / PI()),0) AS INT) % 360 % 180
             BETWEEN 45 and 134 THEN 'WE' ELSE 'NS' END CalcDirection
    FROM Lines
    ) a

Using the following geometries:

    (1, 'WE', Geometry::STGeomFromText('LINESTRING (0 0, 20 5)',0))
    ,(2, 'WE', Geometry::STGeomFromText('LINESTRING (0 10, 20 0)',0))
    ,(3, 'WE', Geometry::STGeomFromText('LINESTRING (0 20, 20 15)',0))
    ,(4, 'NS', Geometry::STGeomFromText('LINESTRING (0 0, 5 20)',0))
    ,(5, 'NS', Geometry::STGeomFromText('LINESTRING (10 0, 10 20)',0))
    ,(6, 'NS', Geometry::STGeomFromText('LINESTRING (15 0, 0 20)',0))
    ,(7, 'WENS', Geometry::STGeomFromText('LINESTRING (0 0, 20 20)',0))
    ,(8, 'NS', Geometry::STGeomFromText('LINESTRING (5 20, 0 0)',0))
    ,(9, 'NS', Geometry::STGeomFromText('LINESTRING (10 20, 10 0)',0))
    ,(10, 'NS', Geometry::STGeomFromText('LINESTRING (0 20, 15 0)',0))
    ,(11, 'NSWE', Geometry::STGeomFromText('LINESTRING (0 20, 20 0)',0))
    ,(12, 'WE', Geometry::STGeomFromText('LINESTRING (20 5, 0 0)',0))
    ,(13, 'WE', Geometry::STGeomFromText('LINESTRING (20 0, 0 10)',0))
    ,(14, 'WE', Geometry::STGeomFromText('LINESTRING (20 15, 0 20)',0))

I get the following results:

ID          Direction CalcDirection CalcID
----------- --------- ------------- --------------------
4           NS        NS            1
11          NSWE      NS            3
10          NS        NS            5
8           NS        NS            7
9           NS        NS            9
5           NS        NS            11
6           NS        NS            13
1           WE        WE            2
7           WENS      WE            4
13          WE        WE            6
12          WE        WE            8
2           WE        WE            10
14          WE        WE            12
3           WE        WE            14

Answer 2

Here is a decent start to what you're asking. It can be more accurate if you somehow gathered the "mean direction" of all the vertices and used that to determine the azimuth for a given road segment. Also, this uses an arcpy.Sort_management function that sorts the feature table spatially starting, in this case, with a particular corner of the data (like "LR" for lower right).

import arcpy
import math
### Change these variables
shp = "original.shp"
newshp1 = "copy.shp"

arcpy.env.overwriteOutput=True

sort_fields = [["Shape","ASCENDING"]]
## Sort fields from Lower Right of data so the iteration grabs the first ones in the southeast first and works to northwest
arcpy.Sort_management(shp, newshp1, sort_fields, "LR")

##The two input fields to be updated and used for naming "NewNum" is resulting street number
arcpy.AddField_management(newshp1, 'NewNum','DOUBLE')
arcpy.AddField_management(newshp1, 'Cardinal','TEXT')

## Since we add 2 everytime, it's necessary to start at -1 for the odd numbered streets
oddNum = -1
evenNum = 0

## Must be UpdateCursor to update the row on copied shp
with arcpy.da.UpdateCursor(newshp1, ["OID@","SHAPE@","NewNum","Cardinal"]) as cursor:
    for row in cursor:
        print("Feature {0}:".format(row[0]))
        startpt = row[1].firstPoint
        endpt = row[1].lastPoint
        azi = 180 + math.atan2(( startpt.X - endpt.X ),(startpt.Y - endpt.Y ) ) * (180 / math.pi)
        if azi >= 315 or azi <= 45:
            card = "N"
            evenNum+=2
            row[2] = evenNum
        elif azi < 315 and azi >= 225:
            card = "W"
            oddNum+=2
            row[2] = oddNum
        elif azi < 225 and azi >= 135:
            card = "S"
            evenNum+=2
            row[2] = evenNum
        else:
            card = "E"
            oddNum+=2
            row[2] = oddNum
        row[3] = card
        cursor.updateRow(row)
        print card

An even better way might be to sort based on latitude and longitude values of the start and/or end point where the southern-most road would have the lowest latitudinal value (or highest in the southern hemisphere) and the eastern-most road would have the highest longitudinal value.