Enhanced results feature in ALMA

CHANGES.rst

@@ -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
 ^^^
 

astroquery/alma/__init__.py

@@ -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'
            ]

astroquery/alma/core.py

@@ -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
 
@@ -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.*']
 
@@ -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
 
@@ -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.
 
@@ -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
 
@@ -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
 
@@ -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')

astroquery/alma/tests/data/alma-shapes.xml

@@ -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>

astroquery/alma/tests/test_alma.py

@@ -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
 
 
@@ -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'),

astroquery/alma/tests/test_alma_remote.py

@@ -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
@@ -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'}