<jupyter_start><jupyter_text>Introduction to Astropy Tables  Astropy tables: columns and units<jupyter_code>import numpy as np
import astropy.units as u
from astropy.table import Table, QTable
from astropy.coordinates import SkyCoord, Angle
import matplotlib.pyplot as plt
a = np.array([1, 4, 5], dtype=np.int32)
b = [2.0, 5.0, 8.5] * u.cm
c = ['x', 'y', 'z']
d = [10, 20, 30] * u.m / u.s

t = QTable([a, b, c, d],
           names=('id', 'length', 'label', 'velocity'),
           meta={'name': 'first table'})
t
t['velocity']
t['velocity'][1]
t['velocity'][0:2]
t[2]
t['velocity'].unit<jupyter_output><empty_output><jupyter_text>Obtaining information and statistics<jupyter_code>t.info
t.info('stats')
stats = t.info('stats', out=None)
stats
stats[3]['mean']<jupyter_output><empty_output><jupyter_text>Time and coordinates<jupyter_code>from astropy.time import Time
from astropy.coordinates import SkyCoord

Time.FORMATS
mjd = Time([56200.25, 56400.33, 57500.66, 58000], format='mjd', scale='utc')
date = mjd.strftime('%Y-%m-%d %H:%M')   # You can also convert to .to_datetime() object

sc = SkyCoord([10, 20, 30, 40], [-45, +40, +55, 33], unit='deg')

tab = QTable([mjd, date,  sc, sc.to_string('hmsdms')], names=['MJD', 'Date', 'skycoord', 'coord'])
tab
tab['coord']<jupyter_output><empty_output><jupyter_text>Managing table columns<jupyter_code>tab.columns
tab.colnames<jupyter_output><empty_output><jupyter_text>Substitute column<jupyter_code>tab['MJD'] = [56100, 56200, 56300, 56400]
tab<jupyter_output><empty_output><jupyter_text>Create new column. Directly or with operations<jupyter_code>tab['x'] = [1,2,3,4]
tab['y'] = [10,20,30,40]
tab
tab['z'] = tab['x']*2 + tab['y']
tab<jupyter_output><empty_output><jupyter_text>Selection and slicing of relevant columns<jupyter_code>tab['x']
type(tab['x'])
tab['x'].data
tab['x'].mean()<jupyter_output><empty_output><jupyter_text>Select multiple columns<jupyter_code>tab[['MJD', 'x', 'z']]

my_columns = ['x', 'y', 'z']
tab[my_columns]<jupyter_output><empty_output><jupyter_text>QTable vs Table ⚠️ `Table` with units have limited functionality. For full unit management, it is recommended to use `QTable`. For most applications, `Table` is enough, but if you have to manage and operate with columns with maximum compatibility with unit management, `QTable` is recommended. The caveat is that `QTables` are much more strict and slightly limited when dealing with other python modules, for example to work with matplotlib you constantly need to use `with quantity_support()`. More info in [Plotting Astropy objects in Matplotlib](https://docs.astropy.org/en/stable/visualization/matplotlib_integration.html).  ⛏ Exercise 2.0- Create a `Table` called `tab1` with one column with values `[np.arange(10)]` named `[velocity]`. Use m/s.- Create a `QTable` called `qtab1` with one column with values `[np.arange(10)]` named `[velocity]`. Use m/s.- Show or print both tables. Do you see any difference?- Compute the mean  value of `tab1['velocity']` using `np.mean()` and also using `.mean()`. Do the same for `qtab1['velocity']`- Compute the square `()**2` of the velocity column, for both tables. What are the differences?- Execute `tab1['velocity'].data`. What do you obtain?- Do you obtain the same with `qtab1['velocity'].data`? Try `np.array(qtab1['velocity'])`<jupyter_code>tab1 = Table([np.arange(10)*u.m/u.s], names=['velocity'])
print(tab1)
qtab1 = QTable([np.arange(10)*u.m/u.s], names=['velocity'])
print(qtab1)
np.mean(tab1['velocity']), tab1['velocity'].mean()
np.mean(qtab1['velocity']), qtab1['velocity'].mean()
tab1['velocity']**2
qtab1['velocity']**2
tab1['velocity'].data
qtab1['velocity'].data
np.array(qtab1['velocity'])<jupyter_output><empty_output><jupyter_text>Data filtering<jupyter_code>selection = tab['y'] > 20
selection
tab[selection]
tab['Date'][selection]
selection2 = (tab['x'] > 0) & (tab['MJD'] < 56400)
selection2
tab[['MJD', 'skycoord', 'x']][selection2]<jupyter_output><empty_output><jupyter_text>Storing tables on diskThe Enhanced Character Separated Values table (`.ecsv`) format is a modified `csv` format that can handle metadata from astronomical tables and still keep the universaility and easy accessibility of `csv` files. It is able to store the data types, physical units and descriptions. It uses YAML, which is easy to read by humans and machines and can be opened with a plain text editor. It is not astronomy-specific. More details can be found in the [following link](https://github.com/astropy/astropy-APEs/blob/main/APE6.rst), including the comparison with other formats.<jupyter_code>tab.write('test_tab.csv', format='ascii.csv', overwrite=True)
tab.write('test_tab.ecsv', format='ascii.ecsv', overwrite=True)
tab.write('test_tab.html', format='jsviewer', overwrite=True)
tab.write('test_tab.tex', format='ascii.latex', overwrite=True)<jupyter_output><empty_output><jupyter_text>All available formats are explained in the documentation `tab.write?`  Load our problem data from last session ⚠️  If you didn't create and save your own `data0` you can use the "official" data0 stored in the backup folder to be able to continue with this tutorial. Note that these commands will don't extract the "official" file if `data0.ecsv` already exists. If you want to discard your file and use the official one, first open a terminal and manually remove your `data0.ecsv`.<jupyter_code>import os
if not os.path.isfile('../data/data0.ecsv'):
    print('Using official data0 file')
    os.system('unzip ../data/backup/data0.ecsv.zip')
    os.system('mv data0.ecsv ../data')
else:
    print('Doing nothing because ../data/data0.ecsv already exists')           
data0 = Table.read('../data/data0.ecsv')
data0.show_in_browser(jsviewer=True)<jupyter_output><empty_output><jupyter_text>⛏ Exercise 2.1- Print the first 5 rows of the table with `data0[0:5]`.- Obtain a list of columns of the table `data0`.- Create a new table `my_table` containing only the columns Right Ascension, Declination, and the two proper motions, and the associated errors for all of them. You should obtain a table with 8 columns.- Convert it to a `QTable` using `my_qtable = QTable(my_table)`. Work with `my_qtable` for the rest of the exercise.- Obtain the description (column name, format, units and description) using `.info()` method.- Print the mean (`.mean()`) R.A. in degress and the median Declination in arcmin.- Print the standard deviation (`.std`) of the R.A. in arcseconds and of the Declination in arcmin.- Compute the minimum and maximum, `.min()`, `.max()`, proper motion in right ascension in units of mas/yr- Compute the minimum and maximum proper motion in declination in units of arcmin/day- Create a new column `pm_error` in `my_qtable` with the total uncertainty in the position that accounts for the quadratic sum (`np.sqrt(()**2 + ()**2)`) of the uncertainty of the proper motion components in units of arcmin per minute.- Create a filter to select stars with `pm_error` < 0.0005 arcsec/yr (remember the units!).  - How many stars comply that condition?  - What is the mean `ra` and mean `de` of that subgroup?<jupyter_code>data0[0:5]
data0.columns
my_columns = ['ra', 'dec', 'pmra', 'pmdec', 'pmra_error', 'pmdec_error']
my_table = data0[my_columns]
my_qtable = QTable(my_table)
my_qtable
my_qtable = QTable(my_table)

print(my_qtable['ra'].mean())
print(my_qtable['dec'].mean().to('arcmin'))
print(my_qtable['ra'].std().to('arcsec'))
print(my_qtable['dec'].std().to('arcmin'))


print(my_qtable['pmra'].min())
print(my_qtable['pmra'].max())

print(my_qtable['pmdec'].min().to(u.arcmin/u.day))
print(my_qtable['pmdec'].max().to(u.arcmin/u.day))

pm_error = np.sqrt(my_qtable['pmra_error']**2 + my_qtable['pmdec_error']**2)
my_qtable['pm_error'] = pm_error.to(u.arcsec/u.yr)
      
selection = my_qtable['pm_error'] < 0.00056 * my_qtable['pm_error'].unit
print(len(my_table[['ra','dec']][selection]))

print('Average values:')
print(my_qtable['ra'][selection].mean())
print(my_qtable['dec'][selection].mean())<jupyter_output><empty_output><jupyter_text>🌪 Additional funRepeat all the previous steps but with `my_table` instead of `my_qtable`. What are the main limitations you find when working with `Table` units?<jupyter_code>my_table.info()
my_table['ra'].mean(), my_table['dec'].mean()
np.std(my_table['ra'].to(u.arcsec)), np.std(my_table['dec'].to(u.arcmin))
np.min(my_table['pmra']), np.max(my_table['pmra'])
np.min(my_table['pmdec'].to(u.arcmin/u.day)), np.max(my_table['pmdec'].to(u.arcmin/u.day))
pm_error = np.sqrt(my_table['pmra_error']**2 + my_table['pmdec_error']**2)
my_table['pm_error'] = pm_error.to(u.arcmin/u.min)
my_table[0:3]
selection = my_table['pm_error'] < 0.00056 

print(len(my_table[['ra','dec']][selection]))

print('Average values:')
print(my_table['ra'][selection].mean())
print(my_table['dec'][selection].mean())<jupyter_output><empty_output><jupyter_text>✨ Exercise 2.2Compute a new column for `data0` with the absolute magnitude in the Gaia filter, g. Call the new column `Mg` and use the formula:  $M_{\rm V} = m + 5log_{\rm 10}p - 10$,  where $M_{\rm V}$ is the absolute magnitude, $m$ is the apparent magnitude and $p$ is the stellar parallax in milliarcseconds.For example:Vega has a parallax $p$ of 0.129 arcsec, and an apparent magnitude $m_{\rm V}$ of 0.03:  $M_{\rm V} = 0.03 + 5log_{\rm 10}129 - 10$ = 0.58The column containing the apparent magnitude in the Gaia band is called `phot_g_mean_mag`, and the parallax is `parallax`. The logarithm is computed as `np.log10()`.You have more details, for instance, here: https://en.wikipedia.org/wiki/Absolute_magnitudeExamples- Print the first few rows of `data0` and make sure the new column is created.- Verify that your solution is right with this code:  `np.abs(data0['Mg'][0] - 7.267679281230688) < 1e-5`<jupyter_code>data0['Mg'] = data0['phot_g_mean_mag'] + 5*np.log10(data0['parallax']) - 10
data0[0:3]
np.abs(data0['Mg'][0] - 7.267679281230688) < 1e-5<jupyter_output><empty_output><jupyter_text>Exploratory analysis Visualization Visualization is much easier when working with a `Table`. Feel free to use a `QTable` instead. In this tutorial we will use `data0`, which was loaded as a `Table`.<jupyter_code>import matplotlib.pyplot as plt
#from astropy.visualization import quantity_support   # If you use QTable, you will need this
plt.plot(data0['ra'],     # X variable
         data0['dec'],    # Y variable
         '.',             # marker, can also be 'o', '+', '>', 's', etc..
         color='k',       # color k:black, g:green, r:red, b:blue, etc
         alpha = 0.2,     # transparency, between 0 and 1
         ms=1);           # marker size<jupyter_output><empty_output><jupyter_text>This is not very informative. Let's create a `figure` and `axes` with `subplots`. Let's make the image bigger and use `scatter` to be able to set a symbol size `s` proportional to the flux of the star. We will also add some axes labels.<jupyter_code>fig, ax = plt.subplots(figsize=(15,10))

# We use scatter here because it allows us to plot points with different sizes/colors
ax.scatter(data0['ra'],
           data0['dec'],
           s = data0['phot_g_mean_flux']/1e5)  # Size of the data points, here proportional to the star flux

ax.set_aspect('equal')

# Here we invert the direction of the right ascension axis
ax.invert_xaxis()

ax.set_xlabel('Right Ascension [deg]')  # Write X axis label
ax.set_ylabel('Declination [deg]');     # Write Y axis label<jupyter_output><empty_output><jupyter_text>We can plot any pair of variables. Here a quick plot without too many configuration for quick lookup.<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0['pmra'], data0['pmdec'], '.k')   # '.k' is used as an abbreviation for the format marker='.', color='k', ls='' (black points without lines)<jupyter_output><empty_output><jupyter_text>Now we manually fix the X and Y limits so we can zoom in the relevant region<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0['pmra'],
        data0['pmdec'],
        '.k')

ax.set_xlim(-80, 60)
ax.set_ylim(-80, 60);<jupyter_output><empty_output><jupyter_text>Also we can use log scale when that is relevant.<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0['phot_g_mean_flux'],
        data0['parallax_error'], '.k', ms=1, alpha=0.5)  # ms is markersize and alpha is the transparency, between 0 and 1

ax.loglog();<jupyter_output><empty_output><jupyter_text>But! Never forget to add labels to identify what you are plotting!<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0['phot_g_mean_flux'],
        data0['parallax_error'], '.k', ms=1, alpha=0.5)

ax.loglog()

ax.set_xlabel('G-band mean flux [e/s]')
ax.set_ylabel('Parallax error [mas]');<jupyter_output><empty_output><jupyter_text>Tracking the axis label and units manually means that it is possible that we make mistakes. We can use the column information instead.<jupyter_code>print(data0['phot_g_mean_flux'].description)
print(data0['parallax_error'].unit)<jupyter_output><empty_output><jupyter_text>⛏ Exercise 2.3- Repeat the plot above (phot_g_mean_flux vs parallax_error) but using the methods `.description` and `.unit` to write the x and y labels automatically.- Second, repeat the plot again but now make the script more generic: start by defining `col1` and `col2` with the names of the relevant columns and use only those variables so you don't need to write the explicit column name more than once:```pythoncol1 = 'phot_g_mean_flux'col2 = 'parallax_error'```<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0['phot_g_mean_flux'],
        data0['parallax_error'], '.k', ms=1, alpha=0.5)

ax.loglog()

ax.set_xlabel(f"{data0['phot_g_mean_flux'].description} [{data0['phot_g_mean_flux'].unit}]")
ax.set_ylabel(f"{data0['parallax_error'].description} [{data0['parallax_error'].unit}]");
col1 = 'phot_g_mean_flux'
col2 = 'parallax_error'

fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0[col1], data0[col2], '.k', ms=1, alpha=0.5)

ax.loglog()

ax.set_xlabel(f"{data0[col1].description} [{data0[col1].unit}]")
ax.set_ylabel(f"{data0[col2].description} [{data0[col2].unit}]");<jupyter_output><empty_output><jupyter_text>⛏ Exercise 2.4- Repeat the plot but now use the columns `phot_g_mean_mag`, `parallax`.- Can you interpret the horizontal line with an overdensity of bright stars?- Compute the average parallax of the stars with magnitude < 10. Define a variable `selection`.- Repeat the plot but add one more `ax.plot()` to show in the X axis `data0[col1][selection]` and in the Y axis `data0[col2][selection]`. Make the symbols larger.- Use the `np.mean` and `np.median` functions to find the mean and median parallax of the selection.Tip: you may not be able to compute the **median** because there are missing values and it produces a masked array. You can use the method `.compressed()` to the column selection to be able to compute `np.median`.<jupyter_code>col1 = 'phot_g_mean_mag'
col2 = 'parallax'

fig, ax = plt.subplots(ncols=1, figsize=(8,6))

ax.plot(data0[col1], data0[col2], '.k', ms=1, alpha=0.5)

ax.loglog()

ax.set_xlabel(f"{data0[col1].description} [{data0[col1].unit}]")
ax.set_ylabel(f"{data0[col2].description} [{data0[col2].unit}]");
col1 = 'phot_g_mean_mag'
col2 = 'parallax'

selection = data0[col1] < 10
data0[col1][selection]

fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0[col1], data0[col2], '.k', ms=1, alpha=0.5)
ax.plot(data0[col1][selection], data0[col2][selection], '.', ms=5, alpha=0.5)

ax.loglog()

ax.set_xlabel(f"{data0[col1].description} [{data0[col1].unit}]")
ax.set_ylabel(f"{data0[col2].description} [{data0[col2].unit}]");

print(np.mean(data0[col2][selection]))
print(np.median(data0[col2][selection].compressed()))<jupyter_output><empty_output><jupyter_text>🌪 Additional funPlot horizontal and vertical lines with ax.hline and ax.vline on the relevant magnitude and parallax computed in the previous exercise. You can use `ax.axvline` and/or `ax.axhline` functions.<jupyter_code>col1 = 'phot_g_mean_mag'
col2 = 'parallax'
mag_lim = 10

selection = data0[col1] < mag_lim
data0[col1][selection]

fig, ax = plt.subplots(figsize=(8,6))

ax.plot(data0[col1], data0[col2], '.k', ms=1, alpha=0.5)
ax.plot(data0[col1][selection], data0[col2][selection], '.', ms=5, alpha=0.5)

ax.loglog()

ax.set_xlabel(f"{data0[col1].description} [{data0[col1].unit}]")
ax.set_ylabel(f"{data0[col2].description} [{data0[col2].unit}]")

ax.axvline(mag_lim, color='k', alpha=0.8)
ax.axhline(np.median(data0[col2][selection].compressed()), color='r', alpha=0.2, lw=5);<jupyter_output><empty_output><jupyter_text>Coordinate plots in sky projections<jupyter_code>fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection="mollweide")

ax.scatter(data0['ra'].to(u.radian),
           data0['dec'].to(u.radian), marker='.', color='r')

ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])
ax.grid(True)<jupyter_output><empty_output><jupyter_text>Let's overplot the galactic plane<jupyter_code>galactic_plane = SkyCoord(l=np.arange(-180, 180),
                          b=np.zeros(360),
                          frame='galactic', unit=u.deg)

galactic_plane_eq = galactic_plane.transform_to('icrs')
gal_ra  = galactic_plane_eq.ra.wrap_at('180d').radian
gal_dec = galactic_plane_eq.dec.radian

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection="mollweide")
ax.scatter(data0['ra'].to(u.radian), data0['dec'].to(u.radian), marker='.', color='r')
plt.plot(gal_ra, gal_dec, 'k.')

ax.set_xticklabels(['14h','16h','18h','20h','22h','0h','2h','4h','6h','8h','10h'])
ax.grid(True)<jupyter_output><empty_output><jupyter_text>🌪 Additional funMake the previous plot but in Galactic coordinates  🌪 Additional funPrepare a plot `pmra` vs `pmdec` including errorbars using the matplotlib function `plt.errorbars`.<jupyter_code>pm = data0[['pmra','pmdec','pmra_error', 'pmdec_error']]

plt.errorbar(pm['pmra'], pm['pmdec'],
             yerr=pm['pmra_error'],
             xerr=pm['pmdec_error'],
             marker='.', ls='', color='k', alpha=0.4)
plt.xlim(-50, 50)
plt.ylim(-50, 50);<jupyter_output><empty_output><jupyter_text>Histogram distributions<jupyter_code>fig, ax = plt.subplots(ncols=1, figsize=(10,8))

ax.hist(data0['parallax'], bins=np.arange(-5, 15, 0.5))

ax.set_xlabel('Parallax [mas]')
ax.set_ylabel('Number of stars');<jupyter_output><empty_output><jupyter_text>⛏ Exercise 2.5Modify the line `ax.set_xlabel('Parallax [mas]')` to automatically find the units of the column being plotted. Do the same we did before when we defined the function `str_label()`.<jupyter_code># Plot the same histogram as before but starting with the generic variable col1. Do not use 'parallax' in any other part
col1 = 'parallax'
fig, ax = plt.subplots(ncols=1, figsize=(10,8))

ax.hist(data0[col1], bins=np.arange(-5, 15, 0.5))

description = data0[col1].description
unit = data0[col1].unit

ax.set_xlabel(f"{description} [{unit}]")
ax.set_ylabel('Number of stars');<jupyter_output><empty_output><jupyter_text>✨ Exercise 2.6First of all, we see that there are negative parallaxes, which don't have physical meaning, but are a consequence of errors in the parallax determination. We can create a second `Table` named `data1` that ignores any negative parallax. We will use this table from now on.- Define the variable `positive_parallaxes` as those entries with parallax greater than 0. `data0['parallax'] > 0`. This is a boolean array with `True` values for rows with positive parallaxes, and `False` otherwise.- Now slice table `data0` by selecting the rows with `True` values: `data0[positive_parallaxes]`. Assign that new table to `data1`.- Print the lenght of `data0` and `data1`. How many stars each one has? Use `len(data0)`, `len(data1)`<jupyter_code>positive_parallaxes = data0['parallax'] > 0
data1 = data0[positive_parallaxes]
print(len(data0))
print(len(data1))<jupyter_output><empty_output><jupyter_text>Now we are going to overlay the `data0` and `data1` distributions<jupyter_code>col1 = 'parallax'

fig, ax = plt.subplots(ncols=1, figsize=(10,8))

ax.hist(data0[col1], bins=np.arange(-5, 15, 0.5), color='grey')
ax.hist(data1[col1], bins=np.arange(-5, 15, 0.5))
    
description = data0[col1].description
unit = data0[col1].unit

ax.set_xlabel(f"{description} [{unit}]")
ax.set_ylabel('Number of stars');<jupyter_output><empty_output><jupyter_text>EquivalenciesLet's work with distances in kpc, which are more familiar to us than parallaxes. We will use an astropy unit transformation as before. However, an angle (mas) cannot be converted to distance (kpc) without knowing how the transformation should occur. We need to parse which equivalency to use to make the transormation.There is a lot of information in section [Equivalencies](https://docs.astropy.org/en/stable/units/equivalencies.html), for example to convert spectral units (nm to Hz) or conversions from wavelength/frequency/energy including doppler effect.<jupyter_code>data1['distance'] = data1['parallax'].to(u.kpc, equivalencies=u.parallax())
data1['distance'].description = 'Distance from Earth'
data1[0:3]<jupyter_output><empty_output><jupyter_text>⛏ Exercise 2.7Compute manually the distance in kpc and check that the transformation has worked. Start by using the numpy array without units `np.array(data1['parallax'])` and compute the distance in kpc  with the formula $d [{\rm kpc}]= \frac{1}{p[{\rm mas}]}$, where $p$ is the parallax. Compute the average of the residual `data1['distance'] - d_kpc`, where `d_kpc` is the manually computed distance in units of `kpc`. What is the average residual?<jupyter_code># We compute the ditance as 1/parallax. When the parallax is in mas, the distance will be in kpc.
p = np.array(data1['parallax']) # Let's get rid of the units because we want to do the conversion manually
d = 1./p

d_kpc = u.Quantity(d, unit='kpc') # Now we set the correct units back

print(f"Astropy distances: {data1['distance'].data}")
print(f"My distances: {d_kpc}")

# We compute how different are the astropy parallax from equivalence and our manually computed distances:
residual = data1['distance'] - d_kpc
print(f"Mean residual {np.nanmean(residual)}")<jupyter_output><empty_output><jupyter_text>✨ Exercise 2.8- Going back to the parallax distribution. Plot again the histogram of the `data[parallax]`, starting with `col1='parallax'`. Make the figure size slightly bigger, and use a smaller bin size. Try different `bins` selection either with `np.arange` (which fixes the step) or `np.linspace` (which fixes the number of steps) until you see some significant structure. Identify, approximately, the range of parallaxes that look interesting to you. (You can use `ax.set_xlim()` to tune the plot range).- Make a similar plot but with `col1='distance'`.<jupyter_code>col1 = 'parallax'

fig, ax = plt.subplots(ncols=1, figsize=(15,10))

ax.hist(data1[col1], bins=np.arange(-0.5, 15, 0.1))

description = data1[col1].description
unit = data1[col1].unit

ax.set_xlabel(f"{description} [{unit}]")
ax.set_ylabel('Number of stars');
col1 = 'distance'

fig, ax = plt.subplots(ncols=1, figsize=(15,10))

ax.hist(data1[col1], np.arange(0, 8, 0.05))

description = data1[col1].description
unit = data1[col1].unit

ax.set_xlabel(f"{description} [{unit}]")
ax.set_ylabel('Number of stars');<jupyter_output><empty_output><jupyter_text>We see a very interesting accummulation of stars at a parallax of approximately 5.2 mas. We can create a filter to select the start in that particular range. We will overplot the distribution of the whole sample and the one of the selected group. Check how many stars we have selected with the filter.<jupyter_code>manual_filter1 = (data1['parallax'] > 5.0*u.mas) & (data1['parallax'] < 5.7*u.mas)

cluster1 = data1[manual_filter1]
cluster1<jupyter_output><empty_output><jupyter_text>There is an alternative way to deal with this problem. Instead of creating a new `Table` we could create a new Boolean column with True/False that indicates if each star is part of the cluster or not. We will not use this new column `cluster`, but it is good so we have all the information together.<jupyter_code>data1['cluster'] = manual_filter1
data1[0:3]<jupyter_output><empty_output><jupyter_text>The get the stars from the cluster we can simply select the rows matching the filter, i.e., selecting the `True` values. We could use `data1[manual_filter1]`, but we do not want to carry the filter around, it is better if the information is embedded in the table.<jupyter_code>data1[data1['cluster']]  # This is exactly the same as cluster1<jupyter_output><empty_output><jupyter_text>We will not use the 'cluster' column. It is created here just for training purposes.  ✨ Exercise 2.9- Use the last two plots with the parallax and the distance distributions. Appart from the `data1`, add now the histograms for the table `cluster1`. You should see the two distributions overlapped.- Set the labels to 'Full sample' and 'Cluster', respectively.- Fine tune the data ranges as needed.- What is the approximate distance of the selected cluster?<jupyter_code>col1 = 'parallax'

fig, ax = plt.subplots(ncols=1, figsize=(15,10))

ax.hist(data1[col1], bins=np.arange(-0.5, 15, 0.1), label='Full sample')
ax.hist(cluster1[col1], bins=np.arange(-0.5, 15, 0.1), label='Cluster')

description = data1[col1].description
unit = data1[col1].unit

ax.set_xlabel(f"{description} [{unit}]")
ax.set_ylabel('Number of stars')

ax.legend();
col1 = 'distance'

fig, ax = plt.subplots(ncols=1, figsize=(15,10))

ax.hist(data1[col1], bins=np.arange(0, 8, 0.05), label='Full sample')
ax.hist(cluster1[col1], bins=np.arange(0, 8, 0.05), label='Cluster')


description = data1[col1].description
unit = data1[col1].unit

ax.set_xlabel(f"{description} [{unit}]")
ax.set_ylabel('Number of stars')

ax.legend();<jupyter_output><empty_output><jupyter_text>✨ Exercise 2.10Create the `phot_g_mean_mag` vs `parallax_error` we did before, but including data from the two tables: `data1` in black and `cluster1` in red.<jupyter_code>col1 = 'phot_g_mean_flux'
col2 = 'parallax_error'

fig, ax = plt.subplots(ncols=1, figsize=(8,6))

ax.plot(data1['phot_g_mean_mag'], data1['parallax'],'.k', ms=1, alpha=0.5)
ax.plot(cluster1['phot_g_mean_mag'], cluster1['parallax'], marker='o', color='r', ms=2, ls='')

ax.loglog()

description1 = data1[col1].description
unit1 = data1[col1].unit
description2 = data1[col2].description
unit2 = data1[col2].unit

ax.set_xlabel(f"{description1} [{unit1}]")
ax.set_ylabel(f"{description2} [{unit2}]");<jupyter_output><empty_output><jupyter_text>We see that the cluster is dominated by bright stars (lower magnitude), and form a sharp cluster in parallax/distance, as we already knew.  Spatial distribution of the clusterWe plot the distribution of start in the sky. First, all the stars in the sample are plotted in grey. The stars of the cluster are plotted in color, with the colorscale representing the distance from the Earth in pc.<jupyter_code>fig, ax = plt.subplots(figsize=(15,10))
ax.set_aspect('equal')

ax.scatter(data1['ra'], data1['dec'], c='gray', s=1, alpha=0.5);
l = ax.scatter(cluster1['ra'], cluster1['dec'], c=cluster1['distance']*1000., s=20);  # We save the variable l to be used in the colorbar

ax.set_xlabel('Right Ascension [deg]')
ax.set_ylabel('Declination [deg]');

# Here we invert the direction of the right ascension axis
ax.invert_xaxis()

# Show the color bar
cb = fig.colorbar(l);
cb.set_label('Distance [kpc]')<jupyter_output><empty_output><jupyter_text>Let's include the photometry information and use a different color for the `cluster1` stars.<jupyter_code>fig, ax = plt.subplots(ncols=1, figsize=(15,10))
ax.scatter(data0['ra'], data0['dec'], s=data0['phot_g_mean_flux']/1e5);
ax.scatter(cluster1['ra'], cluster1['dec'], s=cluster1['phot_g_mean_flux']/1e5);

ax.set_aspect('equal')

# Here we invert the direction of the right ascension axis
ax.invert_xaxis()

ax.set_xlabel('Right Ascension [deg]')
ax.set_ylabel('Declination [deg]');<jupyter_output><empty_output><jupyter_text>🌪 Additional funRepeat the previous plot by making it more general, using generic columns `col1`, `col2`, and automatic labels. There is no apparent pattern of the selected stars, although there seems to be an overdensity at the center, specially in Right Ascension.<jupyter_code>fig, ax = plt.subplots(figsize=(8,6))

ra_range = [np.min(data1['ra']), np.max(data1['ra'])]
                 
                                       
ax.hist(data1['ra'],    bins=np.arange(ra_range[0], ra_range[1], 0.04), label='Full sample')
ax.hist(cluster1['ra'], bins=np.arange(ra_range[0], ra_range[1], 0.12), label='Cluster')

ax.set_xlabel('Right Ascension [deg]')
ax.set_ylabel('Number of stars');

ax.legend();
fig, ax = plt.subplots( figsize=(8,6))

de_range = [np.min(data1['dec']), np.max(data1['dec'])]

ax.hist(data1['dec'],    bins=np.arange(de_range[0], de_range[1], 0.02), label='Full sample')
ax.hist(cluster1['dec'], bins=np.arange(de_range[0], de_range[1], 0.08), label='Cluster')

ax.set_xlabel('Declination [deg]')
ax.set_ylabel('Number of stars');

ax.legend();<jupyter_output><empty_output><jupyter_text>Nothing obvious or much interesting. The overdensity in the centre for R.A. is clear, but in Dec it is not clear, there may be a slope. What it is clear is that the cluster is not easily identified in R.A. and Declination. We will need to explore more variables to characterize the cluster  ✨ Exercise 2.11Let's save the data sets for future use. Like in session 1, we will save them with `format='ascii.ecsv'`- Use method `data1.write()` to save the table in file `'../data/data1.ecsv'` with the correct format.- Use method `cluster1.write()` to save the table in file `'../data/cluster1.ecsv'` with the correct format.<jupyter_code>data1.write('../data/data1.ecsv', format='ascii.ecsv', overwrite=True)
cluster1.write('../data/cluster1.ecsv', format='ascii.ecsv', overwrite=True)<jupyter_output><empty_output><jupyter_text>Pandas tables<jupyter_code>import pandas as pd
df1 = data1.to_pandas()
df1<jupyter_output><empty_output><jupyter_text>Exploring a pandas table<jupyter_code>df1.head()<jupyter_output><empty_output><jupyter_text>The first thing we see is that we have lost the unit information. That is an important problem if we don't track the column operations properly.<jupyter_code>df1.info()
df1.columns
df1['ra']
df1[['ra','dec']]
df1['pmra'].values<jupyter_output><empty_output><jupyter_text>Slicing To slice by index value:<jupyter_code>df1.loc[3]<jupyter_output><empty_output><jupyter_text>To slice by index number:<jupyter_code>df1.iloc[3]<jupyter_output><empty_output><jupyter_text>Can be combined with slices in one or several columns.<jupyter_code>df1['pmra'].iloc[3]
df1.iloc[5:8]
df1[['pmra','pmdec']].iloc[6]
df1[['pmra','pmdec']].iloc[6:10]<jupyter_output><empty_output><jupyter_text>Pandas operations. Aggregate and groupby We can easily do operations to columns. First by aggregating values according to some functions<jupyter_code>df1.aggregate(['sum', 'min'])
df1.aggregate({'ra' : ['mean', 'min', 'max', 'std'],
               'dec' : ['mean', 'min', 'max', 'std'],
               'parallax': 'std'})
df1.describe()
def my_func(row):
    return np.sqrt(row['pmra']**2 + row['pmdec']**2)

df1.apply(lambda row: my_func(row), axis=1)<jupyter_output><empty_output><jupyter_text>NaN values<jupyter_code>df1['dr2_radial_velocity'].isna().sum()
df_rv = df1[['dr2_radial_velocity', 'parallax']].dropna()
df_rv<jupyter_output><empty_output><jupyter_text>Groupby<jupyter_code>df1.groupby('cluster').aggregate('mean')
cluster1['astrometric_matched_transits'].description
df1['astrometric_matched_transits'].describe()<jupyter_output><empty_output><jupyter_text>Let's make many groups according to how many observations a star has.<jupyter_code>df1.groupby('astrometric_matched_transits').aggregate(['mean','min'])[['parallax_error','pmra_error', 'pmdec_error']]
df1.groupby('astrometric_matched_transits').aggregate(['mean','min'])[['parallax_error']].plot()<jupyter_output><empty_output><jupyter_text>Working with datetime series<jupyter_code>N = 100
times = pd.date_range("2021-06-9", periods=N, freq='D')

ts = pd.DataFrame({'v1': np.random.normal(0.5,100,N) + 0.1*np.arange(N)**2,
                   'v2': np.random.normal(0.5,100,N) - 0.1*np.arange(N)**2},
                   index=times)
ts
ts.plot();<jupyter_output><empty_output><jupyter_text>Other plots<jupyter_code>df1.plot.scatter(x='phot_g_mean_flux', y='parallax_error',  marker='.', s=0.5, alpha=0.1)
plt.loglog();
plt.plot(data0['phot_g_mean_flux'], data0['parallax_error'], '.k', ms=1, alpha=0.1)
plt.loglog();
df1[['pmra','pmdec']].plot.hist(bins=np.arange(-80, 60, 3), alpha=0.5);
df1[['phot_g_mean_mag','phot_bp_mean_mag','phot_rp_mean_mag']].plot.box();
df1[['phot_rp_mean_mag', 'phot_bp_mean_mag', 'phot_g_mean_mag']].plot.hist(bins = 80, alpha=0.8);
df1.plot.hexbin(x='ra', y='dec', gridsize=20);
from pandas.plotting import scatter_matrix
scatter_matrix(df1[['pmra_error','pmdec_error','pmra','pmdec']], alpha=0.2, diagonal='kde', figsize=(8,8));<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1500618718498174, 'fraction_numerical': 0.028823227597863158, 'mean_word_length': 5.944875288199539, 'pattern_counts': {'":': 0, '<': 235, '<?xml version=': 0, '>': 231, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '56', 'dir': 'jupyter-structured-clean-dedup', 'id': '172077', 'license': 'MIT', 'n_tokens_mistral': 12152, 'n_tokens_neox': 11595, 'n_words': 3149, 'path': 'tutorials/tutorial_02/pysnack_1_02_solutions.ipynb', 'repo_name': 'spsrc/pysnacks_I-', 'size': '33133'}