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..99e3aea501 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,106 @@ 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:] + points = SkyCoord( + [(float(tokens[ii]), float(tokens[ii + 1])) * u.deg + for ii in range(0, len(tokens), 2)], frame=csys) + 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 = s_region[s_region.index('(') + 1:s_region.rindex( + ')')].strip() + # omit the first char in the string to force it look for the second + # occurrence + last_pos = None + not_operation = False # not operation - signals that the next substring is just the not operation + not_shape = False # not shape - signals that the next substring is a not shape + for shape in re.finditer(r'not|polygon|circle', input_regions): + pos = shape.span()[0] + if last_pos is None: + last_pos = pos + continue # this is the first elem + if shape.group() == 'not': + not_operation = True + elif not_operation: + not_shape = True + not_operation = False + last_pos = pos + continue + if res is not None: + next_shape = _get_region( + input_regions[last_pos:pos - 1].strip(' ()').split()) + if not_shape: + res = (res or next_shape) ^ next_shape + not_shape = False + else: + 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 +477,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. @@ -390,6 +492,9 @@ def query_object_async(self, object_name, *, public=True, science : bool True to return only science datasets, False to return only calibration, None to return both + enhanced_results : bool + True to return a table with quantities instead of just values. It + also returns the footprint as `regions` objects. payload : dict Dictionary of additional keywords. See `help`. """ @@ -398,10 +503,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 @@ -419,6 +526,9 @@ def query_region_async(self, coordinate, radius, *, public=True, calibration, None to return both payload : dict Dictionary of additional keywords. See `help`. + enhanced_results : bool + True to return a table with quantities instead of just values. It + also returns the footprints as `regions` objects. """ rad = radius if not isinstance(radius, u.Quantity): @@ -433,11 +543,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 @@ -458,6 +569,9 @@ def query_async(self, payload, *, public=True, science=True, Flag to indicate whether to simply return the payload. maxrec : integer Cap on the amount of records returned. Default is no limit. + enhanced_results : bool + True to return a table with quantities instead of just values. It + also returns the footprints as `regions` objects. Returns ------- @@ -492,7 +606,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..b5d40b63cd --- /dev/null +++ b/astroquery/alma/tests/data/alma-shapes.xml @@ -0,0 +1,65 @@ + + + + + + + + + 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)
ADS/JAO.ALMA#2016.1.00298.S30.57219940476191412.388698412698416106.99696191911804Union ICRS ( Polygon 24.688037 10.311983 24.686396 10.310080 24.684212 10.308752 24.681715 10.308155 24.679150 10.308346 24.675920 10.309890 24.673492 10.312999 24.669332 10.312311 24.665417 10.313483 24.663126 10.315516 24.661753 10.318233 24.661454 10.320747 24.661960 10.323226 24.663552 10.325814 24.666441 10.327836 24.666411 10.331617 24.667968 10.334789 24.671222 10.337213 24.675213 10.337795 24.677255 10.340016 24.680013 10.341357 24.683594 10.341533 24.686482 10.340470 24.688448 10.338832 24.689759 10.336779 24.692298 10.336273 24.694935 10.334697 24.697356 10.330856 24.697547 10.327333 24.696046 10.323872 24.697052 10.320175 24.696122 10.316262 24.692952 10.313013 24.690552 10.312113 Not (Polygon 24.681021 10.324923 24.681816 10.324257 24.682229 10.324855 24.681434 10.325521) Polygon 36.472346 14.440666 36.470405 14.438548 36.468128 14.437316 36.465567 14.436828 36.462472 14.437285 36.459379 14.439158 36.457576 14.441702 36.454774 14.442332 36.451320 14.444138 36.449514 14.445972 36.448380 14.448255 36.448210 14.452274 36.449149 14.454632 36.450564 14.456378 36.450365 14.459666 36.451674 14.462953 36.454306 14.465382 36.456726 14.466327 36.459320 14.466501 36.461261 14.468619 36.464035 14.470010 36.467148 14.470312 36.470623 14.469241 36.472651 14.467639 36.474091 14.465464 36.476894 14.464834 36.479112 14.463490 36.481645 14.459695 36.481908 14.456175 36.480365 14.452572 36.481385 14.449020 36.480517 14.445090 36.477360 14.441784 36.474940 14.440839 Not (Polygon 36.465340 14.453623 36.465668 14.452977 36.466327 14.453544 36.465998 14.454189))
+ + + diff --git a/astroquery/alma/tests/test_alma.py b/astroquery/alma/tests/test_alma.py index 3c81b18b77..c6c1177ac8 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,102 @@ 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) == 5 + enhanced_result = get_enhanced_table(result) + assert len(enhanced_result) == 5 + # 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) + elif row['obs_publisher_did'] == 'ADS/JAO.ALMA#2016.1.00298.S': + # pick random points inside and outside + inside = SkyCoord('1:38:44 10:18:55', unit=(u.hourangle, u.deg)) + hole = SkyCoord('1:38:44 10:19:31.5', unit=(u.hourangle, u.deg)) + pix_inside = regions.PixCoord.from_sky(inside, ww) + pix_outside = regions.PixCoord.from_sky(hole, ww) + assert pix_region.contains(pix_inside) + # assert not pix_region.contains(pix_outside) + 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..db8362a57e 100644 --- a/astroquery/alma/tests/test_alma_remote.py +++ b/astroquery/alma/tests/test_alma_remote.py @@ -14,7 +14,7 @@ from pyvo.dal.exceptions import DALOverflowWarning from astroquery.exceptions import CorruptDataWarning -from .. import Alma +from astroquery.alma import Alma, get_enhanced_table # ALMA tests involving staging take too long, leading to travis timeouts # TODO: make this a configuration item @@ -62,14 +62,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 = 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..5dd9a86e91 100644 --- a/docs/alma/alma.rst +++ b/docs/alma/alma.rst @@ -211,7 +211,7 @@ One can also query by keyword, spatial resolution, etc: ... "AND science_observation='T'") # doctest: +IGNORE_OUTPUT -Use the ''help_tap'' method to learn about the ALMA 'ObsCore' keywords and +Use the ``help_tap`` method to learn about the ALMA 'ObsCore' keywords and their types. .. doctest-remote-data:: @@ -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 `~astroquery.alma.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 require 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