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

[Jeff W. <js...@fa...>]

On 8/20/12 10:26 PM, Scott Henderson wrote:
> 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

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