astro-datalab / notebooks-latest

Default set of Data Lab notebooks, by DL team and contributed by users
BSD 3-Clause "New" or "Revised" License
57 stars 48 forks source link

A new NB using PHAT data on M31 stellar pops #209

Closed galaxyumi closed 10 months ago

galaxyumi commented 11 months ago

The NB draft is ready for review. With some confusion in ZP values with photometry, I decided to plot stellar SEDs in magnitude in stead of converting them into rates or fluxes. This could be replaced later when the rate columns in the PHAT tables are properly re-ingested.

galaxyumi commented 10 months ago

Thanks @jacquesalice!! I really appreciate all your feedback. I made updates based on your suggestions/feedback. Please take a look at them and let me know what could be improved further.

rnikutta commented 10 months ago

A very nice NB @galaxyumi , thank you for writing it!

My comments below:

result = qc.query(sql=query) # by default the result is a CSV formatted string

instead of:

try:
    result = qc.query(sql=query) # by default the result is a CSV formatted string
except Exception as e:
    print(e.message)
df = qc.query(sql=query,fmt='pandas')
print("Number of columns:",len(df.columns))
print("List of columns:", df.columns)
df

Before each figure, it could just say "We are now plotting the distribution of MS stars in M31". etc.

WHERE ((f275w_gst=1 AND f336w_gst=1) AND 
       (f475w_gst=1 AND f814w_gst=1) AND 
       (f110w_gst=1 AND f160w_gst=1)) AND
       brick=1

with the simpler

WHERE f275w_gst=1 AND f336w_gst=1 AND 
      f475w_gst=1 AND f814w_gst=1 AND 
      f110w_gst=1 AND f160w_gst=1 AND
      brick=1
def make_cmds(brick,starlist=None,cmap='gray_r'):
    """
    brick: Pandas dataframe for a given PHAT Brick
    starlist: list of indices for stars (default=None)
    """

    if starlist is None:
        cmap = 'viridis'

    fig, (ax1, ax2, ax3) = plt.subplots(1,3,figsize=(18,4))

    def plot_panel(ax,x,y,range_=None,xlabel='',ylabel='',title='',starlist=None):
        h = ax.hist2d(x,y,bins=200,range=range_,cmap=cmap,norm=plt.matplotlib.colors.LogNorm())
        ax.set_xlabel(xlabel,fontsize=15)
        ax.set_ylabel(ylabel,fontsize=15)
        ax.set_xlim(h[1].min()-0.5,h[1].max()+0.5)
        ax.set_ylim(h[2].max()+1,h[2].min()-1)
        ax.set_title(title,fontsize=20)
        if starlist is not None:
            colors = cycle(mcolors.TABLEAU_COLORS)  # recurring cycle of defined plot colors
            ax.scatter(x[starlist], y[starlist], s=50, c=[next(colors) for j in range(len(starlist))])

    plot_panel(ax1,brick['f275w_vega']-brick['f336w_vega'],brick['f336w_vega'],((-2,4),(15,27)),'F275W - F336W','F336W','UV CMD',starlist)
    plot_panel(ax2,brick['f475w_vega']-brick['f814w_vega'],brick['f814w_vega'],((-1,6),(14,29)),'F475W - F814W','F814W','Optical CMD',starlist)
    plot_panel(ax3,brick['f110w_vega']-brick['f160w_vega'],brick['f160w_vega'],((-1,3),(11,28)),'F110W - F160W','F160W','IR CMD',starlist)

    plt.show()

Note that I've used cycle from itertools and colors from matplotlib, to get the same colors here and later in the function make_seds. This requires two extra lines in the Imports cell:

In the # std lib section:

from itertools import cycle

and just after import matplotlib.pyplot as plt:

import matplotlib.colors as mcolors

def make_seds(brick,k=3):
    fig, ax = plt.subplots(1, 1, figsize=(6,4))

    stars_6b = good_stars(brick)
    sIDs = random.choices(stars_6b, k=k)
    c = cycle(mcolors.TABLEAU_COLORS)  # recurring cycle of defined plot colors
    for j,i in enumerate(sIDs):
        ax.plot(plambda, [brick['f275w_vega'][i], brick['f336w_vega'][i], 
                          brick['f475w_vega'][i], brick['f814w_vega'][i],
                          brick['f110w_vega'][i], brick['f160w_vega'][i]], marker='.', ls='-', c=next(c))
    ymin = np.min(brick.iloc[sIDs].values.ravel())
    ymax = np.max(brick.iloc[sIDs].values.ravel())
    ax.set_ylim(ymax+0.5, ymin-0.5)
    ax.set_xlabel(r'$\lambda$ [nm]',fontsize=15)
    ax.set_ylabel('magnitude',fontsize=15)

    make_cmds(brick, sIDs)

    plt.show()

Note also that I made 'k' an argument to make_seds() so that users can specify how many random stars they want to get.

make_seds(df_b1,k=3)

, for instance.

"Brick 15 covers a portion of the 10 kpc star-forming ring of the M31's disk."

to

"Brick 15 covers a portion of the 10 kpc star-forming ring in M31's disk."

galaxyumi commented 10 months ago

@rnikutta Thanks for the comments. I implemented all of your comments and also added an html version. I think this NB is ready to be merged, and could be included in the newsletter if you want!

rnikutta commented 10 months ago

Thank you @galaxyumi , NB looks very good! Could you please make one more commit with these tiny changes, before I merge?

galaxyumi commented 10 months ago

@rnikutta All 4 points were addressed.

rnikutta commented 10 months ago

@galaxyumi Sorry to bother you again, but I noticed that your NB is saved without the outputs. We usually store them with the outputs, unless the outputs make the ipynb file very large (multi-MB). Could you please re-run, save with the outputs, and also generate the HTML preview with the outputs?

This will also retain the per-cell-runtime lines (which we want).

Many thanks!

galaxyumi commented 10 months ago

Done! Thank you.

rnikutta commented 10 months ago

Thank you @galaxyumi for the NB, merging now!