Find the nearest point features in a shapefile to a given coordinate

nearest_features.py

import numpy as np
from scipy.spatial import KDTree
from fiona import collection

def nearest_features(datasource, x, y, n=1):
    """
    datasource: The ogr datasource path
    x, y: coordinates to query for nearest
    n: number of features (default = 1 or the single nearest feature)

    returns a (zipable) tuple of lists containing
     ([features, ...], [distances from query point, ...])

    example: looking for the nearest 3 ports to -122 42
    >>> features, distances = nearest_features('world_ports.shp', -122.0, 42.0, n=1)
    >>> print [x['properties']['MAIN_PORT_'] for x in features] # name of closest port
    [u'CRESCENT CITY']
    >>> features, distances = nearest_features('world_ports.shp', -122.0, 42.0, n=5)
    >>> print [x['properties']['MAIN_PORT_'] for x in features] # names of 5 closest ports
    [u'CRESCENT CITY', u'SAMOA', u'EUREKA', u'FIELDS LANDING', u'COOS BAY']
    """
    with collection(datasource, "r") as source:
        features = list(source)

    pts = np.asarray([feat['geometry']['coordinates'] for feat in features])
    tree = KDTree(pts)

    querypoint = np.asarray([[x, y]])
    result = tree.query(querypoint, n)

    if n == 1:
        nearest_features = [features[result[1][0]],]
        distances = list(result[0])
    else:
        nearest_features = [features[x] for x in result[1][0]]
        distances = list(result[0][0])

    return nearest_features, distances

if __name__ == "__main__":
    import doctest
    doctest.testmod()