""" Python script to plot a radiant visibility map for the 2011 Draconids Author: Geert Barentsen, gba-at-arm.ac.uk Date: 2011 August 15 """ from mpl_toolkits.basemap import Basemap import matplotlib.pyplot as plt import matplotlib as mpl import numpy as np import ephem hour = "20" # Which hour to plot? mydate = "2011/10/08 %s:00:00" % hour """ Compute radiant+solar altitude across Earth using PyEphem """ radiant = ephem.FixedBody() radiant._ra = ephem.degrees('263.2') radiant._dec = ephem.degrees('55.8') sun = ephem.Sun() result = [] # list containing radiant altitudes across the world (or NaN if daylight) lon1, lon2, lonstep = -170, 190, 0.5 lat1, lat2, latstep = -40, 74, 0.5 longitudes = np.arange(lon1, lon2+lonstep, lonstep) latitudes = np.arange(lat1, lat2+latstep, latstep) for lat in latitudes: for lon in longitudes: # Create a PyEphem "Observer" at the given lon/lat observer = ephem.Observer() observer.lat = ephem.degrees('%s' % lat) observer.long = ephem.degrees('%s' % lon) observer.date = mydate # Compute radiant and solar altitude radiant.compute(observer) radiant_alt = np.degrees(radiant.alt) sun.compute(observer) sun_alt = np.degrees(sun.alt) if sun_alt > -6 or radiant_alt < 0: result.append( np.NaN ) else: result.append( radiant_alt ) """ Create the plot """ plt.Figure() plt.subplots_adjust(0.01, 0.00, 0.99, 1.0, hspace=0.0, wspace=0.0) ax = plt.subplot(111) # Show earth with Mercator projection m = Basemap(projection='merc', llcrnrlon=lon1, llcrnrlat=lat1, urcrnrlon=lon2, urcrnrlat=lat2, resolution="c", fix_aspect=False) m.drawcoastlines() # Plot the altitudes with an appropriate color map cdict = {'red' : ((0., 1., 1.), (1., 0.1, 0.)), 'green': ((0., 1., 1.), (1., 1., 1.)), 'blue' : ((0., 0., 0.), (1., 0., 0.))} my_cmap = mpl.colors.LinearSegmentedColormap('my_colormap', cdict, 1024) m.imshow(np.array(result).reshape(len(latitudes), len(longitudes)), cmap=my_cmap) # Plot the time in the upper left corner plt.text(0.02, 0.98, '%s:00 UT' % hour, horizontalalignment='left', verticalalignment="top", fontsize=32, backgroundcolor = 'white', transform = ax.transAxes) plt.show()