Find the shortest geometric distance between four observations and 1148 forecast locations

find_distance.py

import numpy as np



def great_circle(lon1, lat1, lon2, lat2):
    delta = lon2 - lon1
    a = np.radians(lat1)
    b = np.radians(lat2)
    c = np.radians(delta)
    x = np.sin(a) * np.sin(b) + np.cos(a) * np.cos(b) * np.cos(c)
    dist = np.arccos(x)         # in radians
    dist = np.degrees(dist)     # in degrees
    dist = dist * 60            # degrees to nautical miles
    dist = dist * 1.15077945    # nautical miles to miles
    return dist

def conversion(degree, minute, second):
    var = []
    for x,y,z in zip(degree, minute, second):
        var.append(x + y/60. + z/3600.)
    return var


data = np.genfromtxt('soundingsites.txt', dtype=None, names=True)
number = data["number"]
names = data["name"]
flons = data["lon"]
flats = data["lat"]
felevations = data["elevation"]


latdeg = [36, 36, 36, 36]
latmin = [34, 29, 46, 36]
latsec = [43.25, 28.39, 3.18, 19.01]

londeg = [-97, -97, -97, -97]
lonmin = [-21, -35, -32, -29]
lonsec = [-49.64, -37.86, -50.93, -8.43]

lats = conversion(latdeg, latmin, latsec)
lons = conversion(londeg, lonmin, lonsec)

distance1 = []
distance2 = []
distance3 = []
distance4 = []

for flon,flat in zip(flons, flats):
    distance1.append(great_circle(lons[0], lats[0], flon, flat))
    distance2.append(great_circle(lons[1], lats[1], flon, flat))
    distance3.append(great_circle(lons[2], lats[2], flon, flat))
    distance4.append(great_circle(lons[3], lats[3], flon, flat))

distance1 = np.array(distance1)
distance2 = np.array(distance2)
distance3 = np.array(distance3)
distance4 = np.array(distance4)

ind1 = np.where(distance1 == distance1.min())
ind2 = np.where(distance2 == distance2.min())
ind3 = np.where(distance3 == distance3.min())
ind4 = np.where(distance4 == distance4.min())

print (lons[0], lats[0]), distance1[ind1], names[ind1], flons[ind1], flats[ind1], felevations[ind1]
print (lons[1], lats[1]), distance2[ind2], names[ind2], flons[ind2], flats[ind2], felevations[ind2]
print (lons[2], lats[2]), distance3[ind3], names[ind3], flons[ind3], flats[ind3], felevations[ind3]
print (lons[3], lats[3]), distance1[ind4], names[ind4], flons[ind4], flats[ind4], felevations[ind4]

I take it that this runs like a slug? Try converting everything to cartesian first (any reference point will do). Then use scipy's kdtree for highly efficient point comparisons.

Yeah, I'm sure it does run slow. However, it's a code sample for an informal Python course I'm teaching some OU seniors. I gave them an assignment and this is the simple solution. You can see the assignment here: http://www.patricktmarsh.com/teaching/python-for-meteorology/. We're only comparing against 1148 points, so it still takes less than a second to complete.