Get distances between points in basemap?

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

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.

Scott: The pyproj Geod methods don't take numpy arrays, as you discovered. It's very accurate and robust, but slow since you have to loop over the points. Since you are assuming the earth is a perfect sphere you could use the less accurate Haversine formula:

# function to compute great circle distance between point lat1 and lon1 and arrays of points
# given by lons, lats
def get_dist(lon1,lons,lat1,lats):
     # great circle distance.
     arg = np.sin(lat1)*np.sin(lats)+np.cos(lat1)*np.cos(lats)*np.cos(lon1-lons)
     arg = np.where(np.fabs(arg) < 1., arg, 0.999999)
     return np.arccos(arg)

-Jeff

···

On 8/20/12 10:26 PM, Scott Henderson wrote:

# --------------------------------------
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)

I'd like to use the same patch to clip two images that share the same axes, and extract values from the un-clipped region of both arrays. Unfortunately this seems harder than expected. Code & questions below, Thanks!

from matplotlib.patches import Polygon
import matplotlib.pyplot as plt
from numpy import random

poly = patches.Circle((5,5), radius=3, fill=False, ec='none')

fig = plt.figure()
ax = fig.add_subplot(121)
ax.autoscale_view(0,0,0)
test = random.rand(10,10)
im = ax.imshow(test)

test1 = random.rand(10,10)
# How to prevent automatic axis scaling?
ax1 = fig.add_subplot(122, sharex=ax, sharey=ax)
im1 = ax1.imshow(test1)

ax.add_patch(poly)
im.set_clip_path(poly)

# Doesn't work b/c poly vertices auto-transformed to display coords by add_patch?
ax1.add_patch(poly)
im1.set_clip_path(poly)

# How to extract non-clipped values instead of full array?
clipped_data = im.get_array()

plt.show()