[AstroPy] all-sky map

Michael Droettboom mdroe at stsci.edu
Thu Apr 7 12:39:41 EDT 2011


I can confirm Erik's non-blank results with the latest git master and 
v1.0.1.  Let us know which version of matplotlib you are using -- I 
would like to determine whether it's a "bug since fixed" or some 
interaction with another library on your platform that is causing this.  
If the latter, we should get to the bottom of it and add a workaround to 
matplotlib.

Cheers,
Mike

On 04/07/2011 02:15 AM, Erik Tollerud wrote:
> Whoops, sorry, you're aboslutely right - I left out a very important
> line!  The corrected version (the key change being the addition of the
> call to "axes" with the projection keyword) is:
>
> import matplotlib
> matplotlib.use('agg') #I tried a few other backends and got the same result
>
> from matplotlib.pyplot import figure,axes,imshow,savefig
> from numpy.random import randn
> from math import pi
>
> img = randn(100,100)
>
> figure(figsize=(10,5))
> axes(projection='mollweide')
> imshow(img,extent=(-pi,pi,-pi/2,pi/2))
> savefig('filename.png')
>
>
>
> Also, when I run your script you linked to, I get the following three
> images (tested both on a Snow Leopard Mac and Ubuntu Lucid, both with
> matplotlib 1.0.1):
> http://dl.dropbox.com/u/8683962/allsky1.png
> http://dl.dropbox.com/u/8683962/allsky2.png
> http://dl.dropbox.com/u/8683962/allsky3.png
>
> And as you can see, 2 and 3 have colored imshow ovals rather than
> blank.  This makes me think this is a bug that was fixed at some point
> (if you're using a version earlier than 1.0.1), or some rendering
> problem that is platform-specific (if you're using something other
> than OS X or Linux)... Definitely very weird behavior, though!
>
>
> On Wed, Apr 6, 2011 at 4:08 PM, Matthew Turk<matthewturk at gmail.com>  wrote:
>> Hi Eric,
>>
>> On Wed, Apr 6, 2011 at 8:49 AM, Erik Tollerud<erik.tollerud at gmail.com>  wrote:
>>>> If anybody knows why the figure size has to be 2:1 to avoid a blank
>>>> image, I'd really appreciate a tip.
>>> Maybe this is a version problem? If I do the following with matplotlib 1.0.1:
>>>
>>> import matplotlib
>>> matplotlib.use('agg') #I tried a few other backends and got the same result
>>>
>>> from matplotlib.pyplot import figure,axes,imshow,savefig
>>> from numpy.random import randn
>>> from math import pi
>>>
>>> img = randn(100,100)
>>>
>>> figure(figsize=(10,5))
>>> imshow(img,extent=(-pi,pi,-pi/2,pi/2))
>>> savefig('filename.png')
>>>
>>> I get this result: http://dl.dropbox.com/u/8683962/moll-agg1.png
>>>
>>> And if I add aspect=.5 to the imshow call, I get
>>> http://dl.dropbox.com/u/8683962/moll-agg12.png
>>>
>>> Clearly it is not blank in either case, although the second case is
>>> presumably is what you want.
>> Yup, that works for me, but it doesn't apply the Mollweide projection.
>>   I put up a paste here:
>>
>> http://dpaste.com/hold/529372/
>>
>> I run this, and allsky1.png has the image.  allsky2.png and
>> allsky3.png both have empty ovals, and it's not clear to me what's
>> causing that.
>>
>> -Matt
>>
>>>
>>>
>>> On Tue, Apr 5, 2011 at 3:56 PM, Matthew Turk<matthewturk at gmail.com>  wrote:
>>>> 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
>>>> array.
>>>>
>>>> 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')
>>>>
>>>> cb.set_label(r"$\mathrm{Column}\/\mathrm{Density}\/[\mathrm{g}/\mathrm{cm}^2]$")
>>>> canvas = matplotlib.backends.backend_agg.FigureCanvasAgg(fig)
>>>> canvas.print_figure("allsky.png")
>>>>
>>>> (a bit more discussion:
>>>> http://blog.enzotools.org/yt-development-all-sky-column-density-calcula
>>>> )
>>>>
>>>> If anybody knows why the figure size has to be 2:1 to avoid a blank
>>>> image, I'd really appreciate a tip.
>>>>
>>>> -Matt
>>>>
>>>>>         # 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
>>>>>
>>>> _______________________________________________
>>>> AstroPy mailing list
>>>> AstroPy at scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/astropy
>>>>
>>>
>>>
>>> --
>>> Erik Tollerud
>>>
>> _______________________________________________
>> AstroPy mailing list
>> AstroPy at scipy.org
>> http://mail.scipy.org/mailman/listinfo/astropy
>>
>
>




More information about the AstroPy mailing list