In [1]:
import numpy as np
import pandas as pd
from astropy import coordinates
from astropy.coordinates import match_coordinates_sky
import astropy.units as u
import astroquery
from astroquery.irsa import Irsa
from astroquery.vizier import Vizier
from astropy.table import Table, join
Irsa.ROW_LIMIT = -1
Vizier.ROW_LIMIT = -1
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
obj = ["PKS J0006-0623", 1.55789, -6.39315, 1.0]
# name, ra, dec, radius of cone
obj_name = obj[0]
obj_ra = obj[1]
obj_dec = obj[2]
cone_radius = obj[3]
In [3]:
obj_coord = coordinates.SkyCoord(ra=obj_ra, dec=obj_dec, unit=(u.deg, u.deg), frame="icrs")
In [4]:
data_2mass = Irsa.query_region(obj_coord, catalog="fp_psc", radius=cone_radius * u.deg)
data_wise = Irsa.query_region(obj_coord, catalog="allwise_p3as_psd", radius=cone_radius * u.deg)
__data_galex = Vizier.query_region(obj_coord, catalog='II/335', radius=cone_radius * u.deg)
data_galex = __data_galex[0]
In [5]:
num_2mass = len(data_2mass)
num_wise = len(data_wise)
num_galex = len(data_galex)
print("Number of object in (2MASS, WISE, GALEX): ", num_2mass, num_wise, num_galex)
In [6]:
ra_2mass = data_2mass['ra']
dec_2mass = data_2mass['dec']
c_2mass = coordinates.SkyCoord(ra=ra_2mass, dec=dec_2mass, unit=(u.deg, u.deg), frame="icrs")
ra_wise = data_wise['ra']
dec_wise = data_wise['dec']
c_wise = coordinates.SkyCoord(ra=ra_wise, dec=dec_wise, unit=(u.deg, u.deg), frame="icrs")
ra_galex = data_galex['RAJ2000']
dec_galex = data_galex['DEJ2000']
c_galex = coordinates.SkyCoord(ra=ra_galex, dec=dec_galex, unit=(u.deg, u.deg), frame="icrs")
In [7]:
sep_min = 1.0 * u.arcsec # minimum separation in arcsec
In [8]:
idx_2mass, idx_wise, d2d, d3d = c_wise.search_around_sky(c_2mass, sep_min)
c_2mass_wise = c_2mass[idx_2mass]
print("Only 2MASS and WISE: ", len(idx_2mass))
In [9]:
# from matching of 2 cats (2MASS and WISE) coordinate
w1 = data_wise[idx_wise]['w1mpro']
j = data_2mass[idx_2mass]['j_m']
w1j = w1-j
cutw1j = -1.7
# match between WISE and 2MASS
data_wise_2mass = data_wise[idx_wise] # WISE dataset
galaxy = data_wise_2mass[w1j<cutw1j] # https://academic.oup.com/mnras/article/448/2/1305/1055284
In [10]:
w1j_galaxy = w1j[w1j<cutw1j]
w1_galaxy = w1[w1j<cutw1j]
j_galaxy = j[w1j<cutw1j]
plt.scatter(w1j, w1, marker='o', color='blue')
plt.scatter(w1j_galaxy, w1_galaxy, marker='.', color="red")
plt.axvline(x=cutw1j) # https://academic.oup.com/mnras/article/448/2/1305/1055284
Out[10]:
In [11]:
# + GALEX
___idx_2mass_wise_galex, idx_galex, d2d, d3d = c_galex.search_around_sky(c_2mass_wise, sep_min)
match_all_coord = c_galex[idx_galex]
print("Match all 3 cats: ", len(match_all_coord))
In [12]:
### GALEX data which match with 2MASS and WISE!
match_galex = data_galex[idx_galex]
match_galex
Out[12]:
In [13]:
idx_2mass, idx, d2d, d3d = match_all_coord.search_around_sky(c_2mass, sep_min)
match_2mass = data_2mass[idx_2mass]
match_2mass
Out[13]:
In [17]:
ra_2mass = match_2mass['ra']
dec_2mass = match_2mass['dec']
c_2mass = coordinates.SkyCoord(ra=ra_2mass, dec=dec_2mass, unit=(u.deg, u.deg), frame="icrs")
idx_wise, idx, d2d, d3d = c_2mass.search_around_sky(c_wise, sep_min)
match_wise = data_wise[idx_wise]
match_wise
Out[17]:
In [18]:
# joindata = Table([match_2mass['j_m'],
# match_2mass['j_m']-match_2mass['h_m'],
# match_2mass['j_m']-match_2mass['k_m'],
# match_2mass['j_m']-match_wise['w1mpro'],
# match_2mass['j_m']-match_wise['w2mpro'],
# match_2mass['j_m']-match_wise['w3mpro'],
# match_2mass['j_m']-match_wise['w4mpro'],
# match_2mass['j_m']-match_galex['NUVmag']],
# names=('J', 'J-H', 'J-K', 'J-W1', 'J-W2', 'J-W3', 'J-W4', 'J-NUV'))
In [19]:
joindata = np.array([match_2mass['j_m'],
match_2mass['j_m']-match_2mass['h_m'],
match_2mass['j_m']-match_2mass['k_m'],
match_2mass['j_m']-match_wise['w1mpro'],
match_2mass['j_m']-match_wise['w2mpro'],
match_2mass['j_m']-match_wise['w3mpro'],
match_2mass['j_m']-match_wise['w4mpro'],
match_2mass['j_m']-match_galex['NUVmag']])
joindata = joindata.T
In [20]:
# joindata = np.array([match_2mass['j_m'], match_2mass['h_m'], match_2mass['k_m'],
# match_wise['w1mpro'], match_wise['w2mpro'], match_wise['w3mpro'], match_wise['w4mpro'],
# match_galex['NUVmag']])
# joindata = joindata.T
In [21]:
len(match_2mass['j_m'])
Out[21]:
In [22]:
len(match_wise['w3mpro'])
Out[22]:
In [23]:
from sklearn import datasets
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
X = scale(joindata)
pca = PCA(n_components=2)
X_r = pca.fit(X).transform(X)
In [24]:
plt.scatter(X_r[:,0], X_r[:,1])
for i, name in enumerate(match_wise['designation']):
for galaxyname in galaxy['designation']:
if name == galaxyname:
plt.scatter(X_r[i,0], X_r[i,1], color="red")
In [32]:
from sklearn.cluster import DBSCAN
from sklearn.preprocessing import StandardScaler
X = scale(joindata)
db = DBSCAN(eps=0.2, min_samples=3).fit(X)
core_samples_mask = np.zeros_like(db.labels_, dtype=bool)
core_samples_mask[db.core_sample_indices_] = True
labels = db.labels_
# Number of clusters in labels, ignoring noise if present.
n_clusters_ = len(set(labels)) - (1 if -1 in labels else 0)
print('Estimated number of clusters: %d' % n_clusters_)
#print(labels)
In [34]:
# Black removed and is used for noise instead.
unique_labels = set(labels)
colors = [plt.cm.Spectral(each) for each in np.linspace(0, 1, len(unique_labels))]
for k, col in zip(unique_labels, colors):
if k == -1:
# Black used for noise.
col = [0, 0, 0, 1]
class_member_mask = (labels == k)
## J vs J-W1
xy = X[class_member_mask & core_samples_mask]
plt.plot(xy[:, 3], xy[:, 0], 'o', markerfacecolor=tuple(col), markeredgecolor='k', markersize=14)
xy = X[class_member_mask & ~core_samples_mask]
plt.plot(xy[:, 3], xy[:, 0], 'o', markerfacecolor=tuple(col), markeredgecolor='k', markersize=8)
for i, name in enumerate(match_wise['designation']):
for galaxyname in galaxy['designation']:
if name == galaxyname:
plt.plot(X[i,3], X[i,0], marker="X", markerfacecolor='red', markeredgecolor='none', markersize=8)
plt.title('Estimated number of clusters: %d' % n_clusters_)
plt.show()
In [35]:
from sklearn.manifold import TSNE
X = scale(joindata)
X_r = TSNE(n_components=2).fit_transform(X)
In [36]:
plt.scatter(X_r[:,0], X_r[:,1], marker='.', color="blue")
for i, name in enumerate(match_wise['designation']):
for galaxyname in galaxy['designation']:
if name == galaxyname:
plt.scatter(X_r[i,0], X_r[i,1], marker='.', color="red")
In [ ]: