[Python] Calculating the closest intersection based upon latitude and longitude

I have a geojson file containing the latitude and longitude coordinates of all of the streets and avenues in New York City - they're all formatted as either LineString and MultiLineString as follows:

{
"type": "FeatureCollection",
"features": [
{ "type": "Feature", "properties": { "STATEFP": "36", "COUNTYFP": "005", "LINEARID": "110391528508", "FULLNAME": "Longwood Ave", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ -73.900023, 40.818595 ], [ -73.899818, 40.81847 ] ], [ [ -73.899818, 40.81847 ], [ -73.899285, 40.818147 ], [ -73.898521, 40.817686 ], [ -73.897786, 40.817246 ], [ -73.897045, 40.816802 ], [ -73.896213, 40.816319 ], [ -73.895801, 40.816066 ], [ -73.895428, 40.815845 ], [ -73.895158, 40.815683 ], [ -73.894986, 40.815581 ], [ -73.894656, 40.815383 ], [ -73.894014, 40.814998 ], [ -73.893071, 40.814438 ], [ -73.891831, 40.813701 ], [ -73.891434, 40.813466 ], [ -73.890975, 40.813213 ] ] ] } },
{ "type": "Feature", "properties": { "STATEFP": "36", "COUNTYFP": "005", "LINEARID": "110391524085", "FULLNAME": "E 149th St", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.917625, 40.816055 ], [ -73.916771, 40.815699 ], [ -73.914954, 40.814937 ], [ -73.912954, 40.814269 ], [ -73.911812, 40.813883 ], [ -73.910948, 40.813621 ], [ -73.910019, 40.813432 ], [ -73.909075, 40.813233 ], [ -73.908159, 40.813044 ], [ -73.907245, 40.812855 ], [ -73.90633, 40.812667 ], [ -73.905413, 40.812476 ], [ -73.904466, 40.812282 ], [ -73.90418, 40.812125 ], [ -73.903417, 40.811507 ] ] } },
{ "type": "Feature", "properties": { "STATEFP": "36", "COUNTYFP": "005", "LINEARID": "110391528025", "FULLNAME": "Timpson Pl", "RTTYP": "M", "MTFCC": "S1400" }, "geometry": { "type": "LineString", "coordinates": [ [ -73.906468, 40.809563 ], [ -73.906252, 40.809719 ], [ -73.90507, 40.810487 ], [ -73.903417, 40.811507 ], [ -73.901093, 40.81218 ], [ -73.899167, 40.812665 ] ] } },


Given latitude and longitude, what would be the best way to find out the closest intersection of two streets (54th street & 8th avenue, 23rd street and Park Avenue, etc.)

I've tried the following using the Haversine formula but am not sure how to conceptualize this at a closer level.

import json
from geopy.distance import geodesic
from shapely.geometry import Point, shape


def find_closest_point(geojson_file, target_lat, target_lon):
with open(geojson_file, "r") as f:
geojson_data = json.load(f)

target_point = (target_lat, target_lon)
min_distance = float("inf")
closest_feature = None

for feature in geojson_data["features"]:
if (
feature["geometry"]["type"] == "MultiLineString"
or feature["geometry"]["type"] == "LineString"
):
coords = feature["geometry"]["coordinates"]
feature_point = (coords[1], coords[0])
distance = geodesic(target_point, feature_point).meters

if distance < min_distance:
min_distance = distance
closest_feature = feature

return closest_feature


geojson_file = "nyc.geojson"
target_latitude = 40.748687497557995 # Times Square
target_longitude = -73.98545265197755

closest_point = find_closest_point(geojson_file, target_latitude, target_longitude)

if closest_point:
print("Closest Feature:", closest_point)
else:
print("No points found in GeoJSON")




Continue reading...