Catalogs¶
In this example we have 3 catalogs containing different informations. We want to select data and merge these catalogs in order to plot the size-mass relation for objects with z<0.5. We will see, how to:
- merge 2 catalogs (Santini et al. 2015), selecting only relevant columns.
- cross-match, on the basis of objects coordinates, the merged catalog with a 3rd catalog (see van der Wel et al. 2012).
- select passive galaxies on the basis of object properties and gather galaxies in 3 separated redshift bins.
- select a redshift bin (z<0.5 in this case) and convert objects size in the effective radius in kpc
- plot the size-mass relation
Keys: Table, catalog operations (selection, match, join, sort, ...), units, sky coordinates, cosmology, matplotlib (histogram, scatter plot)
numpy.loadtxt: "low level" function¶
import numpy as np
filename = 'data/CheckPositions.dat'
#Open an homogenus table
table = np.loadtxt(filename)
table[0:5]
!head -2 $filename
#Select columns and define column types and names
subtable = np.loadtxt(filename, usecols=(0, 1, 2, 4),
dtype={'names': ('ID', 'mag', 'expected', 'measured'),
'formats':('i8', 'f4', 'f4', 'f4')})
subtable[0:5]
subtable[0]
dtype¶
A data type object (numpy.dtype) describes how the bytes, in the fixed-size block of memory, corresponding to an array item, should be interpreted. It describes the following aspects:
- type of the data: integer, float, string, etc.
- size of the data: how many bytes are used to store each item (e.g the integer)
- if the data type is a structured data type (e.g. catalog case), an aggregate of other data types:
- what are the names of the “fields” of the structure, by which they can be accessed
- what is the data-type of each field
- which part of the memory block each field takes.
Table¶
astropy.table provides functionality for storing and manipulating heterogeneous tables of data in a way that is familiar to numpy users. A few notable capabilities of this package are:
- Modify a table by adding or removing columns, changing column names, or adding new rows of data.
- Handle tables containing missing values.
- Include table and column metadata as flexible data structures.
- Specify a description, units and output formatting for columns.
- Create a new table by selecting rows or columns from a table.
- Perform Table operations like joins, concatenation, and binning.
- Manipulate multidimensional columns.
- Methods for reading and writing Table objects to files.
Open the phsys parameters catalog¶
from astropy.table import Table
# https://archive.stsci.edu/missions/hlsp/candels/uds/catalogs/v1/
filename = 'data/hlsp_candels_hst_wfc3_uds_santini_v1_physpar_cat.fits'
physpar = Table.read(filename, format='fits')
physpar
Custum open¶
In case we have to open a custom file, not supported by Astropy/Numpy, we can create our own reader.
Open the file and browse it, looping through lines to retrieve the required informations.
In this weird example values are enclosed in ***, and we are interested in just the first value of each row.
!cat data/custom.dat
tab = []
with open('data/custom.dat', 'r') as custom:
for line in custom:
try:
value = line.split('***')[1]
except:
continue
tab.append(int(value))
print(tab)
Download files¶
Astropy provides facility to download files from external sites
from astropy.utils.data import download_file
download_file('https://archive.stsci.edu/missions/hlsp/candels/uds/catalogs/v1/hlsp_candels_hst_wfc3_uds_santini_v1_physpar_cat.fits')
Columns selection¶
columns can not be retrieved by index (it doesn't make sense in this case...).
They can be retrieved by name
physpar[:,1]
physpar['age_2a_tau']
Subcatalog (columns selection)¶
Using the same syntax, we can select a sub catalog (selection by columns in this case)
physpar_sel=physpar[['Seq', 'UMag_6a_tau', 'VMag_6a_tau', 'JMag_6a_tau']]
physpar_sel
Open the masses catalog¶
#Read coordinates and masses values
# https://archive.stsci.edu/missions/hlsp/candels/uds/catalogs/v1/
filename = 'data/hlsp_candels_hst_wfc3_uds_santini_v1_mass_cat.fits'
mass = Table.read(filename, format='fits')
mass
again perform columns selection on this catalog
mass_sel=mass[['Seq', 'RAdeg', 'DECdeg', 'zbest', 'StarFlag', 'AGNFlag', 'M_med']]
mass_sel
from astropy.table import join
joined = join(mass_sel, physpar_sel, keys='Seq', join_type='inner')
joined
from astropy.table import Table, vstack
obs1 = Table.read("data/obs1.dat", format='ascii')
obs1
obs2 = Table.read("data/obs2.dat", format='ascii')
obs2
outer vstack¶
The obs2 table is missing the mag_b column, so in the stacked output table those values are marked as missing
#Join type outer
v1=vstack([obs1, obs2])
v1
v1['mag_b'].data
inner vstack¶
In the case of join_type='inner', only the common columns (the intersection) are present in the output table.
#Join type inner
vstack([obs1, obs2], join_type='inner')
exact vstack¶
In the case of join_type='exact', the vstack requires that all the input tables have exactly the same column names.
vstack([obs1, obs2], join_type='exact') #Error!
obs3 = Table.read("""name obs_date mag_b logLx
M45 2012-02-03 15.0 40.5""", format='ascii')
vstack([obs1, obs3], 'exact')
Join tables¶
The Table class supports the database join operation. This provides a flexible and powerful way to combine tables based on the values in one or more key columns.
from astropy.table import join
optical = Table.read("data/optical.dat", format='ascii')
optical
xray = Table.read("data/xray.dat", format='ascii')
xray
Inner join¶
Represents the strict intersection of the two tables on the key columns
join(optical, xray)
join(optical, xray, keys='name')
Left join¶
The left join creates a new table which has every row from the left table and includes matching values from the right table when available
join(optical, xray, join_type='left')
join(optical, xray, join_type='left', keys='name')
join(optical, xray, join_type='left')
right join¶
Likewise one can construct a new table with every row of the right table and matching left values (when available) using join_type=right.
outer join¶
To make a table with the union of rows from both tables do an outer join
join(optical, xray, join_type='outer')
join(optical, xray, join_type='outer', keys='name')
from astropy.table import unique
obs = Table.read("""name obs_date mag_b mag_v
M31 2012-01-02 17.0 17.5
M82 2012-02-14 16.2 14.5
M101 2012-01-02 15.1 13.5
M31 2012-01-02 17.1 17.4
M101 2012-01-02 15.1 13.5
M82 2012-02-14 16.2 14.5
M31 2012-02-14 16.9 17.3
M82 2012-02-14 15.2 15.5
M101 2012-02-14 15.0 13.6
M82 2012-03-26 15.7 16.5
M101 2012-03-26 15.1 13.5
M101 2012-03-26 14.8 14.3""", format='ascii')
unique(obs, keys='name')
unique(obs, keys=['name', 'obs_date'])
Sort¶
obs.sort(['name'])
obs
obs.sort(['name', 'mag_b'])
obs
NOTE: unlike the other methods, sort works in place (it modifies the table itself)
... back to the execise¶
joined
Rename columns¶
joined.rename_column('M_med', 'M')
joined.rename_column('UMag_6a_tau', 'U')
joined.rename_column('VMag_6a_tau', 'V')
joined.rename_column('JMag_6a_tau', 'J')
joined.rename_column('zbest', 'z')
joined
Add units¶
from astropy import units as u
joined['DECdeg'].unit=u.deg
joined['RAdeg'].unit=u.deg
joined['U'].unit=u.ABmag
joined['V'].unit=u.ABmag
joined['J'].unit=u.ABmag
joined['M'].unit=u.solMass
joined
Units and quantities¶
astropy.units handles defining, converting between, and performing arithmetic operations with physical quantities, such as meters, seconds, Hz, etc.
It also handles logarithmic units such as magnitude and decibel.
quantities are the combination of a value and a unit.
The easiest way to create a Quantity is to simply multiply or divide a value by one of the built-in units.
It works with scalars, sequences and Numpy arrays
v = 23563.7 * u.kilometer / (1.8 * u.second)
v
v.value
v.unit
x=2.7*u.hour*v
x
x=x.decompose()
x
x=x.to(u.lightyear)
x
no_star = (joined['StarFlag'] == 0)
no_agn = (joined['AGNFlag'] == 0)
just_galaxies = no_star & no_agn
selection = joined[just_galaxies]
Syntax more compact
uds=joined[(joined['StarFlag'] == 0) & (joined['AGNFlag'] == 0)]
len(uds)
numpy.where¶
numpy.where returns elements numpy.where(condition[, x, y]) chosen from x or y depending on condition
a=np.arange(10)
np.where(a<5, 2*a, 3.14)
np.where(a<5)
a<5
NOTE: when only condition is provided, this function is a shorthand for np.nonzero. Using nonzero directly should be preferred, as it behaves correctly for subclasses.
np.nonzero(a<5)
Remove columns¶
star and ANG flags are no more required now
uds.remove_columns(['StarFlag', 'AGNFlag'])
Create and add new columns¶
The selection of the passive galaxies will be performed on the V-J/U-V color/color plane. We want to add these values to the catalog
#Compute colors (add columns like arrays)
uv = uds['U'] - uds['V']
vj = uds['V'] - uds['J']
#Add columns to table
uds.add_column(uv, name='U-V', index=4)
uds.add_column(vj, name='V-J', index=5)
uv = uds['U'] - uds['V']
vj
Define passive galaxies by redfhit bins¶
- create a plot grid with 3 different redshift bins: $z<0.5$, $0.5\leq z<1.0$, $z\ge1.0$
#%matplotlib qt
%matplotlib notebook
#Import moudules
import matplotlib.pyplot as plt
How to display data: matplotlib¶
- matplotlib is a Python plotting library
- produces publication quality figures
- can be used in Python scripts, Python/IPython shells, the Jupyter notebook, ...
matplotlib.pyplotprovides a MATLAB-like plotting framework.import pylabis no longer recommended
%matplotlib¶
Set up matplotlib to work interactively
This function lets you activate matplotlib interactive support at any point during an IPython session. It does not import anything into the interactive namespace.
In this case, where the matplotlib default is TkAgg: %matplotlib
Using matplotlib backend: TkAgg
But you can explicitly request a different GUI backend: %matplotlib qt
You can list the available backends using the -l/–list option: %matplotlib --list
The plot grid¶
Create a plot grid a 3 different redshift bins: $z<0.5$, $0.5\le z<1.0$, $z\ge1.0$
#Add prepare the layout 2x2
fig = plt.figure("color vs color")
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222, sharex=ax1, sharey=ax1)
ax3 = fig.add_subplot(223, sharex=ax1, sharey=ax1)
ax1.set_ylabel('U-V')
ax3.set_ylabel('U-V')
ax2.set_xlabel('V-J')
ax3.set_xlabel('V-J')
Plot z selection¶
zmask1 = uds['z'] < 0.5
gal1 = uds[zmask1]
ax1.scatter(gal1['V-J'], gal1['U-V'], alpha=0.1)
zmask2= (0.5<= uds['z']) & (uds['z'] <1.0)
gal2 = uds[zmask2]
ax2.scatter(gal2['V-J'], gal2['U-V'], alpha=0.1)
zmask3= uds['z'] >=1
gal3 = uds[zmask3]
ax3.scatter(gal3['V-J'], gal3['U-V'], alpha=0.1)
fig.show() #no interactive mode
gal1
"#Select all passive galaxies, redshift independet
cmask = (uds['U-V'] > 1.3) & (uds['V-J'] < 1.6)
# z < 0.5 passive galaxies
pmask1 = (uds['U-V']>0.88*uds['V-J']+0.69) & cmask & zmask1
pas1 = uds[pmask1]
ax1.scatter(pas1['V-J'], pas1['U-V'], alpha=0.1, color='orange')
# 0.5 <= z < 1.0 passive galaxies
pmask2 = (uds['U-V']>0.88*uds['V-J']+0.59) & cmask & zmask2
pas2 = uds[pmask2]
ax2.scatter(pas2['V-J'], pas2['U-V'], alpha=0.1, color='orange')
# z >= 1.0 all galaxies
pmask3 = (uds['U-V']>0.88*uds['V-J']+0.49) & cmask & zmask3
pas3 = uds[pmask3]
ax3.scatter(pas3['V-J'], pas3['U-V'], alpha=0.1, color='orange')
plt.show()
plt.close(fig)
Pretty plot for papers¶
fig = plt.figure("paper plot")
fig.set_figwidth(13.0)
fig.set_figheight(9.0)
fig.subplots_adjust(hspace=0, wspace=0)
#0<z<0.5
ax1=fig.add_subplot(221, ylabel='$U-V_{rest}$')
ax1.xaxis.set_visible(False)
ax1.text(2, 0.25, '$0\leq z_{best}\leq 0.5$', horizontalalignment='left', verticalalignment='bottom')
ax1.minorticks_on()
ax1.set_xlim([0.,2.6])
ax1.set_ylim([0.,2.9])
ax1.scatter(gal1['V-J'], gal1['U-V'], alpha=0.1, color='blue', label='all galaxies')
ax1.scatter(pas1['V-J'], pas1['U-V'], alpha=0.1, color='orange', label='passive galaxies')
#ax1.legend()
ax2=fig.add_subplot(222, xlabel='$V-J_{rest}$')
ax2.yaxis.set_visible(False)
ax2.text(2, 0.25, '$0.5 < z_{best} \leq 1.0$', horizontalalignment='left', verticalalignment='bottom')
ax2.minorticks_on()
ax2.set_xlim([0.,2.6])
ax2.set_ylim([0.,2.9])
ax2.scatter(gal2['V-J'], gal2['U-V'], alpha=0.1, color='blue', label='all galaxies')
ax2.scatter(pas2['V-J'], pas2['U-V'], alpha=0.1, color='orange', label='passive galaxies')
#ax2.legend()
ax3=fig.add_subplot(223, xlabel='$V-J_{rest}$', ylabel='$U-V_{rest}$')
ax3.text(2, 0.25, '$z_{best} > 1.0$', horizontalalignment='left', verticalalignment='bottom')
ax3.minorticks_on()
ax3.set_xlim([0.,2.6])
ax3.set_ylim([0.,2.9])
ax3.scatter(gal3['V-J'], gal3['U-V'], alpha=0.1, color='blue', label='all galaxies')
ax3.scatter(pas3['V-J'], pas3['U-V'], alpha=0.1, color='orange', label='passive galaxies')
leg3=ax3.legend()
for lh in leg3.legendHandles: lh.set_alpha(1)
fig.savefig('tmp/paper_plot.png')
plt.close(fig)
Add the passive column¶
passive = Table.Column(data=(pmask1 | pmask2 | pmask3), name="passive")
uds.add_column(passive)
save the table¶
- add the extetion name
uds.meta['EXTNAME']='UDS GAL'
- write the file
uds.write('tmp/uds_galaxies.dat', format='ascii', overwrite=True)
uds.write('tmp/uds_galaxies.fits', format='fits', overwrite=True) #ABmag raise an error
BinTableHDU.from_columns does not convert ABmag
for col in uds.columns:
if uds[col].unit == u.ABmag:
uds[col].unit = u.mag
uds.write('tmp/uds_galaxies.fits', format='fits', overwrite=True)
The 3rd catalog: morphological parameters¶
filename = 'data/uds_3dhst.v4.1_f160w.galfit'
# http://www.mpia-hd.mpg.de/homes/vdwel/allfields.tar
galfit = Table.read(filename, format='ascii')
galfit
from astropy.coordinates import SkyCoord
galfit_coords=SkyCoord(galfit['RA']*u.deg, galfit['DEC']*u.deg)
uds_coords=SkyCoord(uds['RAdeg'], uds['DECdeg'])
idx, d2d, _ = uds_coords.match_to_catalog_sky(galfit_coords)
idx: Indexes of galfit_coords objectd2d: 2d separation of the matches- the galfit sub catalog matching with uds catalog
match = galfit[idx]
match = galfit[idx]
Create new columns¶
- the matching ids column
new_col1 = Table.Column(idx, name="match_indexs")
- the separation in match
new_col2 = Table.Column(data=d2d.arcsec, name="separation", unit=u.arcsec)
Add columns to UDS table¶
match["re"].unit=u.arcsec
uds.add_columns((match["re"], match["q"], new_col1, new_col2))
uds
Check separation values¶
fig = plt.figure("separation")
plt.clf()
_ = plt.hist(uds["separation"], color='red', bins=50)
plt.clf()
_ = plt.hist(uds["separation"][uds["separation"]<0.5], color='orange', bins=50)
plt.close(fig)
Final selection¶
good separation
flags = uds["separation"] < 0.3passive
flags&= uds["passive"]suitable z bin
flags&= uds["z"]<=0.5suitable in M
flags&= uds["M"]>10**10and suitable re
flags&= uds["re"]>0.03
flags = uds["separation"] < 0.3
flags&= uds["passive"]
flags&= uds["z"]<=0.5
flags&= uds["M"]>10**10
flags&= uds["re"]>0.03
good=uds[flags]
good
cosmologies¶
The astropy.cosmology subpackage contains classes for representing cosmologies, and utility functions for calculating commonly used quantities that depend on a cosmological model: distances, ages and lookback times corresponding to a measured redshift or the transverse separation corresponding to a measured angular separation.
cosmology.parameters.available
Compute distances¶
Here WMAP9 is a built-in object describing a cosmology with the parameters from the 9-year WMAP results
from astropy.cosmology import parameters
parameters.available
from astropy.cosmology import WMAP9
re_arcmin = good["re"].to(u.arcmin)*np.sqrt(good["q"])
re_kpc = WMAP9.kpc_proper_per_arcmin(good["z"])*re_arcmin
#Clear the figure
fig = plt.figure()
plt.clf()
plt.title("Size mass relation")
plt.errorbar(good['M'], re_kpc, fmt='o')#re_kpc)
#Set axis scales
plt.semilogx()
plt.semilogy()
#Set axis labels
plt.xlabel("mass ($M_{\odot}$)")
plt.ylabel("size (kpc)")
#Display dashed grid
plt.grid(which="both", ls='--')
Useless return parameters¶
Syntax used to ignore returned function parameters
lines=["2012-01-02 M31 17.0 16.0 16.3", "2012-10-29 M82 16.2 15.2", "2012-10-31 M101"]
for line in lines:
_, name, *_ = line.split()
print(name)
pandas¶
pandas is an open source, BSD-licensed library providing high-performance, easy-to-use data structures and data analysis tools for the Python programming language.
Python has long been great for data munging and preparation, but less so for data analysis and modeling. pandas helps fill this gap, enabling you to carry out your entire data analysis workflow in Python without having to switch to a more domain specific language like R.
pandas=physpar.to_pandas()
pandas
astroquery¶
Astroquery is an Astropy affiliated package that contains a collection of tools to access online Astronomical data. Each web service has its own sub-package
Importing astroquery on its own doesn't get you much: you need to import each sub-module specifically
$ conda install -c astropy astroquery
Query results¶
We can obtain Tables instances
from astroquery.esasky import ESASky
# Available catalogs
ESASky.list_catalogs()
#List of tables
results=ESASky.query_object_catalogs("M31","all")
results[2]
or list of HDUImages
from astroquery.ukidss import Ukidss
Ukidss.get_images(SkyCoord(49.489*u.deg, -0.27*u.deg, frame='galactic'), image_width=5 * u.arcmin)