Street Centerline "From and To Streets"

1849
28
Jump to solution
06-20-2017 11:52 AM
KevinCross2
New Contributor III

I am trying to update our street center line information that does not have the "From Street" or "To Street" field populated.  All our data has is "On Street".   Is there an "easy" way to write a script or is there a set of tools that can be used to help populate this information automatically?  

0 Kudos
1 Solution

Accepted Solutions
XanderBakker
Esri Esteemed Contributor

Hi Kevin Cross  and Joe Borgione ,

So here goes... As I mentioned before I used the MultiPart to SinglePart tool to create a featureclass without multipart features. The reason for this is that there can be several line parts in a single feature, and when I extract the firstPoint and lastPoint of the Polyline it will not represent the start and end of the street.

After doing that I noticed that there were some cases (actually quite a number) where the same street continuous with another feature. Dissolve with the unsplit and without creating multiparts could help but this would make you loose some attributes that might be relevant.

So, I ordered the features on the FullName and checked the cases where in an edit session I could merge the feature into one feature. There are some very challenging cases in your data Kevin.

Some of them will simply have not result:

In this case Webster Ln in the middle would have from and to street with the same name and those results are filtered out.

There where some multipart features created in the manual process and if the parts are not in the correct order (and in some case there aren't) it will not find the correct start and end points of the line and because of that not find a from and/or to feature within the defined search tolerance.

Also the fact that there are a substantial number of streets without a FullName makes things a bit more complicated, When a candidate is found that has no name it is filtered out of the results and may result in the intersecting street not to have a from or to street.

There are however, a number of case that should have a result, but do not have any. This would require some serious debugging to understand what is happening in those cases in order to understand what needs to be changed in the code to get results there without the possibility to loose results elsewhere. 

After the editing the data manually I ended up with 289 streets and and 212 have a from and/or to street as result. Some of the 77 cases may not have a from or to street, but there are cases, as mentioned before, that should have a result. 

For now let's have a look at the code and I will attache the resulting featureclass for you to revise.

#-------------------------------------------------------------------------------
# Name:        from_to_street.py
# Purpose:
#
# Author:      xbakker
#
# Created:     21/06/2017
#-------------------------------------------------------------------------------
import arcpy

def main():
    fc = r'C:\GeoNet\Streets\GeoNet Street Sample.gdb\LebStreetSample_MP2SP'
    fld_fullname = 'Fullname'
    fld_from = 'From_Street'
    fld_to = 'To_Street'
    tolerance = 5  # 5 feet distance tolerance for finding intersections
    min_angle = 45  # to detect angle intersections
    sr = arcpy.Describe(fc).spatialReference

    flds = ('OID@', 'SHAPE@', fld_fullname)
    dct = {r[0]: [r[1], r[2]] for r in arcpy.da.SearchCursor(fc, flds)}
    print "dict filled..."

    cnt = 0
    dct_res = {}  # fill with results
    for oid, lst in dct.items():
        polyline = lst[0]
        full_name = lst[1]
        if full_name != ' ':
            cnt += 1
            if cnt % 50 == 0:
                print "Processing:", cnt
            from_pnt = polyline.firstPoint
            to_pnt = polyline.lastPoint
            # get a list of candidates for from and to points
            from_oids = SelectFromToCandidates(oid, full_name, from_pnt, dct, tolerance, sr)
            to_oids = SelectFromToCandidates(oid, full_name, to_pnt, dct, tolerance, sr)
            print
            print "  - before filter: ",  full_name, oid, from_oids, to_oids
            if len(from_oids) > 0:
                from_oids = FilterListOnAngles(dct, oid, from_oids, min_angle)
            if len(to_oids) > 0:
                to_oids = FilterListOnAngles(dct, oid, to_oids, min_angle)

            print "  - after filter: ", full_name, oid, from_oids, to_oids
            if len(from_oids) == 0:
                from_streets = ''
            else:
                from_streets = ', '.join([dct[o][1] for o in from_oids])
            if len(to_oids) == 0:
                to_streets = ''
            else:
                to_streets = ', '.join([dct[o][1] for o in to_oids])

            if from_streets in [' ' , ' ,  ']:
                from_streets = ''
            if to_streets in [' ' , ' ,  ']:
                to_streets = ''

            print oid, " - From:", from_streets, from_oids
            print oid, " - To  :", to_streets, to_oids
            dct_res[oid] = [oid, from_streets, to_streets]

    print "update cursor..."
    cnt = 0
    updates = 0
    flds = ('OID@', fld_from, fld_to)
    with arcpy.da.UpdateCursor(fc, flds) as curs:
        for row in curs:
            cnt += 1
            oid = row[0]
            if oid in dct_res:
                result = dct_res[oid]
                curs.updateRow(tuple(result))
                updates += 1

            if cnt % 50 == 0:
                print"processing row:", cnt, " - updates:", updates


def FilterListOnAngles(dct, oid, candidates, angle_tolerance):
    polyline = dct[oid][0]
    oids_ok = []
    for candidate in candidates:
        polyline_cand = dct[candidate][0]
        if ValidateAngle(polyline, polyline_cand, angle_tolerance, oid):
            oids_ok.append(candidate)
    return oids_ok


def ValidateAngle(polyline1, polyline2, angle_tolerance, oid):
    try:
        # get start segment and end segment
        start1 = polyline1.segmentAlongLine(0, 1, False)
        end1 = polyline1.segmentAlongLine(polyline1.length-1, polyline1.length, False)
        # sr = polyline1.spatialReference

        # determine distance start segment and end segment to intersecting line
        d_start = start1.distanceTo(polyline2)
        d_end = end1.distanceTo(polyline2)
        if d_start < d_end:
            # start is closer to intersecting line
            pntg = polyline2.queryPointAndDistance(polyline1.firstPoint, False)[0]
            seg1 = start1
        else:
            # end is closer to intersecting line
            pntg = polyline2.queryPointAndDistance(polyline1.lastPoint, False)[0]
            seg1 = end1

        # get position of projected point on intersecting line
        distance_along = polyline2.measureOnLine (pntg, False)
        # and extract segment
        seg2 = polyline2.segmentAlongLine(max([distance_along -1, 0]),
                            min([distance_along +1, polyline2.length]), False)

        print "seg1.length", seg1.length
        print "seg2.length", seg2.length

        # compare
        compare_angles = [seg1, seg2]
        angle_segs = GetAngleSegments(compare_angles)
        print "angle_segs 1", angle_segs

        while angle_segs < 0:
            angle_segs += 180

        while angle_segs > 360:
            angle_segs -= 360

        print "angle_segs 2", angle_segs
        if IsBetween(angle_segs, angle_tolerance, 180 - angle_tolerance) or IsBetween(angle_segs, angle_tolerance + 180, 360 - angle_tolerance):
            return True
        else:
            return False
    except Exception as e:
        print "ValidateAngle ERROR @ OID:", oid
        return False

def IsBetween(val, val_min, val_max):
    if val >= val_min and val <= val_max:
        return True
    else:
        return False

def GetAngleSegments(segments):
    try:
        angles = []
        for seg in segments:
            pntg1 = arcpy.PointGeometry(seg.firstPoint, seg.spatialReference)
            pntg2 = arcpy.PointGeometry(seg.lastPoint, seg.spatialReference)
            angle = GetAngle(pntg1, pntg2)
            angles.append(angle)

        #print "angles", angles
        angle_between = angles[0] - angles[1]
        return angle_between
    except:
        return None

def GetAngle(pntg1, pntg2):
    '''determine angle of line based on start and end points geometries'''
    return pntg1.angleAndDistanceTo(pntg2, method='PLANAR')[0]

def SelectFromToCandidates(cur_oid, cur_full_name, pnt, dct, tolerance, sr):
    oids = []
    pntg = arcpy.PointGeometry(pnt, sr)
    for oid, lst in dct.items():
        full_name = lst[1]
        if oid != cur_oid:
            # don't include the current polyline in results
            if full_name != cur_full_name:
                # don't include street with the same name
                polyline = lst[0]
                if polyline.distanceTo(pntg) <= tolerance:
                    oids.append(oid)
    return list(set(oids))

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

So let's see what happens in the script (and to be honest, I'm not very happy with it as it is, since there is a lot of room for enhancement).

  • Lines 12 to 17 set the location of the featureclass that will be updated, the names of the relevant fields used in the process and there are two configurations regarding the intersecting streets. First we have a search distance "tolerance" that will look in this case within 5 map units of the start and end point to see if there are any intersecting lines. The second is the minimum angle "min_angle". This is used to see if the lines found near the start and end point intersect (are heading towards perpendicular lines or are extensions in the same direction=. In this case I require a minimum angle of 45° for a line to be considered an intersecting line. Play with these settings as you wish.
  • Lines 20 and 21 will create a master dictionary using the OID as key value and a list with the polyline and the full name as value. Since we are working with a small set of streets and the procedure will loop through the entire set of streets to find near streets, this is faster. In case you want to run this on a large dataset, it might be wise to change this and take advantage of the SelectByLocation capabilities to make things run faster.
  • At line 26 a loop through the dictionary items is started and for each feature the following steps are performed:
    • get from and to point from the polyline
    •  use the function SelectFromToCandidates to generate a list of candidate features that are within the tolerance of the from point and another list for the to point. In this process it will filter out the current polyline and streets with the same name as the current polyline
    • Next the function FilterListOnAngles will be used to trigger the function ValidateAngle that will do the following:
      • extract a segment at start and end
      • determine which is nearest to the candidate (start of end segment)
      • project the start or end point from the current line onto the candidate intersecting line
      • extract a segment from the intersecting polyline around the point snapped to that line
      • get the angles from both segments and determine the angle between
      • if the angle is within the defined tolerance it remains in the list of candidates, else it is removed
    • Finally the name(s) of the intersecting streets are determined and stored in a directionary.
    • The dictionary will be used in an update cursor to write the results back to the featureclass.

This is basically what goes on in the script.

Attached you will find the result.

View solution in original post

28 Replies
ChrisDonohue__GISP
MVP Frequent Contributor

To help solve this, what sort of information is in "On Street"?  Is the "To" and "From" address range there? 

Addressing 

Chris Donohue, GISP

0 Kudos
DanPatterson_Retired
MVP Esteemed Contributor

a blog from a few years ago

https://blogs.esri.com/esri/arcgis/2013/02/11/streetintersectionlist/

I am sure there are other references but the mechanics will undoubtedly depend on your existing data

0 Kudos
KevinCross2
New Contributor III

Below is a sample of what the data looks like.  I inherited this and am trying to fix it!  Currently nothing has been filled out in the "From Street or the To Street".  

0 Kudos
ChrisDonohue__GISP
MVP Frequent Contributor

All our data has is "On Street"

There does not appear to be an "On Street" field in the example.  Can you elaborate on what is in the "On Street" field currently?  Does it contain the street address ranges?  If so, it may be possible to write a script to extract out the "To Street" and "From Street" values from the "On Street" field.  However, we will need to be able to see the data and how it is arranged in the "On Street" field to evaluate this and come up with the correct code.

If the "On Street" field does not contain the range and there are no other fields in the data that contain the street address range information, then the challenge will instead be to derive the ranges.  This may be possible if you have other data.  For example, if you have address range data in an Intersection points layer, it can be spatial-joined to the lines and then calculated across.  Again, it will depend on what data you have and how it is arranged.

Chris Donohue, GISP

0 Kudos
KevinCross2
New Contributor III

My bad, the "On Street" Field would be the "Full Name" Field in the example above.

0 Kudos
XanderBakker
Esri Esteemed Contributor

kevin.cross , if you are willing to share a sample of your data, I'm sure we can have a look at it to see what can be done.

0 Kudos
KevinCross2
New Contributor III

Xander,  What would be the best way to get you the data?

0 Kudos
XanderBakker
Esri Esteemed Contributor

Hi Kevin Cross , you could select a subset of your data, export it to a new File Geodatabase, and ZIP that fgdb. Make sure it is a representative subset of the entire dataset.

To attach the ZIP you can either edit the original post or to this comment and use the advanced editor (link upper right corner when replying to the thread) which will reveal the Attach option in the lower left corner of the advanced editor.

KevinCross2
New Contributor III

Xander,  hopefully I attached this correctly!   I appreciate your help on this!

0 Kudos