SELECT TOP 100000 * FROM public.gaia_source WHERE CONTAINS(POINT('ICRS',public.gaia_source.ra,public.gaia_source.dec),CIRCLE('ICRS',64.1171208,28.1266139,1.0))=1
In [1]:
import pandas as pd
import numpy as np
import seaborn as sns
%config InlineBackend.figure_format = 'retina'
%matplotlib inline
In [2]:
import matplotlib.pyplot as plt
In [3]:
from astropy.io import votable
In [4]:
out = votable.parse_single_table('../data/result.vot_01.gz')
In [5]:
out2 = out.to_table()
In [6]:
out3 = out2.to_pandas()
In [7]:
out2[0:5]
Out[7]:
In [10]:
plt.plot(out3.ra, out3.dec, '.', alpha=0.1)
Out[10]:
In [29]:
gaia = out3
In [30]:
gpos = SkyCoord(gaia.ra*u.deg, gaia.dec*u.deg, frame='icrs')
In [11]:
from astropy.nddata import Cutout2D
from astropy import units as u
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
sns.set_style('ticks')
In [13]:
fn = '../data/LkCa4_DSS.fits'
hdu_list = fits.open(fn)
hdu_list.info()
In [14]:
image_data = hdu_list[0].data
hdu_list.close()
wcs = WCS('../data/LkCa4_DSS.fits')
In [15]:
position = SkyCoord('04h16m28.109s +28d07m35.81s', frame='icrs')
In [25]:
cutout = Cutout2D(image_data, position, 50.0*u.arcsecond, wcs=wcs)
cutout2 = Cutout2D(image_data, position, 150.0*u.arcsecond, wcs=wcs)
In [26]:
plt.imshow(image_data, origin='lower', interpolation='none')
cutout.plot_on_original(color='red', lw=2)
Out[26]:
In [27]:
from matplotlib.patches import Rectangle
In [52]:
ra_vec = np.arange(-5, 5)*u.arcsecond
d_ra = ra_vec.to(u.degree).value
mult_ra = position.ra.value + d_ra
dec_vec = np.arange(-5, 5)*u.arcsecond
d_dec = dec_vec.to(u.degree).value
mult_dec = position.dec.value + d_dec
fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(131, projection=wcs)
plt.imshow(image_data.data, origin='lower', cmap='cubehelix_r', interpolation='none')
plt.xlabel('RA')
plt.ylabel('Dec')
cutout.plot_on_original(color='red', lw=2)
lon1 = ax1.coords[0]
lon1.set_major_formatter('hh:mm:ss')
#lon1.set_separator((':', ":"))
ax1.set_title('DSS overview')
xlims, ylims = ax1.get_xlim(), ax1.get_ylim()
ax1.set_xlim(xlims)
ax1.set_ylim(ylims)
ax1.scatter([gpos.ra.value], [gpos.dec.value], transform=ax1.get_transform('icrs'), s=20, c='g', marker='o', lw=1)
ax2 = fig.add_subplot(132, projection=cutout.wcs)
plt.imshow(cutout.data, origin='lower', cmap='Greys', interpolation='none')
lon = ax2.coords[0]
lat = ax2.coords[1]
lon.set_ticklabel_visible(False)
lat.set_ticklabel_visible(False)
lon.set_ticks(mult_ra * u.degree)
lat.set_ticks(mult_dec * u.degree)
ax2.scatter([position.ra.value], [position.dec.value], transform=ax2.get_transform('icrs'), s=60, c='r', marker='+', lw=2)
ra1 = position.ra.value + (12.0*u.arcsecond).to(u.degree).value
dec1 = position.dec.value + (10.0*u.arcsecond).to(u.degree).value
rect_post = (ra1, dec1)
ras, decs = (7.8*u.arcsecond).to(u.degree).value, (7.8*u.arcsecond).to(u.degree).value
r = Rectangle(rect_post, ras, decs, edgecolor='blue', facecolor='none', lw=2,
transform=ax2.get_transform('icrs'))
#ax2.add_patch(r)
ax2.set_title(r"DSS 50''$\times$50'' ")
xlims, ylims = ax2.get_xlim(), ax2.get_ylim()
ax2.set_xlim(xlims)
ax2.set_ylim(ylims)
ax2.scatter([gpos.ra.value], [gpos.dec.value], transform=ax2.get_transform('icrs'), s=20, c='g', marker='o', lw=1)
ax3 = fig.add_subplot(133, projection=cutout2.wcs)
plt.imshow(cutout2.data, origin='lower', cmap='Greys', interpolation='none')
lon3 = ax3.coords[0]
lat3 = ax3.coords[1]
lon3.set_ticklabel_visible(False)
lat3.set_ticklabel_visible(False)
lon3.set_ticks(mult_ra * u.degree)
lat3.set_ticks(mult_dec * u.degree)
ax3.scatter([position.ra.value], [position.dec.value], transform=ax3.get_transform('icrs'), s=60, c='r', marker='+', lw=2)
xlims, ylims = ax3.get_xlim(), ax3.get_ylim()
ax3.set_xlim(xlims)
ax3.set_ylim(ylims)
ax3.scatter([gpos.ra.value], [gpos.dec.value], transform=ax3.get_transform('icrs'), s=20, c='g', marker='o', lw=1)
ax3.set_title(r"DSS 150''$\times$150'' ")
plt.savefig('../results/LkCa4_DSS_gaia.png', dpi=300)
In [51]:
!mkdir ../results
In [53]:
ra_vec = np.arange(-5, 5)*u.arcsecond
d_ra = ra_vec.to(u.degree).value
mult_ra = position.ra.value + d_ra
dec_vec = np.arange(-5, 5)*u.arcsecond
d_dec = dec_vec.to(u.degree).value
mult_dec = position.dec.value + d_dec
fig = plt.figure(figsize=(10, 5))
ax1 = fig.add_subplot(131, projection=wcs)
plt.imshow(image_data.data, origin='lower', cmap='cubehelix_r', interpolation='none')
plt.xlabel('RA')
plt.ylabel('Dec')
cutout.plot_on_original(color='red', lw=2)
lon1 = ax1.coords[0]
lon1.set_major_formatter('hh:mm:ss')
#lon1.set_separator((':', ":"))
ax1.set_title('DSS overview')
xlims, ylims = ax1.get_xlim(), ax1.get_ylim()
ax1.set_xlim(xlims)
ax1.set_ylim(ylims)
#ax1.scatter([gpos.ra.value], [gpos.dec.value], transform=ax1.get_transform('icrs'), s=20, c='g', marker='o', lw=1)
ax2 = fig.add_subplot(132, projection=cutout.wcs)
plt.imshow(cutout.data, origin='lower', cmap='Greys', interpolation='none')
lon = ax2.coords[0]
lat = ax2.coords[1]
lon.set_ticklabel_visible(False)
lat.set_ticklabel_visible(False)
lon.set_ticks(mult_ra * u.degree)
lat.set_ticks(mult_dec * u.degree)
ax2.scatter([position.ra.value], [position.dec.value], transform=ax2.get_transform('icrs'), s=60, c='r', marker='+', lw=2)
ra1 = position.ra.value + (12.0*u.arcsecond).to(u.degree).value
dec1 = position.dec.value + (10.0*u.arcsecond).to(u.degree).value
rect_post = (ra1, dec1)
ras, decs = (7.8*u.arcsecond).to(u.degree).value, (7.8*u.arcsecond).to(u.degree).value
r = Rectangle(rect_post, ras, decs, edgecolor='blue', facecolor='none', lw=2,
transform=ax2.get_transform('icrs'))
#ax2.add_patch(r)
ax2.set_title(r"DSS 50''$\times$50'' ")
xlims, ylims = ax2.get_xlim(), ax2.get_ylim()
ax2.set_xlim(xlims)
ax2.set_ylim(ylims)
#ax2.scatter([gpos.ra.value], [gpos.dec.value], transform=ax2.get_transform('icrs'), s=20, c='g', marker='o', lw=1)
ax3 = fig.add_subplot(133, projection=cutout2.wcs)
plt.imshow(cutout2.data, origin='lower', cmap='Greys', interpolation='none')
lon3 = ax3.coords[0]
lat3 = ax3.coords[1]
lon3.set_ticklabel_visible(False)
lat3.set_ticklabel_visible(False)
lon3.set_ticks(mult_ra * u.degree)
lat3.set_ticks(mult_dec * u.degree)
ax3.scatter([position.ra.value], [position.dec.value], transform=ax3.get_transform('icrs'), s=60, c='r', marker='+', lw=2)
xlims, ylims = ax3.get_xlim(), ax3.get_ylim()
ax3.set_xlim(xlims)
ax3.set_ylim(ylims)
ax3.scatter([gpos.ra.value[vec]], [gpos.dec.value[vec]], transform=ax3.get_transform('icrs'), s=50, c='r', marker='o', lw=1)
ax3.set_title(r"DSS 150''$\times$150'' ")
plt.savefig('../results/LkCa4_DSS_reg.png', dpi=300)
In [47]:
vec = out3.parallax == out3.parallax
vec.sum(), len(vec)
Out[47]:
In [ ]:
In [ ]: