Accepted Solutions
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.
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".
