[Matplotlib-users] Get distances between points in basemap?

[Scott H. <st...@co...>]

I'm trying to efficiently get the distances of all points on a map to a 
specified point. If the map is in projected coordinates, what is the 
best way of going about this? Is there is a 'standard' way to get the 
distance between points through internal basemap functions? After some 
scrounging around what I have below works for individual points, but is 
there a way to avoid the big for loop?

Any advice would be appreciated.
Thanks.


# --------------------------------------
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import pyproj
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)

LL = (-69, -26)
UR = (-66, -21)
map = Basemap(projection='merc',
                llcrnrlat=LL[1],
                urcrnrlat=UR[1],
                llcrnrlon=LL[0],
                urcrnrlon=UR[0],
                resolution='i',
                suppress_ticks=False,
                ax=ax)

lons, lats, xs, ys = map.makegrid(200, 200, returnxy=True)
lon1, lat1 = (-67.28, -23.78)

gc = pyproj.Geod(a=map.rmajor, b=map.rminor)
#azis12, azis21, distances = gc.inv(lon1,lat1,lons,lats) #doesn't work 
with vectors or matrices?
distances = np.zeros(lons.size)
for i, (lon, lat) in enumerate(zip(lons.flatten(),lats.flatten())):
     azi12, azi21, distances[i] = gc.inv(lon1,lat1,lon,lat)
distances = distances.reshape(200,200) / 1000.0 #in km

# Plot perimeters of equal distance
levels = [50] #[50,100,150]
map.drawcountries()
x1,y1 = map(lon1,lat1)
map.plot(x1,y1,'m*')
cs = map.contour(xs, ys, distances, levels)

-- 
---------------
Scott T. Henderson
http://www.geo.cornell.edu/eas/gstudent/sth54/contact.html