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