[GIS] Shapely python spatial join – point with line

Question

I have two datasets, one with address points (lat,long) coordinates. And another one with line string geometry (each linestring represents a road). I am trying to find the road name for each address point. I'm trying to do spatial join or find the road line with shortest distance to the address point using shapely line.distance(point). How can I achieve this in python? One thing I can think of is to find the distance between the address point and every line in road dataset. And get the one with minimal distance. Is there an alternate approach to this?

Both the dataset has lat long points in wgs84 projection system.

Line coordinate example :

LINESTRING (168.0473, -44.40220, 168.0485, -44.40117, 168.0486, -44.40085, 168.0487, -44.40067)

Point coordinate example:

(168.04665153,-44.3990252)

Answer 1

If you have datasets with (long,lat) coordinates, you need to project them prior to distance calculations. For instance, next code opens a point and a line layers previous to these calculations. Results (all distances and minimal distances and its points) are printed by console.

import fiona
from shapely.geometry import shape, Point, LineString
import pyproj

srcProj = pyproj.Proj(init='EPSG:4326')
dstProj = pyproj.Proj(init='EPSG:32612')

path1 = '/home/zeito/pyqgis_data/points_test_4326.shp'
path2 = '/home/zeito/pyqgis_data/new_lines_4326.shp'

points = fiona.open(path1)
line = fiona.open(path2)

points = [ (shape(feat["geometry"]).xy[0][0], shape(feat["geometry"]).xy[1][0]) 
           for feat in points ]

lines = [ zip(shape(feat["geometry"]).coords.xy[0], shape(feat["geometry"]).coords.xy[1]) 
          for feat in line ]

proj_lines = [ [] for i in range(len(lines)) ]

for i, item in enumerate(lines):
    for element in item:
        x = element[0]
        y = element[1]
        x, y = pyproj.transform(srcProj, dstProj, x, y)
        proj_lines[i].append((x, y))

proj_points = []

for point in points:
    x = point[0]
    y = point[1]
    x, y = pyproj.transform(srcProj, dstProj, x, y)    
    proj_points.append(Point(x,y))

distances = [ [] for i in range(len(lines)) ]

for i, line in enumerate(proj_lines):
    for point in proj_points:
        distances[i].append(LineString(line).distance(point))

print distances

for i, list in enumerate(distances):
    print "minimal distance: {:.2f} to this point: {}".format(min(list), points[i])

After running the code, it was gotten:

[[240.76129848041424, 82.17367496992877, 534.8119728316819], [180.2963692590486, 601.7275241365959, 1155.1004277619784]]
minimal distance: 82.17 to this point: (-112.17112916813936, 40.17311785635686)
minimal distance: 180.30 to this point: (-112.1560813471297, 40.170581085323384)

Layers look like at next image. QuickWKT plugin of QGIS helped to corroborate that these determined points are corresponding to minimal distance.