What algorithm would be most efficient when trying to find the nearest city given a set of coordinates?

I have a dataset which contains the longitude and latitude of the 1000 largest US cities. I'm designing an API which returns the user's nearest city, given an input of the user's longitude/latitude.

What is the most efficient algorithm I can use to calculate the nearest city? I know that I can use the haversine formula to calculate the distance between the user's coordinate and each cities, but it seems inefficient to have to do this for all 1000 cities. I've previously used a k-d tree to solve nearest neighbour problems on a plane - is there a similar solution that can be used in the context of a globe?

Edit: keeping this simple - distance I'm looking for is as the crow flies. Not taking roads or routes into account at this stage.

First, https://gisgeography.com/python-libraries-gis-mapping/ has a bunch of Python libraries for dealing with geography. Or you can push the problem to a database, for example with https://www.gaia-gis.it/fossil/libspatialite/index, and solve it there. I'd recommend looking at those ideas before wrapping your own.

That said, in the past what I've done for US cities is to stick everything into a lat/long 2-D structure, take the city, construct a lat/long bounding box (using the fact that lat/long can be turned into approximate distances), use that to return a handful of candidate closest cities, and then do the hard calculation from the candidate cities. I've personally seen a database successfully perform this calculation several hundred thousands times/hour while handling a lot of other traffic. And this was 20 years ago - I'd expect more now.

For a worldwide dataset I'd put the cities into a handful of such structures, with different choices of poles for lat/long. For a given city I'd pick the axis which puts the city farthest from either pole, then proceed as before.

You can split the map into squares that do not overlap and they cover the whole US map (i.e., you will have a grid). You will number the squares using the coordinates of their upper left corner (i.e., each one will have a unique ID) and you will do a preprocessing where each city will be assigned with the ID of the square where it belongs. You will find the square where the user lies into and then you will check only the cities that lie into this square and the ones that are one step from this (total: 9 squares). If these are empty of cities, you will check the ones that are two steps of it etc. In this way, on average you will check much less cities to find the closest

This answer is very similar to that of ckc.

First, spilt the 1000 cities in 2 groups : a big one located located between Canada and Mexico and the few others cities located outside this rectangle (i.e Alaska, Hawai, ...).

When processing coordinates, check if they belong to the small group : in this case, no optimisation needed.

To optimize the other case, you may divide the map in rectangles (example 5°lat x 7° lon) and associate to each rectangle the list of cities belonging to each rectangle.

To find the nearest city, consider the rectangle R containing the point.
Compute the distance to the cities of the rectangle.
Process the 8 rectangles adjacent to R by computing the distance of the point to each rectangle : you may then eliminate the adjacent rectangles whose distance is greater than the best distance already found.

Iterate the process to a next level, i.e. the next crown (rectangles located on the outside of the area composed of 5x5 rectangles whose center is R).

There is some confusion in the discussion, here is an example using Haversine, with sklearn.

import pandas as pd
import numpy as np

from sklearn.neighbors import BallTree

cities = pd.read_json("https://gist.githubusercontent.com/Miserlou/c5cd8364bf9b2420bb29/raw/2bf258763cdddd704f8ffd3ea9a3e81d25e2c6f6/cities.json")

cities_gps = cities[['latitude','longitude']].values
cities_radians = np.radians(cities_gps)

The above code loads 1000 cities I found on some git. Pandas is just being used to parse the JSON, it is not needed to calculate.

You can build a tree as follows (leafsize can be tweaked a little, has minor effect)

tree = BallTree(cities_radians, leaf_size=15, metric='haversine')

Now say we have 10000 lat/longs, you can query as follows:

random_geo = np.random.normal(loc=(30,-80), scale=(8,8), size=(10000,2))
random_geo_radians = np.radians(random_geo)

If you just query 1, make sure not to re-build the tree for each query. That probably takes more time than query for just 1 lat/long pair. SO keep the tree alive and fire queries at it.

distances, idx = tree.query( random_geo_radians, k=1)

The distances contains the unit circle distance, you need to convert to miles/km by multiplying with the earth radius. The mean radius of the USA would be best.

You could get the names with

cities.city[ idx[:,0] ]