Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Produce a RA vs DEC equatorial coordinates plot with python

I'm trying to generate an equatorial coordinates plot that should look more or less like this one:

enter image description here

(The figure is taken from this article, and it shows the position of the Large and Small MCs in equatorial coordinates)

Important things to notice about this plot:

  • The theta axis (ie: the right ascension) is in h:m:s (hours, minutes, seconds) as it is accustomed in astronomy, rather than in degrees as the default polar option does in matplotlib.
  • The r axis (ie: the declination) increases outward from -90º and the grid is centered in (0h, -90º).
  • The plot is clipped, meaning only a portion of it shows as opposed to the entire circle (as matplotlib does by default).

Using the polar=True option in matplotlib, the closest plot I've managed to produce is this (MWE below, data file here; some points are not present compared to the image above since the data file is a bit smaller):

enter image description here

I also need to add a third column of data to the plot, which is why I add a colorbar and color each point accordingly to a z array:

enter image description here

So what I mostly need right now is a way to clip the plot. Based mostly on this question and this example @cphlewis came quite close with his answer, but several things are still missing (mentioned in his answer).

Any help and/or pointers with this issue will be greatly appreciated.


MWE

(Notice I use gridspec to position the subplot because I need to generate several of these in the same output image file)

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

def skip_comments(f):
    '''
    Read lines that DO NOT start with a # symbol.
    '''
    for line in f:
        if not line.strip().startswith('#'):
            yield line

def get_data_bb():
    '''RA, DEC data file.
    '''

    # Path to data file.
    out_file = 'bb_cat.dat'

    # Read data file
    with open(out_file) as f:
        ra, dec = [], []

        for line in skip_comments(f):
            ra.append(float(line.split()[0]))
            dec.append(float(line.split()[1]))

    return ra, dec

# Read RA, DEC data from file.
ra, dec = get_data_bb()
# Convert RA from decimal degrees to radians.
ra = [x / 180.0 * 3.141593 for x in ra]

# Make plot.
fig = plt.figure(figsize=(20, 20))
gs = gridspec.GridSpec(4, 2)
# Position plot in figure using gridspec.
ax = plt.subplot(gs[0], polar=True)
ax.set_ylim(-90, -55)

# Set x,y ticks
angs = np.array([330., 345., 0., 15., 30., 45., 60., 75., 90., 105., 120.])
plt.xticks(angs * np.pi / 180., fontsize=8)
plt.yticks(np.arange(-80, -59, 10), fontsize=8)
ax.set_rlabel_position(120)
ax.set_xticklabels(['$22^h$', '$23^h$', '$0^h$', '$1^h$', '$2^h$', '$3^h$',
    '$4^h$', '$5^h$', '$6^h$', '$7^h$', '$8^h$'], fontsize=10)
ax.set_yticklabels(['$-80^{\circ}$', '$-70^{\circ}$', '$-60^{\circ}$'],
    fontsize=10)

# Plot points.
ax.scatter(ra, dec, marker='o', c='k', s=1, lw=0.)

# Use this block to generate colored points with a colorbar.
#cm = plt.cm.get_cmap('RdYlBu_r')
#z = np.random.random((len(ra), 1))  # RGB values
#SC = ax.scatter(ra, dec, marker='o', c=z, s=10, lw=0., cmap=cm)
# Colorbar
#cbar = plt.colorbar(SC, shrink=1., pad=0.05)
#cbar.ax.tick_params(labelsize=8)
#cbar.set_label('colorbar', fontsize=8)

# Output png file.
fig.tight_layout()
plt.savefig(ra_dec_plot.png', dpi=300)
like image 588
Gabriel Avatar asked Apr 08 '15 21:04

Gabriel


2 Answers

Getting the colorbar can be done with a merging of the OP code with @cphlewis's excellent answer. I've posted this as a turnkey solution on the request of the OP in chat. The first version of code simply adds a color bar, the final version (under EDIT 2) does an axes affine translation and corrects a few parameters / simplifies the code to suit OP spec exactly.

"""
An experimental support for curvilinear grid.
"""
import numpy as np
import  mpl_toolkits.axisartist.angle_helper as angle_helper
import matplotlib.cm as cmap
from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D

from mpl_toolkits.axisartist import SubplotHost

from mpl_toolkits.axisartist import GridHelperCurveLinear


def curvelinear_test2(fig):
    """
    polar projection, but in a rectangular box.
    """
    global ax1

    # see demo_curvelinear_grid.py for details
    tr = Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform()

    extreme_finder = angle_helper.ExtremeFinderCycle(10, 60,
                                                     lon_cycle = 360,
                                                     lat_cycle = None,
                                                     lon_minmax = None,
                                                     lat_minmax = (0, np.inf),
                                                     )

    grid_locator1 = angle_helper.LocatorHMS(12) #changes theta gridline count
    tick_formatter1 = angle_helper.FormatterHMS()

    grid_locator2 = angle_helper.LocatorDMS(6)
    tick_formatter2 = angle_helper.FormatterDMS()

    grid_helper = GridHelperCurveLinear(tr,
                                        extreme_finder=extreme_finder,
                                        grid_locator1=grid_locator1,
                                        tick_formatter1=tick_formatter1,
                                        grid_locator2=grid_locator2,
                                        tick_formatter2=tick_formatter2
                                        )


    ax1 = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)

    # make ticklabels of right and top axis visible.
    ax1.axis["right"].major_ticklabels.set_visible(True)
    ax1.axis["top"].major_ticklabels.set_visible(True)
    ax1.axis["bottom"].major_ticklabels.set_visible(True) #Turn off? 
    # let right and bottom axis show ticklabels for 1st coordinate (angle)
    ax1.axis["right"].get_helper().nth_coord_ticks=0
    ax1.axis["bottom"].get_helper().nth_coord_ticks=0



    fig.add_subplot(ax1)

    grid_helper = ax1.get_grid_helper()

    ax1.set_aspect(1.)
    ax1.set_xlim(-4,15) # moves the origin left-right in ax1
    ax1.set_ylim(-3, 20) # moves the origin up-down

    ax1.set_ylabel('90$^\circ$ + Declination')
    ax1.set_xlabel('Ascension')
    ax1.grid(True)
    #ax1.grid(linestyle='--', which='x') # either keyword applies to both
    #ax1.grid(linestyle=':', which='y')  # sets of gridlines

    return tr

import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(5, 5))
fig.clf()

tr = curvelinear_test2(fig) # tr.transform_point((x, 0)) is always (0,0)
                            # => (theta, r) in but (r, theta) out...
r_test =   [0, 1.2, 2.8, 3.8, 5,  8,  10, 13.3, 17]  # distance from origin
deg_test = [0,  -7, 12,  28,  45, 70, 79, 90,   100] # degrees ascension
out_test = tr.transform(zip(deg_test, r_test))

sizes = [40, 30, 10, 30, 80, 33, 12, 48, 45]
#hues = [.9, .3, .2, .8, .6, .1, .4, .5,.7] # Oddly, floats-to-colormap worked for a while.
hues = np.random.random((9,3)) #RGB values

# Use this block to generate colored points with a colorbar.
cm = plt.cm.get_cmap('RdYlBu_r')
z = np.random.random((len(r_test), 1))  # RGB values

SC = ax1.scatter(out_test[:,0], #ax1 is a global
            out_test[:,1],
            s=sizes,
            c=z,
            cmap=cm,
            zorder=9) #on top of gridlines
            
# Colorbar
cbar = plt.colorbar(SC, shrink=1., pad=0.05)
cbar.ax.tick_params(labelsize=8)
cbar.set_label('colorbar', fontsize=8)


plt.show()

EDIT

Bit of tidying parameters, adding in OP data, removing redundancy yields the following plot. Still need to centre the data on -90 instead of 0 - at the moment this is hacked, but I'm sure curvelinear_test2() can be changed to account for it...

PIcture matching OP desired format

EDIT 2

Following OP comment on intermediate version in this answer, a final version as below gives the plot at the very end of the post - with -90 on the dec axis and subplot demo

"""
An experimental support for curvilinear grid.
"""
import numpy as np
import  mpl_toolkits.axisartist.angle_helper as angle_helper
import matplotlib.cm as cmap
from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D

from mpl_toolkits.axisartist import SubplotHost

from mpl_toolkits.axisartist import GridHelperCurveLinear


def curvelinear_test2(fig, rect=111):
    """
    polar projection, but in a rectangular box.
    """

    # see demo_curvelinear_grid.py for details
    tr = Affine2D().translate(0,90) + Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform()

    extreme_finder = angle_helper.ExtremeFinderCycle(10, 60,
                                                     lon_cycle = 360,
                                                     lat_cycle = None,
                                                     lon_minmax = None,
                                                     lat_minmax = (-90, np.inf),
                                                     )

    grid_locator1 = angle_helper.LocatorHMS(12) #changes theta gridline count
    tick_formatter1 = angle_helper.FormatterHMS()

    grid_helper = GridHelperCurveLinear(tr,
                                        extreme_finder=extreme_finder,
                                        grid_locator1=grid_locator1,
                                        tick_formatter1=tick_formatter1
                                        )


    ax1 = SubplotHost(fig, rect, grid_helper=grid_helper)

    # make ticklabels of right and top axis visible.
    ax1.axis["right"].major_ticklabels.set_visible(True)
    ax1.axis["top"].major_ticklabels.set_visible(True)
    ax1.axis["bottom"].major_ticklabels.set_visible(True) #Turn off? 
    # let right and bottom axis show ticklabels for 1st coordinate (angle)
    ax1.axis["right"].get_helper().nth_coord_ticks=0
    ax1.axis["bottom"].get_helper().nth_coord_ticks=0



    fig.add_subplot(ax1)

    grid_helper = ax1.get_grid_helper()

    # You may or may not need these - they set the view window explicitly rather than using the
    # default as determined by matplotlib with extreme finder.
    ax1.set_aspect(1.)
    ax1.set_xlim(-4,25) # moves the origin left-right in ax1
    ax1.set_ylim(-3, 30) # moves the origin up-down

    ax1.set_ylabel('Declination')
    ax1.set_xlabel('Ascension')
    ax1.grid(True)
    #ax1.grid(linestyle='--', which='x') # either keyword applies to both
    #ax1.grid(linestyle=':', which='y')  # sets of gridlines

    return ax1,tr
    
    
def skip_comments(f):
    '''
    Read lines that DO NOT start with a # symbol.
    '''
    for line in f:
        if not line.strip().startswith('#'):
            yield line
            
def get_data_bb():
    '''RA, DEC data file.
    '''

    # Path to data file.
    out_file = 'bb_cat.dat'

    # Read data file
    with open(out_file) as f:
        ra, dec = [], []

        for line in skip_comments(f):
            ra.append(float(line.split()[0]))
            dec.append(float(line.split()[1]))

    return ra, dec


import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(5, 5))
fig.clf()

ax1, tr = curvelinear_test2(fig,121) # tr.transform_point((x, 0)) is always (0,0)
                            # => (theta, r) in but (r, theta) out...             

# Read RA, DEC data from file.
ra, dec = get_data_bb()
out_test = tr.transform(zip(ra, dec))

# Use this block to generate colored points with a colorbar.
cm = plt.cm.get_cmap('RdYlBu_r')
z = np.random.random((len(ra), 1))  # RGB values

SC = ax1.scatter(out_test[:,0], #ax1 is a global
            out_test[:,1],
            marker = 'o',
            c=z,
            cmap=cm,
            lw = 0.,
            zorder=9) #on top of gridlines
            
# Colorbar
cbar = plt.colorbar(SC, shrink=1., pad=0.1)
cbar.ax.tick_params(labelsize=8)
cbar.set_label('colorbar', fontsize=8)

ax2, tr = curvelinear_test2(fig,122) # tr.transform_point((x, 0)) is always (0,0)
                            # => (theta, r) in but (r, theta) out...             

# Read RA, DEC data from file.
ra, dec = get_data_bb()
out_test = tr.transform(zip(ra, dec))

# Use this block to generate colored points with a colorbar.
cm = plt.cm.get_cmap('RdYlBu_r')
z = np.random.random((len(ra), 1))  # RGB values

SC = ax2.scatter(out_test[:,0], #ax1 is a global
            out_test[:,1],
            marker = 'o',
            c=z,
            cmap=cm,
            lw = 0.,
            zorder=9) #on top of gridlines
            
# Colorbar
cbar = plt.colorbar(SC, shrink=1., pad=0.1)
cbar.ax.tick_params(labelsize=8)
cbar.set_label('colorbar', fontsize=8)

plt.show()

Final plot:

Picture showing proper subplot

like image 74
J Richard Snape Avatar answered Oct 23 '22 18:10

J Richard Snape


Chewing on the AxisArtist example is actually pretty promising (this combines two AxisArtist examples -- I wouldn't be surprised if AxisArtist was written with RA plots in mind):

enter image description here

Still to do:

  1. Declination should run from -90 at the origin to 0
  2. Be able to use and add a colorbar
  3. adjust limits if plotting outside them

aesthetic:

  1. Serif font in axis labels
  2. Dashed gridlines for ascension

anything else?

"""
An experimental support for curvilinear grid.
"""
import numpy as np
import  mpl_toolkits.axisartist.angle_helper as angle_helper
import matplotlib.cm as cmap
from matplotlib.projections import PolarAxes
from matplotlib.transforms import Affine2D

from mpl_toolkits.axisartist import SubplotHost

from mpl_toolkits.axisartist import GridHelperCurveLinear


def curvelinear_test2(fig):
    """
    polar projection, but in a rectangular box.
    """
    global ax1

    # see demo_curvelinear_grid.py for details
    tr = Affine2D().scale(np.pi/180., 1.) + PolarAxes.PolarTransform()

    extreme_finder = angle_helper.ExtremeFinderCycle(10, 60,
                                                     lon_cycle = 360,
                                                     lat_cycle = None,
                                                     lon_minmax = None,
                                                     lat_minmax = (0, np.inf),
                                                     )

    grid_locator1 = angle_helper.LocatorHMS(12) #changes theta gridline count
    tick_formatter1 = angle_helper.FormatterHMS()

    grid_locator2 = angle_helper.LocatorDMS(6)
    tick_formatter2 = angle_helper.FormatterDMS()

    grid_helper = GridHelperCurveLinear(tr,
                                        extreme_finder=extreme_finder,
                                        grid_locator1=grid_locator1,
                                        tick_formatter1=tick_formatter1,
                                        grid_locator2=grid_locator2,
                                        tick_formatter2=tick_formatter2
                                        )


    ax1 = SubplotHost(fig, 1, 1, 1, grid_helper=grid_helper)

    # make ticklabels of right and top axis visible.
    ax1.axis["right"].major_ticklabels.set_visible(True)
    ax1.axis["top"].major_ticklabels.set_visible(True)
    ax1.axis["bottom"].major_ticklabels.set_visible(True) #Turn off? 
    # let right and bottom axis show ticklabels for 1st coordinate (angle)
    ax1.axis["right"].get_helper().nth_coord_ticks=0
    ax1.axis["bottom"].get_helper().nth_coord_ticks=0



    fig.add_subplot(ax1)

    grid_helper = ax1.get_grid_helper()

    ax1.set_aspect(1.)
    ax1.set_xlim(-4,15) # moves the origin left-right in ax1
    ax1.set_ylim(-3, 20) # moves the origin up-down

    ax1.set_ylabel('90$^\circ$ + Declination')
    ax1.set_xlabel('Ascension')
    ax1.grid(True)
    #ax1.grid(linestyle='--', which='x') # either keyword applies to both
    #ax1.grid(linestyle=':', which='y')  # sets of gridlines

    return tr

import matplotlib.pyplot as plt
fig = plt.figure(1, figsize=(5, 5))
fig.clf()

tr = curvelinear_test2(fig) # tr.transform_point((x, 0)) is always (0,0)
                            # => (theta, r) in but (r, theta) out...
r_test =   [0, 1.2, 2.8, 3.8, 5,  8,  10, 13.3, 17]  # distance from origin
deg_test = [0,  -7, 12,  28,  45, 70, 79, 90,   100] # degrees ascension
out_test = tr.transform(zip(deg_test, r_test))

sizes = [40, 30, 10, 30, 80, 33, 12, 48, 45]
#hues = [.9, .3, .2, .8, .6, .1, .4, .5,.7] # Oddly, floats-to-colormap worked for a while.
hues = np.random.random((9,3)) #RGB values

ax1.scatter(out_test[:,0], #ax1 is a global
            out_test[:,1],
            s=sizes,
            c=hues,
            #cmap=cmap.RdYlBu_r,
            zorder=9) #on top of gridlines


plt.show()
like image 7
cphlewis Avatar answered Oct 23 '22 18:10

cphlewis