diff --git a/CHANGES.rst b/CHANGES.rst index 191c945556..8d6480219a 100644 --- a/CHANGES.rst +++ b/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 ^^^ diff --git a/astroquery/alma/__init__.py b/astroquery/alma/__init__.py index 552a073a8c..e30235f2a4 100644 --- a/astroquery/alma/__init__.py +++ b/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' ] diff --git a/astroquery/alma/core.py b/astroquery/alma/core.py index a1a7e4b47d..a52c962580 100644 --- a/astroquery/alma/core.py +++ b/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') diff --git a/astroquery/alma/tests/data/alma-shapes.xml b/astroquery/alma/tests/data/alma-shapes.xml new file mode 100644 index 0000000000..5b2e2071da --- /dev/null +++ b/astroquery/alma/tests/data/alma-shapes.xml @@ -0,0 +1,58 @@ + + + + + + + + + publisher dataset identifier + + + RA of central coordinates + + + DEC of central coordinates + + + Observed (tuned) reference frequency on the sky. + + + region bounded by observation + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
ADS/JAO.ALMA#2017.1.00358.S337.25073645418405-69.17507615555645103.4188615997524Circle ICRS 337.250736 -69.175076 0.008223
ADS/JAO.ALMA#2013.1.01014.S266.3375206596382-29.04661351679843102.53500890383512Union 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)
ADS/JAO.ALMA#2019.1.00458.S83.83095833330503-5.260619444445853218.0000053819096Polygon ICRS 83.831185 -5.264207 83.827356 -5.260845 83.830732 -5.257032 83.834561 -5.260394
ADS/JAO.ALMA#2021.1.00547.S53.11983108973377-27.807152863261976104.63751498393975Polygon 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)
+ + + diff --git a/astroquery/alma/tests/test_alma.py b/astroquery/alma/tests/test_alma.py index 3c81b18b77..d15086dc7b 100644 --- a/astroquery/alma/tests/test_alma.py +++ b/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'), diff --git a/astroquery/alma/tests/test_alma_remote.py b/astroquery/alma/tests/test_alma_remote.py index 3db8969899..ce07ba6964 100644 --- a/astroquery/alma/tests/test_alma_remote.py +++ b/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'} diff --git a/docs/alma/alma.rst b/docs/alma/alma.rst index ea97bb777f..6e5ac5793d 100644 --- a/docs/alma/alma.rst +++ b/docs/alma/alma.rst @@ -294,6 +294,65 @@ their types. velocity_resolution double m/s Estimated velocity resolution from all the spectral windows, from frequency resolution. +Query Results +============= + +Results of queries are returned in tabular format. For convenience, +the ''get_enhanced_table'' function can be used to have the initial result +in a more useful format, i.e turn values into quantities, footprint into +shape, etc. (Note: this requires the ''regions'' Python package to be installed. + + +.. doctest-remote-data:: + + >>> from astroquery.alma import Alma, get_enhanced_table + >>> alma = Alma() + >>> alma.archive_url = 'https://almascience.eso.org' # optional to make doctest work + >>> res = alma.query_tap("select top 1 * from ivoa.ObsCore where obs_publisher_did='ADS/JAO.ALMA#2011.0.00087.S'") + >>> enhanced_res = get_enhanced_table(res) + >>> enhanced_res[0]['s_ra'] + + >>> enhanced_res[0]['s_region'] + )> + +To further draw the footprint: + +.. doctest-skip:: + + >>> from astropy import wcs + >>> import matplotlib.pyplot as plt + >>> # Create a WCS; for plotting, all that matters is that it is centered on our target region + >>> 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'] + >>> ww.wcs.crval = [enhanced_res[0]['s_ra'].value, enhanced_res[0]['s_dec'].value] + >>> pix_region = enhanced_res[0]['s_region'].to_pixel(ww) + >>> artist = pix_region.as_artist() + >>> axes = plt.subplot(projection=ww) + >>> axes.set_aspect('equal') + >>> axes.add_artist(artist) + >>> axes.axis(pix_region.bounding_box.extent) + >>> pix_region.plot() + >>> plt.show() + + +.. image:: footprint.png + :align: center + :scale: 75% + :alt: observation footprint + + +The above footprint could be transformed into a pixel region and have the mask +extracted or combined with other regions. Refer to the Astropy affiliated +''regions'' package for more details. + + Downloading Data ================ diff --git a/docs/alma/footprint.png b/docs/alma/footprint.png new file mode 100644 index 0000000000..d4066c4e5f Binary files /dev/null and b/docs/alma/footprint.png differ