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