[AstroPy] all-sky map

Matthew Turk matthewturk at gmail.com
Tue Apr 5 14:56:35 EDT 2011

Hi Marsall,

On Tue, Apr 5, 2011 at 2:50 PM, Marshall Perrin <mperrin at stsci.edu> wrote:
> On Apr 5, 2011, at 12:32 PM, Jonathan Slavin wrote:
>> I'm looking for a way to plot an all-sky map of modeled data using a
>> Hammer-Aitoff projection.  The way I've done this in IDL is to create a
>> uniform x-y grid, translate that l, b using the proper conversion for an
>> Aitoff projection and generate the data on that grid.  I then display
>> the image and overlay an Aitoff grid (which I also generate).  So the
>> image is rectangular and extends beyond the plot edges.  That is fine,
>> but then all the labeling, etc. has to be done by hand.  Is there an
>> easier way?  Any help would be appreciated.
> Check out the matplotlib mpl_toolkits.basemap module.   The Basemap class implements a user-selectable map projection, and yields a callable object which handles the translation between projection coordinate systems and x,y positions for plotting.   Here is some code I recently wrote to do a similar task, plotting positions of objects on an all-sky Mollweide projection:

In theory one should be able to get this going with just the standard
Matplotlib projections, without basemap, which has a huge dependency
set.  There seems to be a bug if you don't set the figure to have
aspect ratio 2:1, but this code works for me, where img is a square

import matplotlib.figure
import matplotlib.backends.backend_agg

fig = matplotlib.figure.Figure((10, 5))
ax = fig.add_subplot(1,1,1,projection='mollweide')
image = ax.imshow(img, extent=(-pi,pi,-pi/2,pi/2), clip_on=False, aspect=0.5)
cb = fig.colorbar(image, orientation='horizontal')

canvas = matplotlib.backends.backend_agg.FigureCanvasAgg(fig)

(a bit more discussion:

If anybody knows why the figure size has to be 2:1 to avoid a blank
image, I'd really appreciate a tip.


>        # define base map class.
>        map = Basemap(projection='moll', lat_0 = 0, lon_0 = 0,
>                              resolution = None)  # do *NOT* draw Earth continents at any resolution!
>        map.drawmapboundary()
>        p.title("Equatorial coordinates J2000")
>        # draw and label ra/dec grid lines every 30 degrees.
>        degtoralabel = lambda deg : "%+d$^h$" % int(deg/15)
>        degtodeclabel = lambda deg : "%+d$^\circ$" % deg
>        map.drawparallels(n.arange(-90, 90, 30), fmt=degtodeclabel, labels=[1,0,0,0])
>        map.drawmeridians(n.arange(0, 360, 30) )  # label these manually since I don't like the default label positions:
>                                                                                        # this also demonstrates how to overplot text on map coordinates...
>        for h in [0,6,12,18]:
>            x,y = map(h*15,0)
>            p.text(x,y, degtoralabel(h*15))
>        # draw data points
>        px, py = map(data.radeg, data.dedeg)
>        map.plot(px, py, "o", color="red")
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy

More information about the AstroPy mailing list