Can I optimize two for loops that look for the closest zip code based on lat/lon?

Question

I am new to python and I had the task to find the US zip code based on latitude and longitude. After messing with arcgis I realized that this was giving me empty values for certain locations. I ended up coding something that accomplishes my task by taking a dataset containing all US codes and using Euclidean distance to determine the closest zip code based on their lat/lon. However, this takes approximately 1.3 seconds on average to compute which for my nearly million records will take a while as a need a zip code for each entry. I looked that vectorization is a thing on python to speed up tasks. But, I cannot find a way to apply it to my code. Here is my code and any feedback would be appreciated:

for j in range(len(myFile)):
    p1=0
    p1=0
    point1 = np.array((myFile["Latitude"][j], myFile["Longitude"][j]))  # This is the reference point
    i = 0
    resultZip = str(usZips["Zip"][0])
    dist = np.linalg.norm(point1 - np.array((float(usZips["Latitude"][0]), float(usZips["Longitude"][0]))))
    for i in range(0, len(usZips)):
        lat = float(usZips["Latitude"][i])
        lon = float(usZips["Longitude"][i])
        point2 = np.array((lat, lon))  # This will serve as the comparison from the dataset
        temp = np.linalg.norm(point1 - point2)
        if (temp <= dist):  # IF the temp euclidean distance is lower than the alread set it will:
            dist = temp  # set the new distance to temp and...
            resultZip = str(usZips["Zip"][i])  # will save the zip that has the same index as the new temp
            # p1=float(myFile["Latitude"][58435])
            # p2=float(myFile["Longitude"][58435])
        i += 1

I am aware Google also has a reverse geocoder API but it has a request limit per day. The file called myFile is a csv file with the attributes userId, latitude, longitude, timestamp with about a million entries. The file usZips is public dataset with information about the city, lat, lon, zip and timezone with about 43k records of zips across the US.

Answer 1

for j in range(len(some_list)) is always a red flag. If you need the index, you can use for j,element in enumerate(some_list). If you don't need the index (like here), you should just use for element in some_list. (Speaking of which, csv files don't have lengths that you can iterate over. Numpy arrays or Pandas DataFrames do.)

You have i=0, for i in range and i+=1. You only want the middle of those (and i+=1 will cause a bug if it's even executed). But you don't even want the middle, because of the previous paragraph. Also, you need to guard against returning that your closest zip code is yourself.

You flip around between floats and strings. I find my life is much easier the sooner I convert my data to the appropriate format. This means that you should convert your latitude and longitude to floats as soon as possible (probably when you load the csv file). I don't see a problem leaving the zip code as an int until later. Just be aware that you will need to format it at some point: some zip codes start with 0.

Instead of pre-seeding dist with the first zip code, it might be clearer to start dist off as float('inf'). But:

The whole point of the inner for loop is to find the closest point. We can use min to do that, as long as we tell Python what we want to minimize, and remove the inner loop entirely (from the code, not the runtime).

What is p1? What do you do with resultZip?

Your distance calculation is assuming a flat earth. That breaks down more with higher latitudes (for example, all longitudes at latitude 89 are fairly close together). Also, calculating the square of the distance is usually faster, since most involve taking a square root at the end. (Also, letting the distance function take rows instead of points will allow us to simplify some other parts.) A good approximation for a spherical earth (instead of ellipsoidal) would be

def dist_sq(row1,row2):
    if row1==row2:
        return float('inf')
    dx = (row1['longitude']-row2['longitude'])*math.cos((row1['latitude']+row2['latitude'])/2)
    dy = row1['latitude']-row2['latitude']
    return dx*dx+dy*dy

That leaves us with your logic boiled down to:

import functools
for row in myFile:
    resultZip = min(myFile,key=functools.partial(dist_sq,row))['Zip']

Answer 2

The vectorized version of your code is

def distances(point1,i):
    lat = float(usZips["Latitude"][i])
    lon = float(usZips["Longitude"][i])
    point2 = np.array((lat, lon))
    temp = np.linalg.norm(point1 - point2)
    
[min(range(0, len(usZips)), key = lambda i: distance(
        np.array((myFile["Latitude"][j], myFile["Longitude"][j])) , 
        i))
 for j in range(len(myFile)) ]

However, that's not very Pythonic, as it's iterating over indices rather than elements. More importantly, it's still rather inefficient.

Consider the following algorithm: you take every combination of integer lat/long in the US. That's somewhere around 1000. For each one, you create a list of nearby zip codes. Then, when you're given a lat/long pair, you find the four integer lattice points around it, and take only the zip codes that you previously found are "close" to those. (You should also separate out Alaska and Hawaii zip codes and handle them separately). This takes more initial computation to calculate closest zip code for each lattice point, but you'll probably make it up again when going through all the locations.

You can expand this algorithm to more lattice points, and you can also iterate. For instance, start with 40 lattice points. Then when you go to 1000 lattice points, use the four closest of the 40 lattice points to find the zip code candidates for the 1000 lattice points. Then go to 100,000 lattice points and use the four closest of the 1000 lattice points for the candidates for the 100,000 lattice points.

Something else to keep in mind find "closest zip code" based on distance is not going to necessarily give you the zip code a location is in.

Answer 3

You need a spatial structure to speed up lookup for closest zip code. I would approach this problem in two steps:

1. Create a point quad tree and add all known zip codes to it. Latitude and longitude can be as used as coordinates of the tree node and zip code as value.
2. For every zip entry in your myFile, query the quad tree for closest zip code.

There are many implementations of quad tree available as Python libraries. You can pick one of them.

Answer 4

So a better approach may be to :

1. covert your myFile to a Numpy array [[lat0, long0], [lat1, long10], ...].

lat_long = np.asarray(myFile)

2. Sort this array on latitude

lat_long = np.sort(lat_long, axis=0)

3. find the closest western points to the target point

closest_point_to_west = lat_long[lat_long[:, 0] >= myPoint[0]]

Do the same for the other wind directions, this will give you the four closest points using Numpy's efficient slicing operations and you can do a distance measurement on only 4 of the points and not all of them.