#!/usr/bin/env python
# coding: utf-8

# # Astropy tutorial @ ASTRON
# # Using PyVO to explore LOTSS and TGSS data
# 
# Author: [Tammo Jan Dijkema](t.j.dijkema@gmail.com)
# 
# This tutorial first given at ASTRON on 13 March 2019

# In[3]:


import pyvo as vo


# In[24]:


from astropy.coordinates import SkyCoord
import astropy.units as u
import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')


# We will extract some sources from a more or less random cutout of the LoTSS HETDEX field.

# In[41]:


lotss_center = SkyCoord(ra=195.89191*u.deg, dec=49.74803*u.deg)


# The capabilities of all tables on the ASTRON VO-server are listed on https://vo.astron.nl/ . In particular, for LoTSS, the list of fields is documented [here](https://vo.astron.nl/hetdex/lotss-dr1/cone/info).

# In[42]:


lotss_query = vo.dal.scs.SCSQuery('https://vo.astron.nl/hetdex/lotss-dr1/cone/scs.xml')
lotss_query['RA'] = lotss_center.ra.to(u.deg).value
lotss_query['DEC'] = lotss_center.dec.to(u.deg).value
lotss_query.radius = 1.5*u.deg


# Unfortunately `pyvo` does not understand astropy units: the ra and dec need to be given in degrees without unit.

# In[44]:


lotss_table = lotss_query.execute().to_table()


# In[45]:


lotss_table[:5]


# In[46]:


len(lotss_table)


# **Exercise**: Query the same region in TGSS (TIFR GMRT Sky Survey) alternative data release, available from `https://vo.astron.nl/tgssadr/q/cone/scs.xml`, store the result in `tgss_table`

# In[ ]:





# In[49]:


tgss_table[:5]


# In[50]:


len(tgss_table)


# To perform a cross match, we only look at brighter sources in LoTSS:

# In[51]:


lotss_selection = lotss_table[lotss_table['Total_flux']>20]


# Cross matching is done on SkyCoord objects. We make one SkyCoord object for a lot of coordinates:

# In[52]:


lotss_sources = SkyCoord(ra=lotss_selection['RA'], dec=lotss_selection['DEC'])


# In[53]:


tgss_sources = SkyCoord(ra=tgss_table['RA'], dec=tgss_table['DEC'])


# Matching is done with `match_to_catalog_sky`. This gives for every TGSS source the closest LoTSS source, and the distance to this. (A more thorough crossmatch would use either ADQL queries or `search_around_sky`, which gives a many-to-many relationship.

# In[54]:


lotss_idx, dist_2d, _ = tgss_sources.match_to_catalog_sky(lotss_sources)


# In[55]:


lotss_idx


# **Exercise**: Are all these matches good? Have a look at `dist_2d`.

# In[ ]:





# **Exercise**: plot `lotss_selection[lotss_idx]["Total flux"]` against `tgss_table["Sint"]`. Filter the results to use only the points where `dist2d<10*u.arcsec`.

# In[ ]:





# ## Image access

# In[27]:


tgss_service = vo.dal.sia.SIAService('https://vo.astron.nl/tgssadr/q_fits/imgs/siap.xml')


# In[28]:


lotss_service = vo.dal.sia.SIAService('https://vo.astron.nl/hetdex/lotss-dr1-img/imgs/siap.xml')


# In[29]:


res = lotss_service.search(pos=lotss_center, size=0.01)


# In[30]:


res['imageTitle']


# In[31]:


fig, ax = plt.subplots()
for num, cov in enumerate(res['coverage']):
    ax.plot(cov[::2], cov[1::2], label=str(num)+res['imageTitle'][num].decode('utf-8'));
ax.plot([lotss_center.ra.to(u.deg).value], [lotss_center.dec.to(u.deg).value], 'ro');
ax.legend();


# In[32]:


import astropy.io.fits as fits


# In[33]:


res['accref'][2].decode('utf-8')


# In[34]:


f = fits.open(res['accref'][2].decode('utf-8'))


# In[35]:


fig, ax = plt.subplots(figsize=(10,10))
ax.imshow(f[0].data.squeeze(),
          vmin=0,
          vmax=np.percentile(f[0].data[np.logical_not(np.isnan(f[0].data))], 99.9)
         );


# It would of course be better to add proper axes.

# In[36]:


from astropy.wcs import WCS


# In[37]:


wcs = WCS(f[0])


# In[38]:


wcs


# In[58]:


def plot_lotss(tgss_dots=False, lotss_dots=False, lotss_flux_min=100):
    fig = plt.figure(figsize=(10,10))
    ax = fig.add_subplot(1, 1, 1, projection=wcs)
    ax.imshow(f[0].data.squeeze(),
              vmin=0,
              vmax=np.percentile(f[0].data[np.logical_not(np.isnan(f[0].data))], 99.9)
             );
    if lotss_dots:
        lotss_subset = lotss_table[lotss_table["Total_flux"]>lotss_flux_min]
        ax.plot(lotss_subset["RA"], lotss_subset["DEC"], 'ro', transform=ax.get_transform('world'),
               label="LoTSS catalog > {}mJy".format(lotss_flux_min));            
    if tgss_dots:
        ax.plot(tgss_table["RA"], tgss_table["DEC"], 'y.', transform=ax.get_transform('world'),
               label="TGSS catalog");
    ax.legend()
    return ax


# In[59]:


plot_lotss(tgss_dots=True, lotss_dots=True, lotss_flux_min=10);


# In[102]:


from ipywidgets import interact, interactive, fixed, interact_manual, FloatSlider


# In[103]:


interact(plot_lotss, lotss_flux_min=FloatSlider(min=0, max=150, step=5, continuous_update=False));