# -*- coding: utf-8 -*-
"""Python resources for VO.ipynb

Automatically generated by Colab.

Original file is located at
    https://colab.research.google.com/drive/1NN_hFl463g7Fd0exRwZL2DvMSmWEIRpI

---
# Python resources for the Virtual Observatory
## Some use cases of astroquery, Vizier and PyVO
---

**Authored by** Ricardo Rizzo (ISDEFE, SVO collaborator)

**Revised/approved:** [Spanish Virtual Observatory](https://svo.cab.inta-csic.es)

**Creation date:** October 25, 2025

**Updated:** November, 2025

This tutorial presents an introductory use of Python tools aimed to profit the Virtual Observatory (VO) capabilities.

The tutorial assumes that the user:
- has sufficient knowledge of Python 3.
- is able to manage python notebooks and Google colab platform.
- has a general idea about Virtual Observatory (VO).
- does **not** need previous knowledge about existing catalogs or other VO resources.
---

# Overview: the tools

The tutorial focuses on the use of the packages ```Vizier```, ```astroquery```, and ```PyVO```. These packages are considered a bridge between the VO resources and the powerful and versatile Python language.

It follows a "_learning by doing_" approach and is structured around a series of examples which make use of the necessary packages.

Therefore, it is **NOT** a complete reference to all the capabilities of the packages, but a starting guide to use them. Some external references are included when appropriate.

## Why using VizieR
- VizieR (https://vizier.cds.unistra.fr/) provides the most complete library of published astronomical catalogues (25800+ up to November 2025).
- It is a "live" repository, with verified and enriched data from surveys, catalogs and tables.
- Its content is accessible by multiple online interfaces and desktop applications.
- There are a number of query tools that allow the user to select relevant data tables and to extract records matching a given criteria.
- It has several mirrors worldwide.
- Most of its capabilities are integrated in ```astroquery.vizier```, an affiliated package of ```astroquery```.
- Reference: [https://astroquery.readthedocs.io/en/latest/vizier/vizier.html](https://astroquery.readthedocs.io/en/latest/vizier/vizier.html)

## What about astroquery
- astroquery is **THE** query tool in Python for astronomy.
- It has an excellent management of tables in all their formats, particularly VOTables.
- astroquery includes easy access to VizieR functionalities through its convenience functions, which are added to the astroquery.vizier methods and functions.
- Reference: [https://astroquery.readthedocs.io/en/latest/](https://astroquery.readthedocs.io/en/latest/)

## When to go to PyVO
- PyVO works as an "extension" of astropy (https://www.astropy.org), the most popular collection of software packages written in Python for use in astronomy.
- PyVO includes several methods and functions to facilitate an easy access to the mostly commonly used VO protocols.
- When used in combination with numpy, pandas, matplotlib, etc., other high functionalities are possible.
- Reference: [https://github.com/astropy/pyvo](https://github.com/astropy/pyvo)

# Overview: the protocols
The examples developed in this tutorial are based on a client-server structure.

The communication between the tools (on the client side) and the servers (which provide access to data and metadata on the server side) is done through protocols.

IVOA has defined a number of protocols which are necessary to successfully communicate servers and clients.

In this tutorial we use some of the most popular VO protocols, briefly described below.

## Simple Cone Search protocol (SCS)
This protocol is used to provide results from queries to one or more **tables**, based on a given position on the sky.

The parameters needed are therefore very simple: position (RA, Dec) and a search radius (SR) in degrees around it.

Some packages and utilities may include wrappers to use SCS with other parameters, like a search based on a source name.

Reference: [https://www.ivoa.net/documents/latest/ConeSearch.html](https://www.ivoa.net/documents/latest/ConeSearch.html)

## Simple Image Access protocol (SIA)
This protocol provides capabilities for the discovery, description, access, and retrieval of multi-dimensional **image** datasets, including 2D images as well as datacubes of three or more dimensions.

The queries using this protocols need, _as a minimum_, the same parameters as the SCS. In most cases, some of the fields in the metadata are also searchable, which adds flexibility (and complexity!) to the queries.

The resulting VOTables are lists of metadata on images on the queried service/s.

**Reference:** [https://www.ivoa.net/documents/SIA/20151223/index.html](https://www.ivoa.net/documents/SIA/20151223/index.html)

## Simple Spectral Access protocol (SSA)

This protocol works very similar to SIA. Results from a query are VOTable whose metadata (some of them searchable) describe information about **spectra** which can be analised and downloaded.

The resulting spectra from a SSA query are VOTables, which can later be converted to other formats.

**Reference:** [https://www.ivoa.net/documents/SSA/](https://www.ivoa.net/documents/SSA/)

## Table Access Protocol (TAP)

This protocol also allow the client to access tables and catalogs, but
with more powerful capabilities in comparison with the precedent ones.

Its robustness relies in highly flexible specifications to the query, which are derived from the implementation of **ADQL** [(Astronomical Data Query Language)](https://www.ivoa.net/documents/ADQL/20180112/PR-ADQL-2.1-20180112.html).

ADQL is a tailored query language based on SQL. The client can consult to the server adding more conditions to the query. In this way, the queries are faster than using the previous protocols because filtering is made on the server side.

**Reference:** [https://www.ivoa.net/documents/TAP/](https://www.ivoa.net/documents/TAP/)

---
# Install the necessary packages

The packages needed for this tutorial are indicated in the file ```requirements.txt```.
"""

# Download the list of needed packages from the Zenodo tutorial
! wget -qO requirements.txt https://zenodo.org/records/17533386/files/requirements.txt?download=1
# Load requirements.txt from owner's drive
! pip install -r requirements.txt

"""The following packages are necessary:
```
astropy
astroquery
pyvo
numpy
pandas
matplotlib
rich
requests
```
Most likely, the colab session will only install ```pyvo``` and ```astroquery```. The other packages are included here in order to guarantee a correct run in case of exporting the notebook to another environment.

# Using VizieR in Python: a cone search

Below we develop two examples which perform a Cone Search query to VizieR using astroquery.

## Packages, definitions, setups
Before doing the queries, we import the necessary packages, define the target source and create the query instance.
"""

# Import necessary modules
from astroquery.vizier import Vizier
from astropy.coordinates import SkyCoord
import astropy.units as u
from astropy.table import Table
from rich import print

# Define the target (play and change to your favorite source)
TARGET_NAME = "NGC2359"  # a Wolf Rayet Nebula

# Clear the local cache
Vizier.clear_cache()

# Create a Vizier instance
vizier = Vizier(columns=['*'])  # Retrieve all columns

# Optionally, it is possible to select a mirror
us_vizier = Vizier(vizier_server="vizier.cfa.harvard.edu/viz-bin/VizieR-2")
# [[See help(Vizier) for other parameters]]

"""> **⚠️ Warning**

> As it is specified, the instance ```vizier``` will do the query in all the catalogs (more than 20k!). This may take several seconds, or even minutes, to finish. The user may want to specify the catalogs to be searched by specifying them explicitly or by taking advantage of the [CDS categories](https://vizier.cds.unistra.fr/vizier/doc/catstd-2.htx?utm_source=chatgpt.com) (for example, adding ```, catalog="II/*"``` to the instance in case of looking for photometry surveys).

## Example 1: query the available catalogs around the target
| Key function/method  | Description |
|---------------|-------------|
| ```vizier.query_object``` | Queries the service by internally resolving coordinates of the input source name.
||Returns a table object.

In this example we query all the catalogs available in VizieR around a source, with a search radius of 5 arcsec. The search radius is an optional parameter, whose default is 2 arcmin. Note that ```radius``` should be an astropy ```Quantity```.
"""

# Do the query
SEARCH_RADIUS = 5 * u.arcsec
catalogs = vizier.query_object(TARGET_NAME, radius=SEARCH_RADIUS)
print(f"Number of catalogs found around {TARGET_NAME}: {len(catalogs)}")

# Print the catalog codes and number of entries
for catalog in catalogs:
    print(f"Catalog: {catalog.meta['name']}, Number of entries: {len(catalog)}")
# [[Note the difference w.r.t. pyvo cases]]

"""## Example 2: perform a Cone Search around the target coordinates
| Key function/method  | Description |
|---------------|-------------|
| ```vizier.query_region``` | Queries the service in a region centered at some coordinates and with a given radius.
||Returns a table object.
"""

# Do the query
target_coords = SkyCoord.from_name(TARGET_NAME)
SEARCH_RADIUS = 5 * u.arcsec
catalogs = vizier.query_region(target_coords, radius=SEARCH_RADIUS)
if catalogs != None:
    print(f"Number of catalogs found in Cone Search around {TARGET_NAME}: {len(catalogs)}")
else:
    print(f"No catalogs found in Cone Search around {TARGET_NAME}.")

# Print the catalog codes and number of entries
if catalogs != None:
    for catalog in catalogs:
        print(f"Catalog: {catalog.meta['name']}, Number of entries: {len(catalog)}")

"""# More VizieR: Query keywords, table downloading and writing

In this example we do a query by specifying descriptors which will match the keywords provided for the catalogs.

Additionally, we select, download and inspect one of the catalogs. Finally the catalog is written as a VOTable.

## Example 3: work with catalogs containing data related to pulsar rotation

| Key function/method  | Description |
|----------------------|-------------|
| ```vizier.find_catalogs``` | Queries all catalogs matching some keywords.
| ```vizier.get_catalogs```  | Retrieves a list of catalogs.
| ```table.write```          | Writes a table in several formats.

Find all the catalogues matching the keywords.
"""

# Clear the local cache
Vizier.clear_cache()

# Create a Vizier instance
vizier = Vizier(columns=['*'])  # Retrieve all columns

# Search for catalogs related to pulsar rotation
SEARCH_TERM = 'pulsars rotation' # Play and change it by your favourite topic!
pulsar_list = vizier.find_catalogs(SEARCH_TERM)

# Number of catalogs found
print("-------------------------------------------------------------")
print(f"Searching for catalogs with term: '{SEARCH_TERM}'")
print(f"Number of catalogs found: {len(pulsar_list)}")

# Print catalog codes and their descriptions
print("-------------------------------------------------------------")
for k,v in pulsar_list.items():
   print(k, ":", v.description)

"""Select one of the catalogs and inspect their tables."""

# Get a specific catalog by its code
CATALOG_CODE = 'J/A+A/684/A208'

# Limit the number of rows to retrieve (optional)
vizier.ROW_LIMIT = 1 # Set to the minimum for faster retrieval

# Retrieve the catalog
catalog = vizier.get_catalogs(CATALOG_CODE) # Retrieving may take some time
print("-------------------------------------------------------------")
print(f"Retrieving catalog: {CATALOG_CODE}")
print(f"Number of tables in catalog: {len(catalog)}")
print(catalog)

"""Work with the first table of the catalog."""

# Access to the first table (still with just one row)
table = catalog[0]
print("-------------------------------------------------------------")
print(f"First table in catalog '{CATALOG_CODE}':")
print(table)

# Access to the columns
print("-------------------------------------------------------------")
print("Columns in the table:")
print(table.colnames)

# Select a specific column and print its description
print("-------------------------------------------------------------")
print(f"Description of column 'FeRASS':\n{table['FeRASS'].description}")

"""Access to the data in the table (download and write as VOTable)."""

# Access to data (all rows)
vizier.ROW_LIMIT = -1 # All available rows
#OPTION vizier.TIMEOUT = 500  # Increase timeout to 500 sec
catalog = vizier.get_catalogs(CATALOG_CODE)
table = catalog[0]
print("-------------------------------------------------------------")
print(f"Number of rows: {len(table)}")
print(f"Data in selected table, from catalog {CATALOG_CODE}:")
print(table)

# Write table to a local VOtable
outfile = f"{SEARCH_TERM.replace(' ','_')}.vot"
table.write(outfile, format='votable', overwrite=True)
print("-------------------------------------------------------------")
print(f"Catalog saved to local file/s {outfile}")

"""# Life beyond VizieR: a Cone Search using PyVO

In the following example we mine the Gaia DR3 catalog around the nearby Pleiades star cluster.

Through basic table manipulation, we construct a color-magnitude diagram, plot the proper motion vectors and select candidate Pleiades members based on them.

> **📝 Note:**

> _This science case will be fully developed in this School using Topcat._

## Packages, service, coordinates

Before doing the query, we import the necessary packages, define the service and get the target coordinates.
"""

# Import packages
import pyvo
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from astropy.table import Table
from astropy import units as u
from astropy.coordinates import SkyCoord, get_icrs_coordinates

# Definition of service (protocol and URL)
ACCESS_URL = "https://gea.esac.esa.int/tap-server/conesearch?TABLE=gaiadr3.gaia_source&IDCOL=source_id&RACOL=ra&DECCOL=dec&"
service = pyvo.dal.SCSService(ACCESS_URL)

# Find coordinates of a source from its name
pleiades_coord = get_icrs_coordinates('M45')

# This convenience function may also be used
pleiades_coord = SkyCoord.from_name('M45')
print(f"Pleiades coordinates (ICRS): RA={pleiades_coord.ra.degree}, Dec={pleiades_coord.dec.degree}")

"""## Example 4: query the Gaia DR3 main catalog and manipulate results

| Key function/method  | Description |
|----------------------|-------------|
| ```pyvo.dal.SCSService``` | Definition of a Cone Search service. |
| ```service.search``` | Do a Cone Search query around a target position. |
| ```pyvo.conesearch``` | Convenience function to avoid defining the service. |

Do the query and print the number of stars found.
"""

# Do the query
target_ra = pleiades_coord.ra.degree
target_dec = pleiades_coord.dec.degree
radius = 1.    # also in [deg]
results = service.search((target_ra, target_dec), size=(radius, radius))

# This convenience function do the search in a circle:
#results = pyvo.conesearch(ACCESS_URL, pos=[target_ra, target_dec], radius=radius)

# Results works similarly to an astropy Table
print(f"Number of results: {len(results)}")

"""Table manipulation and conversion."""

# Access to fields/columns and data
results_table = results.to_table()
print(results_table.colnames)  # List of column names
print("\n--------------------------")
print(results_table['source_id'])  # Access to a specific column

# Convert to VOTable
results_vot = results.votable

# Print all the fields (i.e. columns)
for field in results_vot.iter_fields_and_params():
    print(f"Field/Param: {field.name}")

"""Save results as a VOTable."""

# Save results to a local VOTable file
results_vot.to_xml('pleiades_gaia_dr3.vot')

"""Create a color-magnitude diagram."""

# Convert to pandas DataFrame
df = results.to_table().to_pandas()
print(df.head())

# Plot 1: G magnitude vs. BP-RP color
plt.figure(figsize=(8,6))
plt.scatter(df['bp_rp'], df['phot_g_mean_mag'], s=1, color='blue')
plt.gca().invert_yaxis()
plt.xlabel('BP - RP [mag]')
plt.ylabel('G mag [mag]')
plt.title('Color-Magnitude Diagram of Pleiades (Gaia DR3). N={}'.format(len(df)))
plt.savefig('pleiades_cmd.pdf', dpi=300)
plt.show()

"""Plot proper motion vectors in the sky."""

# Plot 2: Proper motion vector field (using quiver with scaling)
plt.figure(figsize=(8,8))
scale_factor = 3e-4  # Adjust this factor to scale the arrows
plt.quiver(df['ra'], df['dec'], df['pmra']*scale_factor, df['pmdec']*scale_factor, angles='xy', scale_units='xy', scale=1, color='red', width=0.002)
plt.xlabel('RA [deg]')
plt.ylabel('Dec [deg]')
plt.gca().invert_xaxis()
plt.title('Proper Motion Vector Field of Pleiades (Gaia DR3). N={}'.format(len(df)))
plt.axis('equal')
plt.show()

"""Now we do a zoom around the center of Pleiades to better distinguish the vectors."""

# Plot 2 repeated: a zoom around the center of Pleiades
fig, ax = plt.subplots(figsize=(8, 8))
scale_factor = 3e-4  # Scale factor for proper motion arrows

# Conversion: 1 arcmin = 1/60 deg
half_size_deg = (10 / 60) / 2   # 10 arcmin total → ±5 arcmin around the center

# Plot the proper motion vectors
ax.quiver(df['ra'], df['dec'],
          df['pmra'] * scale_factor, df['pmdec'] * scale_factor,
          angles='xy', scale_units='xy', scale=1, color='red', width=0.002)

# Labels and title
ax.set_xlabel('RA [deg]')
ax.set_ylabel('Dec [deg]')
ax.set_title(f'Proper Motion Vector Field of Pleiades (10×10 arcmin). N={len(df)}')

# Define the field center
ra_center = pleiades_coord.ra.degree
dec_center = pleiades_coord.dec.degree

# Set axis limits (RA decreases to the right, Dec increases upward)
ax.set_xlim(ra_center + half_size_deg, ra_center - half_size_deg)  # RA inverted
ax.set_ylim(dec_center - half_size_deg, dec_center + half_size_deg)

# Keep the aspect ratio square (so the field is not distorted)
ax.set_aspect('equal')

plt.show()

"""Apparently, there is a number of stars with higher proper motions towards the southeast.

Let's see a diagram of proper motion _values_ in RA, Dec.

_Note the group of points around (20, -45) mas/yr._
"""

# Plot 3: Proper motions in RA vs. Dec
plt.figure(figsize=(8,6))
plt.scatter(df['pmra'], df['pmdec'], s=1, color='blue')
plt.xlabel('pmra [mas/yr]')
plt.ylabel('pmdec [mas/yr]')
plt.xlim(-50, 70)
plt.ylim(-80, 40)
plt.title('Proper Motions of Pleiades (Gaia DR3). N={}'.format(len(df)))
plt.grid()
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.savefig('pleiades_pm.pdf', dpi=300)
plt.show()

"""Create a group of comoving sources around (20, -45) mas/yr."""

# Select comoving sources
pm_mask = (
    (df['pmra'] > 15) &
    (df['pmra'] < 25) &
    (df['pmdec'] > -50) &
    (df['pmdec'] < -40)
)

comoving = df[pm_mask]
print(f"Number of comoving sources: {len(comoving)}")
print(comoving[['source_id', 'ra', 'dec', 'pmra', 'pmdec', 'phot_g_mean_mag']])

"""Repeat the color-magnitude diagram for the whole sample and the comoving stars."""

# Repeat plot 1 (a CMD) for df and comoving sources
plt.figure(figsize=(8,6))
plt.scatter(df['bp_rp'], df['phot_g_mean_mag'], s=1, color='blue', label="All sources")
plt.scatter(comoving['bp_rp'], comoving['phot_g_mean_mag'], s=1, color='red', label="Comoving sources")
plt.gca().invert_yaxis()
plt.xlabel('BP - RP [mag]')
plt.ylabel('G mag [mag]')
plt.title('Color-Magnitude Diagram of Comoving Sources in Pleiades (Gaia DR3). N={}'.format(len(comoving)))
plt.savefig('pleiades_comoving_cmd.pdf', dpi=300)
plt.show()

"""# More PyVO: access to images

In this example we use pyVO to perform a Simple Image Access (SIA) query.

## Example 5: search and download an image

| Key function/method  | Description |
|----------------------|-------------|
| ```pyvo.regsearch``` | Do a query consulting the VO Registry. |
| ```pyvo.dat.SIAService``` | Definition of a Simple Image Access query |

Import the necessary packages
"""

# Import packages
import pyvo
from pyvo.dal.sia import SIAResults
import requests
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from astropy.table import Table
from astropy import units as u
from astropy.coordinates import SkyCoord, get_icrs_coordinates
from astropy.io import fits
from astropy.wcs import WCS
from astropy.io.votable import from_table
import subprocess

"""Define a function to plot an image from a 2D FITS file."""

# Function to plot a 2D FITS image

def plot_FITS(filename):
    from pathlib import Path
    from astropy.io import fits
    from astropy.wcs import WCS
    import matplotlib.pyplot as plt
    from astropy.visualization import ZScaleInterval, ImageNormalize

    # Open FITS file
    hdu = fits.open(filename)[0]
    data = hdu.data
    header = hdu.header

    # Get WCS information
    wcs = WCS(header)

    # Compute zscale normalization (automatic intensity range)
    interval = ZScaleInterval()
    vmin, vmax = interval.get_limits(data)
    norm = ImageNormalize(vmin=vmin, vmax=vmax)

    # Plot the FITS image with WCS projection
    fig = plt.figure()
    ax = fig.add_subplot(111, projection=wcs)
    im = ax.imshow(data, origin='lower', cmap='viridis', norm=norm)

    # Add coordinate grid and labels
    #ax.coords.grid(color='cyan', ls='--', lw=0.7, alpha=0.7)
    ax.set_xlabel('Right Ascension')
    ax.set_ylabel('Declination')

    # Custom formatter for the bottom-right status bar
    def format_coord(x, y):
        try:
            ra, dec = wcs.wcs_pix2world([[x, y]], 0)[0]
            return f"RA={ra:.6f}°  Dec={dec:.6f}°   Pixel=({x:.1f},{y:.1f})"
        except Exception:
            return f"Pixel=({x:.1f},{y:.1f})"

    ax.format_coord = format_coord

    plt.colorbar(im, ax=ax, orientation='vertical', label='Intensity')
    plt.title(Path(filename).name)
    plt.savefig(filename.replace('fits','pdf'))
    print(f"Image saved to {filename.replace('fits','pdf')}")
    plt.show()

"""Look for available services in the VO Registry."""

# It is possible to look for all the available (i.e. registered) SIA services
services = pyvo.regsearch(servicetype='sia', keywords=['image'])
print(f"Number of SIA services found: {len(services)}")

# Too many results. Let's refine the search to the collections hosted by NOIRLab Data Lab
services = pyvo.regsearch(servicetype='sia', keywords=['noirlab'])
print(f"Number of SIA services found for NOIRLab: {len(services)}")

for pos, row in enumerate(services.to_table()):
    print(f"\nService {pos}: {row['res_title']}")
    if 'access_urls' in row.colnames:
        print(f"URL: {row['access_urls']}")
    else:
        print("No access_url found.")

"""> **⚠️ Warning**

> The field found as ```access_urls``` in colab may appear as ```access_url``` if you run the code in a terminal.

> _Take it into account in your daily work._

Select the "coadd images" service (10th in the list)
"""

# Select coadded images from NOIRLab Data Lab
coadd_service = services.to_table()[9]
print(f"Selected service: {coadd_service['res_title']}")
SERVICE_URL = coadd_service['access_urls']
service = pyvo.dal.SIAService(SERVICE_URL)
print(f"URL: {SERVICE_URL}")

"""Find the coordinates and do the query."""

# Find coordinates of a source from its name
TARGET = 'Horsehead Nebula'
target_coords = get_icrs_coordinates(TARGET)
# Convenience:
target_coords = SkyCoord.from_name(TARGET)
target_ra = target_coords.ra.degree
target_dec = target_coords.dec.degree
radius = 0.25*u.deg    # radius may be specified using astropy units
print("\n-------------------------------------------------------------")
print(f"Target coordinates: (RA, Dec)=({target_coords.ra.degree:5f}, Dec={target_coords.dec.degree:5f}). Search radius: {radius}")

# Do the query
images = service.search(target_coords, size=(radius,radius))

# results works similarly to an astropy Table
print(f"Number of images: {len(images)}")

"""Inspect the results."""

# Inspect the column names
images_table = images.to_table()
print(images_table.colnames)  # List of column names

# Note that not all the 'prodtype' are 'image'.
prodtypes = set(images_table['prodtype'])
print(f"Product types in the results: {prodtypes}")

"""Convert to VOTable and save results."""

# Convert to VOTable
images_vot = images.votable

# Print all the fields (i.e. columns)
for field in images_vot.iter_fields_and_params():
    print(f"Field/Param: {field.name}")

# Save query results to a local VOTable file
images_table.write('images.vot', format='votable', overwrite=True)
print("\n-------------------------------------------------------------")
print("Results saved to images.vot")

"""Download and plot one of the images."""

# Access and download one image
horsehead = images[21]
image_url = horsehead.getdataurl()
print(f"Selected image data URL: {image_url}")
response = requests.get(image_url)
with open('horsehead.fits', 'wb') as f:
    f.write(response.content)
print("\n-------------------------------------------------------------")
print("Selected image downloaded as horsehead.fits")

plot_FITS('horsehead.fits')

"""# Diving deep inside VO: the Registry
Most of the Virtual Observatory resources are available through the VO Registry: a compilation of metadata and applications to access them.

In this example we do a query of spectroscopic data by a direct request to the Registry.

## Example 6: explore PyVO capabilites searching the Registry for spectra.

| Key function/method  | Description |
|----------------------|-------------|
| ```pyvo.regsearch``` | Convenience function to do a query to the VO Registry. |
| ```pyvo.dat.SSAService``` | Definition of a Simple Spectral Access query |

Import the necessary packages
"""

import pyvo
from pyvo import registry
from rich import print
import requests
import matplotlib.pyplot as plt
from astropy.table import Table
from astropy.coordinates import SkyCoord

"""Start by searching all the catalogs available through the Simple Spectral Access (SSA) protocol."""

# Search for Simple Spectral Access (SSA) services using pyVO regsearch convenience function
services = pyvo.regsearch(servicetype='ssa')
print(f"Number of Simple Spectral Access services found: {len(services)}")

"""Too many results.

Let's refine the search to services related to LAMOST
"""

# Search only to LAMOST using PyVO registry method
ssa_lamost_services = registry.search(servicetype="ssa", keywords="LAMOST")
type(ssa_lamost_services)

# Print number of services found
print(f"Number of SSA services related to LAMOST: {len(ssa_lamost_services)}")

"""Look at some details about the services found."""

# Title of the services and URLs to access them
pos = 0
for s in ssa_lamost_services:
    print(f"Service {pos}: {s['res_title']}")
    try:
        print(f"Access URLs: {s['access_urls']}")
    except KeyError:
        print("No access_urls found.")
    pos += 1

"""Select the last service (LAMOST DR6)"""

# Select the last service (LAMOST DR6)
lamost_dr6 = ssa_lamost_services[-1]
print(f"Selected LAMOST DR6 MRS service title: {lamost_dr6['res_title']}")

for col,val in lamost_dr6.items():
    print(f"{col}: {val}")

"""Build the service and do the query"""

# Get the URL from the previous request
# query_url = lamost_dr6['access_urls']
query_url = lamost_dr6['access_urls'][0] # Extract the first URL from the list
print(f"Query URL: {query_url}")

# Build the service object
lamost_service = pyvo.dal.SSAService(query_url)
print(f"SSAService object created: {lamost_service}")

# Performing the query
spectra = lamost_service.search(pos=(180.0, 0.0), diameter=3.)

# Alternatively, try with a source name
#pos = SkyCoord.from_name('VY CMa')
#spectra = lamost_service.search(pos=pos, diameter=3.)
print(f"Number of results from LAMOST DR6 MRS query: {len(spectra)}")

"""Manipulate the results and save the table"""

# Convert results to an Astropy Table
spectra_table = spectra.to_table()
print(spectra_table)

# Convert to VOTable and write to file
spectra_vot = spectra.votable
spectra_vot.to_xml('lamost_dr6_spectra.vot')
print("\n-------------------------------------------------------------")
print("Results saved to lamost_dr6_spectra.vot")

"""Select and save one spectrum from the results."""

# Select one spectrum and download data
k = 690  # Index of the spectrum to download
spectrum = spectra[k]
spectrum_url = spectrum.getdataurl()
print(f"Spectrum data URL: {spectrum_url}")
response = requests.get(spectrum_url)

# response is natively a VOTable file
with open('lamost_spectrum.vot', 'wb') as f:
    f.write(response.content)
print("\n-------------------------------------------------------------")
print("Spectrum saved to lamost_spectrum.vot")

"""Load the spectrum, plot it and save the figure."""

# Load the spectrum and plot it
spectrum_table = Table.read('lamost_spectrum.vot')
plt.plot(spectrum_table['spectral'], spectrum_table['flux'])
plt.xlabel('Wavelength (Angstrom)')
plt.ylabel('Flux (counts)')
plt.title('LAMOST DR6 MRS Spectrum')
plt.savefig('lamost_spectrum.pdf', dpi=300)
print("\n-------------------------------------------------------------")
print("Spectrum saved to lamost_spectrum.pdf")
plt.show()

"""# A sophisticated query using TAP

In this example we query tha Gaia DR3 dataset to download those records matching _several_ conditions, instead of filtering a large table after downloading it.

## Example 7: Tailored query of the Gaia DR3 archive.

| Key function/method  | Description |
|----------------------|-------------|
| ```pyvo.dal.TAPService``` | Class to instantiate a TAP service and include its methods and functions.

**Reference:** this example is adapted from the Mas-Buitrago's et al. (2024) tutorial in [zenodo](https://zenodo.org/records/14226117).

Import the necessary packages
"""

import pyvo
import matplotlib.pyplot as plt
from rich import print
from astropy import units as u
from astropy.coordinates import SkyCoord

"""Create the instance to the Gaia TAP service"""

tap_gaia = pyvo.dal.TAPService("https://gea.esac.esa.int/tap-server/tap")

"""Start with a generic query"""

# Define a query to the Gaia endpoint to get information about the available tables (schemas, types)
query = '''
SELECT DISTINCT schema_name, table_type
FROM tap_schema.tables
'''

# Execute the query and get the results as an Astropy Table
gaia_info = tap_gaia.search(query).to_table()
print(gaia_info)

"""Define a query to the selected Gaia endpoint to get information about tables in the 'gaiadr3' schema related to 'source'

"""

query = '''
SELECT table_name, schema_name, size, table_type
FROM tap_schema.tables
WHERE schema_name = 'gaiadr3'
AND table_name LIKE '%source%'
AND table_type = 'table'
'''

"""Execute the query and get the results as an Astropy Table"""

gaia_info = tap_gaia.search(query).to_table()
print(gaia_info)

"""Define and execute a query to get data from the gaiadr3.gaia_source table.

Coordinates selected correspond to the open cluster M67. Note the conditions applied to the query.
"""

# query for M67
OBJECT_NAME = 'M67'
SEARCH_RADIUS = 0.5*u.deg
coords = SkyCoord.from_name(OBJECT_NAME)
print(f"\nCoordinates of {OBJECT_NAME}: RA={coords.ra.deg} deg, Dec={coords.dec.deg} deg")
query = '''
SELECT TOP 5000 *
FROM gaiadr3.gaia_source
WHERE parallax_over_error > 10
AND ABS(pmra_error/pmra) < 0.1
AND ABS(pmdec_error/pmdec) < 0.1
AND pmra < -10.44
AND pmra > -11.44
AND pmdec < -2.39
AND pmdec > -3.39
AND 1=CONTAINS(
  POINT(ra,dec),
  CIRCLE(132.848,11.815,0.5))
'''

# Execute the query and get the results as an Astropy Table
gaia_data = tap_gaia.search(query).to_table()
print(gaia_data)

"""Plot a color-magnitude diagram and save data"""

# Plot a CMD using matplotlib
plt.figure(figsize=(8, 6))
plt.scatter(gaia_data['bp_rp'], gaia_data['phot_g_mean_mag'], s=1, color='blue')
plt.gca().invert_yaxis()  # Invert y-axis for magnitude
plt.xlabel('BP - RP [mag]')
plt.ylabel('G mean magnitude [mag]')
plt.title(f"CMD of {OBJECT_NAME} from Gaia DR3. Search radius: {SEARCH_RADIUS.to(u.arcmin)} ")
plt.savefig('cmd_gaia_m67.pdf', dpi=300)
plt.show()

# Save data to a VOTable file
gaia_data.write('gaia_m67.vot', format='votable', overwrite=True)
print("\n-------------------------------------")
print("Gaia data saved to gaia_m67.vot")