Skip to content

Commit

Permalink
Enhanced results feature in ALMA
Browse files Browse the repository at this point in the history
  • Loading branch information
Adrian Damian authored and Adrian Damian committed May 30, 2024
1 parent 25e6521 commit d841ca5
Show file tree
Hide file tree
Showing 8 changed files with 334 additions and 13 deletions.
5 changes: 5 additions & 0 deletions CHANGES.rst
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,11 @@ New Tools and Services
Service fixes and enhancements
------------------------------

alma
^^^^

- Added method to return quantities instead of values and regions footprint in alma [#2855]

mpc
^^^

Expand Down
4 changes: 2 additions & 2 deletions astroquery/alma/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,8 @@ class Conf(_config.ConfigNamespace):

conf = Conf()

from .core import Alma, AlmaClass, ALMA_BANDS
from .core import Alma, AlmaClass, ALMA_BANDS, get_enhanced_table

__all__ = ['Alma', 'AlmaClass',
'Conf', 'conf', 'ALMA_BANDS'
'Conf', 'conf', 'ALMA_BANDS', 'get_enhanced_table'
]
108 changes: 101 additions & 7 deletions astroquery/alma/core.py
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
from astropy.utils.console import ProgressBar
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import SkyCoord

from pyvo.dal.sia2 import SIA2_PARAMETERS_DESC, SIA2Service

Expand All @@ -32,7 +33,7 @@
from . import conf, auth_urls
from astroquery.exceptions import CorruptDataWarning

__all__ = {'AlmaClass', 'ALMA_BANDS'}
__all__ = {'AlmaClass', 'ALMA_BANDS', 'get_enhanced_table'}

__doctest_skip__ = ['AlmaClass.*']

Expand Down Expand Up @@ -207,6 +208,92 @@ def _gen_sql(payload):
return sql + where


def get_enhanced_table(result):
"""
Returns an enhanced table with quantities instead of values and regions for footprints.
Note that this function is dependent on the ``astropy`` - affiliated ``regions`` package.
"""
try:
import regions
except ImportError:
print(
"Could not import astropy-regions, which is a requirement for get_enhanced_table function in alma"
"Please refer to http://astropy-regions.readthedocs.io/en/latest/installation.html for how to install it.")

def _parse_stcs_string(input):
csys = 'icrs'

def _get_region(tokens):
if tokens[0] == 'polygon':
if csys == tokens[1].lower():
tokens = tokens[2:]
else:
tokens = tokens[1:]
try:
points = SkyCoord(
[(float(tokens[ii]), float(tokens[ii + 1])) * u.deg
for ii in range(0, len(tokens), 2)], frame=csys)
except Exception as ex:
print(ex)
return regions.PolygonSkyRegion(points)
elif tokens[0] == 'circle':
if csys == tokens[1].lower():
tokens = tokens[2:]
else:
tokens = tokens[1:]
return regions.CircleSkyRegion(
SkyCoord(float(tokens[0]), float(tokens[1]), unit=u.deg),
float(tokens[2]) * u.deg)
else:
raise ValueError("Unrecognized shape: " + tokens[0])
s_region = input.lower().strip()
if s_region.startswith('union'):
res = None
# skip the union operator
input_regions = re.split(r'\(|\)', s_region)[1]
# omit the first char in the string to force it look for the second
# occurrence
last_pos = None
for shape in re.finditer(r'polygon|circle', input_regions):
pos = shape.span()[0]
if last_pos is None:
last_pos = pos
continue # this is the first elem
if res is not None:
next_shape = _get_region(input_regions[last_pos:pos-1].strip().split())
res = res | next_shape
else:
res = _get_region(input_regions[last_pos:pos-1].strip().split())
last_pos = pos
if last_pos:
next_shape = _get_region(input_regions[last_pos:].strip().split())
res = res | next_shape
return res
elif "not" in s_region:
# shape with "holes"
comps = s_region.split('not')
result = _get_region(comps[0].strip(' ()').split())
for comp in comps[1:]:
hole = _get_region(comp.strip(' ()').split())
result = (result or hole) ^ hole
return result
else:
return _get_region(s_region.split())
prep_table = result.to_qtable()
s_region_parser = None
for field in result.resultstable.fields:
if ('s_region' == field.ID) and \
('obscore:Char.SpatialAxis.Coverage.Support.Area' == field.utype):
if 'adql:REGION' == field.xtype:
s_region_parser = _parse_stcs_string
# this is where to add other xtype parsers such as shape
break
if (s_region_parser):
for row in prep_table:
row['s_region'] = s_region_parser(row['s_region'])
return prep_table


class AlmaAuth(BaseVOQuery, BaseQuery):
"""Authentication session information for passing credentials to an OIDC instance
Expand Down Expand Up @@ -376,7 +463,8 @@ def tap_url(self):
return self._tap_url

def query_object_async(self, object_name, *, public=True,
science=True, payload=None, **kwargs):
science=True, payload=None, enhanced_results=False,
**kwargs):
"""
Query the archive for a source name.
Expand All @@ -398,10 +486,12 @@ def query_object_async(self, object_name, *, public=True,
else:
payload = {'source_name_resolver': object_name}
return self.query_async(public=public, science=science,
payload=payload, **kwargs)
payload=payload, enhanced_results=enhanced_results,
**kwargs)

def query_region_async(self, coordinate, radius, *, public=True,
science=True, payload=None, **kwargs):
science=True, payload=None, enhanced_results=False,
**kwargs):
"""
Query the ALMA archive with a source name and radius
Expand Down Expand Up @@ -433,11 +523,12 @@ def query_region_async(self, coordinate, radius, *, public=True,
payload['ra_dec'] = ra_dec

return self.query_async(public=public, science=science,
payload=payload, **kwargs)
payload=payload, enhanced_results=enhanced_results,
**kwargs)

def query_async(self, payload, *, public=True, science=True,
legacy_columns=False, get_query_payload=False,
maxrec=None, **kwargs):
maxrec=None, enhanced_results=False, **kwargs):
"""
Perform a generic query with user-specified payload
Expand Down Expand Up @@ -492,7 +583,10 @@ def query_async(self, payload, *, public=True, science=True,
result = self.query_tap(query, maxrec=maxrec)

if result is not None:
result = result.to_table()
if enhanced_results:
result = get_enhanced_table(result)
else:
result = result.to_table()
else:
# Should not happen
raise RuntimeError('BUG: Unexpected result None')
Expand Down
58 changes: 58 additions & 0 deletions astroquery/alma/tests/data/alma-shapes.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,58 @@
<?xml version="1.0" encoding="UTF-8"?>
<VOTABLE xmlns="http://www.ivoa.net/xml/VOTable/v1.3" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.4">
<RESOURCE type="results">
<INFO name="QUERY_STATUS" value="OK" />
<INFO name="QUERY_TIMESTAMP" value="2023-11-21T22:54:11.627" />
<INFO name="QUERY" value="select obs_publisher_did, s_ra, s_dec, frequency, s_region from ivoa.obscore where obs_id='uid://A001/X1465/X2f66.source.HOPS_56.spw.33' or obs_id='uid://A001/X1284/X146e.source.J2229-6910.spw.17' or obs_id='uid://A001/X145/X3d6.source.CMZ_Far-Side.spw.29' or obs_id='uid://A001/X158f/X745.source.CANDELS_-_GOODS_South.spw.29'" />
<TABLE>
<FIELD name="obs_publisher_did" datatype="char" arraysize="33*" ucd="meta.ref.ivoid" utype="obscore:Curation.PublisherDID">
<DESCRIPTION>publisher dataset identifier</DESCRIPTION>
</FIELD>
<FIELD name="s_ra" datatype="double" ucd="pos.eq.ra" unit="deg" utype="obscore:Char.SpatialAxis.Coverage.Location.Coord.Position2D.Value2.C1" xtype="adql:DOUBLE">
<DESCRIPTION>RA of central coordinates</DESCRIPTION>
</FIELD>
<FIELD name="s_dec" datatype="double" ucd="pos.eq.dec" unit="deg" utype="obscore:Char.SpatialAxis.Coverage.Location.Coord.Position2D.Value2.C2" xtype="adql:DOUBLE">
<DESCRIPTION>DEC of central coordinates</DESCRIPTION>
</FIELD>
<FIELD name="frequency" datatype="double" ucd="em.freq;obs;meta.main" unit="GHz" utype="obscore:Char.SpectralAxis.Coverage.Location.refval" xtype="adql:DOUBLE">
<DESCRIPTION>Observed (tuned) reference frequency on the sky.</DESCRIPTION>
</FIELD>
<FIELD name="s_region" datatype="char" arraysize="*" ucd="pos.outline;obs.field" utype="obscore:Char.SpatialAxis.Coverage.Support.Area" xtype="adql:REGION">
<DESCRIPTION>region bounded by observation</DESCRIPTION>
</FIELD>
<DATA>
<TABLEDATA>
<TR>
<TD>ADS/JAO.ALMA#2017.1.00358.S</TD>
<TD>337.25073645418405</TD>
<TD>-69.17507615555645</TD>
<TD>103.4188615997524</TD>
<TD>Circle ICRS 337.250736 -69.175076 0.008223</TD>
</TR>
<TR>
<TD>ADS/JAO.ALMA#2013.1.01014.S</TD>
<TD>266.3375206596382</TD>
<TD>-29.04661351679843</TD>
<TD>102.53500890383512</TD>
<TD>Union ICRS ( Polygon 266.149398 -29.290586 266.136226 -29.288373 266.119371 -29.256173 266.139594 -29.227875 266.179424 -29.227495 266.200095 -29.255356 266.200602 -29.262616 266.180920 -29.290650 Polygon 266.201178 -29.180809 266.185467 -29.157212 266.203140 -29.128221 266.210443 -29.124729 266.248282 -29.126396 266.267989 -29.156792 266.250347 -29.186328 266.243653 -29.189773 266.209701 -29.189725 Polygon 266.277718 -29.102223 266.255292 -29.101217 266.236836 -29.068844 266.257370 -29.039510 266.295961 -29.038762 266.318687 -29.073177 266.299686 -29.101619 Polygon 266.675227 -28.729991 266.662115 -28.727829 266.645189 -28.695695 266.665158 -28.667319 266.704771 -28.666784 266.725471 -28.694565 266.726013 -28.701823 266.706579 -28.729933)</TD>
</TR>
<TR>
<TD>ADS/JAO.ALMA#2019.1.00458.S</TD>
<TD>83.83095833330503</TD>
<TD>-5.260619444445853</TD>
<TD>218.0000053819096</TD>
<TD>Polygon ICRS 83.831185 -5.264207 83.827356 -5.260845 83.830732 -5.257032 83.834561 -5.260394</TD>
</TR>
<TR>
<TD>ADS/JAO.ALMA#2021.1.00547.S</TD>
<TD>53.11983108973377</TD>
<TD>-27.807152863261976</TD>
<TD>104.63751498393975</TD>
<TD>Polygon ICRS 53.095155 -27.862094 53.079461 -27.868350 53.070554 -27.864641 53.071634 -27.854716 53.063795 -27.848756 53.065662 -27.840191 53.057830 -27.834233 53.059697 -27.825669 53.051864 -27.819711 53.053728 -27.811149 53.045898 -27.805188 53.047766 -27.796624 53.040165 -27.791131 53.040727 -27.783321 53.168803 -27.740304 53.176556 -27.741168 53.180746 -27.747006 53.178683 -27.754479 53.186636 -27.761019 53.183259 -27.768949 53.192133 -27.774067 53.189238 -27.783465 53.198118 -27.788584 53.196264 -27.797148 53.204102 -27.803100 53.202244 -27.811668 53.209857 -27.817148 53.209307 -27.824959 Not (Polygon 53.160470 -27.801462 53.170993 -27.798574 53.173032 -27.789459 53.183944 -27.785173 53.175851 -27.781416 53.177962 -27.770657 53.168833 -27.764876 53.148836 -27.770820 53.146803 -27.779935 53.136929 -27.783386 53.144761 -27.789340 53.143267 -27.798785 53.153398 -27.795816)</TD>
</TR>
</TABLEDATA>
</DATA>
</TABLE>
<INFO name="placeholder" value="ignore" />
</RESOURCE>
</VOTABLE>
93 changes: 92 additions & 1 deletion astroquery/alma/tests/test_alma.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,15 @@

from astropy import units as u
from astropy import coordinates as coord
from astropy import wcs
from astropy.table import Table
from astropy.io import votable
from astropy.coordinates import SkyCoord
from astropy.time import Time
import pyvo

from astroquery.alma import Alma
from astroquery.alma.core import _gen_sql, _OBSCORE_TO_ALMARESULT
from astroquery.alma.core import _gen_sql, _OBSCORE_TO_ALMARESULT, get_enhanced_table
from astroquery.alma.tapsql import _val_parse


Expand Down Expand Up @@ -354,6 +357,94 @@ def test_query():
)


@pytest.mark.filterwarnings("ignore::astropy.utils.exceptions.AstropyUserWarning")
def test_enhanced_table():
pytest.importorskip('regions')
import regions # to silence checkstyle
data = votable.parse(os.path.join(DATA_DIR, 'alma-shapes.xml'))
result = pyvo.dal.DALResults(data)
assert len(result) == 4
enhanced_result = get_enhanced_table(result)
assert len(enhanced_result) == 4
# generic ALMA WCS
ww = wcs.WCS(naxis=2)
ww.wcs.crpix = [250.0, 250.0]
ww.wcs.cdelt = [-7.500000005754e-05, 7.500000005754e-05]
ww.wcs.ctype = ['RA---SIN', 'DEC--SIN']
for row in enhanced_result:
# check other quantities
assert row['s_ra'].unit == u.deg
assert row['s_dec'].unit == u.deg
assert row['frequency'].unit == u.GHz
ww.wcs.crval = [row['s_ra'].value, row['s_dec'].value]
sky_center = SkyCoord(row['s_ra'].value, row['s_dec'].value, unit=u.deg)
pix_center = regions.PixCoord.from_sky(sky_center, ww)
s_region = row['s_region']
pix_region = s_region.to_pixel(ww)
if isinstance(s_region, regions.CircleSkyRegion):
# circle: https://almascience.org/aq/?mous=uid:%2F%2FA001%2FX1284%2FX146e
assert s_region.center.name == 'icrs'
assert s_region.center.ra.value == 337.250736
assert s_region.center.ra.unit == u.deg
assert s_region.center.dec.value == -69.175076
assert s_region.center.dec.unit == u.deg
assert s_region.radius.unit == u.deg
assert s_region.radius.value == 0.008223
x_min = pix_region.center.x - pix_region.radius
x_max = pix_region.center.x + pix_region.radius
y_min = pix_region.center.y - pix_region.radius
y_max = pix_region.center.y + pix_region.radius
assert pix_region.contains(pix_center)
elif isinstance(s_region, regions.PolygonSkyRegion):
# simple polygon: https://almascience.org/aq/?mous=uid:%2F%2FA001%2FX1296%2FX193
assert s_region.vertices.name == 'icrs'
x_min = pix_region.vertices.x.min()
x_max = pix_region.vertices.x.max()
y_min = pix_region.vertices.y.min()
y_max = pix_region.vertices.y.max()
assert pix_region.contains(pix_center)
elif isinstance(s_region, regions.CompoundSkyRegion):
x_min = pix_region.bounding_box.ixmin
x_max = pix_region.bounding_box.ixmax
y_min = pix_region.bounding_box.iymin
y_max = pix_region.bounding_box.iymax
if row['obs_publisher_did'] == 'ADS/JAO.ALMA#2013.1.01014.S':
# Union type of footprint: https://almascience.org/aq/?mous=uid:%2F%2FA001%2FX145%2FX3d6
# image center is outside
assert not pix_region.contains(pix_center)
# arbitrary list of points inside each of the 4 polygons
inside_pts = ['17:46:43.655 -28:41:43.956',
'17:45:06.173 -29:04:01.549',
'17:44:53.675 -29:09:19.382',
'17:44:38.584 -29:15:31.836']
for inside in [SkyCoord(coords, unit=(u.hourangle, u.deg)) for coords in inside_pts]:
pix_inside = regions.PixCoord.from_sky(inside, ww)
assert pix_region.contains(pix_inside)
else:
# polygon with "hole": https://almascience.org/aq/?mous=uid:%2F%2FA001%2FX158f%2FX745
assert pix_region.contains(pix_center)
# this is an arbitrary point in the footprint "hole"
outside_point = SkyCoord('03:32:38.689 -27:47:32.750',
unit=(u.hourangle, u.deg))
pix_outside = regions.PixCoord.from_sky(outside_point, ww)
assert not pix_region.contains(pix_outside)
else:
assert False, "Unsupported shape"
assert x_min <= x_max
assert y_min <= y_max

# example of how to plot the footprints
# artist = pix_region.as_artist()
# import matplotlib.pyplot as plt
# axes = plt.subplot(projection=ww)
# axes.set_aspect('equal')
# axes.add_artist(artist)
# axes.set_xlim([x_min, x_max])
# axes.set_ylim([y_min, y_max])
# pix_region.plot()
# plt.show()


def test_sia():
sia_mock = Mock()
empty_result = Table.read(os.path.join(DATA_DIR, 'alma-empty.txt'),
Expand Down
20 changes: 17 additions & 3 deletions astroquery/alma/tests/test_alma_remote.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
from pyvo.dal.exceptions import DALOverflowWarning

from astroquery.exceptions import CorruptDataWarning
import astroquery
from .. import Alma

# ALMA tests involving staging take too long, leading to travis timeouts
Expand Down Expand Up @@ -62,14 +63,27 @@ def test_public(self, alma):
for row in results:
assert row['data_rights'] == 'Proprietary'

@pytest.mark.filterwarnings(
"ignore::astropy.utils.exceptions.AstropyUserWarning")
def test_s_region(self, alma):
pytest.importorskip('regions')
import regions # to silence checkstyle
alma.help_tap()
result = alma.query_tap("select top 3 s_region from ivoa.obscore")
enhanced_result = astroquery.alma.core.get_enhanced_table(result)
for row in enhanced_result:
assert isinstance(row['s_region'], (regions.CircleSkyRegion,
regions.PolygonSkyRegion,
regions.CompoundSkyRegion))

@pytest.mark.filterwarnings(
"ignore::astropy.utils.exceptions.AstropyUserWarning")
def test_SgrAstar(self, tmp_path, alma):
alma.cache_location = tmp_path

result_s = alma.query_object('Sgr A*', legacy_columns=True)
result_s = alma.query_object('Sgr A*', legacy_columns=True, enhanced_results=True)

assert '2013.1.00857.S' in result_s['Project code']
# "The Brick", g0.253, is in this one
# assert b'2011.0.00217.S' in result_c['Project code'] # missing cycle 1 data

def test_freq(self, alma):
payload = {'frequency': '85..86'}
Expand Down
Loading

0 comments on commit d841ca5

Please sign in to comment.