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)
>
------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Matplotlib-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
