Road points with Shapely

import sys; sys.version.split()[0]
0.2s
Python
Python 3
'3.7.5'

I uploaded our road and point shapefiles:

29_11_2019 - Master Shapefile.rar
Ward 106.rar

I use command line tool 'unrar' to open your RAR files, putting the shapefile in local directory. I use NextJournal's Ctrl+E to put the custom file names:

unrar -y e NJ__REF_
0.8s
Bash in Python
Python 3
unrar -y e NJ__REF_
2.2s
Bash in Python
Python 3
ls *.shp
0.5s
Bash in Python
Python 3

I install dependencies using pip (this is still in the command line)

pip install shapely pyshp
4.5s
Bash in Python
Python 3

In Python, let's load all of our data:

import shapefile
import shapely
by_id = {}
sf = shapefile.Reader("W106 Points.shp")
for f in sf.shapeRecords():
    lnglat = f.shape.points[0]
    #print(f.record)
    id = f.record.id
    by_id[id] = (lnglat[0], lnglat[1])
from shapely.geometry import LineString
sf2 = shapefile.Reader("Road OSM.shp")
roads = []
road_data = []
for road in sf2.shapeRecords():
    ln = []
    for pt in road.shape.points:
        ln.append((pt[0], pt[1]))
    roads.append(LineString(ln))
    road_data.append(road.record)
3.4s
Python
Python 3

This is the part where we enter the start and end points, and draw a box connecting them,

segment_start = 63
segment_end = 65
radius = 0.00001
my_segment_box = LineString([by_id[segment_start], by_id[segment_end]]).buffer(radius)
0.2s
Python
Python 3

Now let's look at all of the roads which intersect with our box, and return the one that has the most overlap:

rdindex = 0
best_inter = 0
best_road = None
for road in roads:
    road_box = road.buffer(radius)
    inter = road_box.intersection(my_segment_box).area
    if inter > best_inter:
        best_inter = inter
        best_road = road_data[rdindex]
    rdindex += 1
print(best_road)
print(best_road[4])
5.6s
Python
Python 3
Runtimes (1)