Calculating the closest intersection based upon latitude and longitude

Question

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")

enter image description here

Answer 1

Iterate the points on the line calculating the distance between then. This should be quite straighforward using UTM coordinates, but you probably can find some ready-to-go code in a lib.

Once you have a distance that's below a certain threshold, say 5 meters, you have your intersection.

Remember GPS in real world usage often max out the precision within 10m. Using multiple constelations 2-3 meters.

Answer 2

Not sure if you want to find the nearest road or the nearest intersection.

This is how you can find the nearest road to an input point using SpatialIndex.nearest

import geopandas as gpd
from shapely.geometry import Point
road_df = gpd.read_file(r"C:\Users\bera\Downloads\nyc-streets.geojson")
#Reproject to a projected coordinate system in meters
road_df = road_df.to_crs(road_df.estimate_utm_crs()) 
#Create a point Geoseries in the same crs
target_latitude = 40.748687497557995
target_longitude = -73.98545265197755
pnt = gpd.GeoSeries(data=[Point(target_longitude, target_latitude)], 
 crs=4326).to_crs(road_df.crs)
#Create a spatial index on the roads
si = road_df.sindex
#Find nearest road 
result = si.nearest(geometry=pnt, return_distance=True)
# (array([[ 0],
# [8045]]),
# array([11.80314789]))
#The nearest road has index 8045 and is 11.8 m away
road_df.iloc[8045]["FULLNAME"]
#'W 34th St'