Created
October 25, 2011 14:06
-
-
Save pmarshwx/1312843 to your computer and use it in GitHub Desktop.
Find the shortest geometric distance between four observations and 1148 forecast locations
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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] |
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.
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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.