diff --git a/astroquery/desi/__init__.py b/astroquery/desi/__init__.py new file mode 100644 index 0000000000..40466fbc12 --- /dev/null +++ b/astroquery/desi/__init__.py @@ -0,0 +1,42 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +""" +DESI LegacySurvery + +https://www.legacysurvey.org/ +----------------------------- + +:author: Gabriele Barni (Gabriele.Barni@unige.ch) +""" + +# Make the URL of the server, timeout and other items configurable +# See +# for docs and examples on how to do this +# Below is a common use case +from astropy import config as _config + + +class Conf(_config.ConfigNamespace): + """ + Configuration parameters for `astroquery.desi`. + """ + legacysurvey_service_url = _config.ConfigItem( + ['https://www.legacysurvey.org/viewer/fits-cutout', + ], + 'url for the LegacySurvey service') + + tap_service_url = _config.ConfigItem( + ['https://datalab.noirlab.edu/tap', + ], + 'url for the TAP service') + + +conf = Conf() + +# Now import your public class +# Should probably have the same name as your module +from .core import DESILegacySurvey, DESILegacySurveyClass + +__all__ = ['DESILegacySurvey', 'DESILegacySurveyClass', + 'Conf', 'conf', + ] diff --git a/astroquery/desi/core.py b/astroquery/desi/core.py new file mode 100644 index 0000000000..5c45e90962 --- /dev/null +++ b/astroquery/desi/core.py @@ -0,0 +1,103 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +import urllib.error +import requests +import warnings + +import pyvo as vo +import numpy as np +import astropy.coordinates as coord + +from astropy import units as u +from astroquery.exceptions import NoResultsWarning +from astroquery.query import BaseQuery +from astroquery.utils import commons +from astroquery.desi import conf + +__all__ = ['DESILegacySurvey', 'DESILegacySurveyClass'] + + +class DESILegacySurveyClass(BaseQuery): + + def query_region(self, coordinates, *, width=0.5*u.arcmin, data_release=9): + """ + Queries a region around the specified coordinates. + + Parameters + ---------- + coordinates : `~astropy.coordinates.SkyCoord` + coordinates around which to query. + width : `~astropy.units.Quantity`, optional + the width of the region. If missing, set to default + value of 0.5 arcmin. + data_release : int + the data release of the LegacySurvey to use. + + Returns + ------- + response : `~astropy.table.Table` + """ + + tap_service = vo.dal.TAPService(conf.tap_service_url) + coordinates_transformed = coordinates.transform_to(coord.ICRS) + + qstr = (f"SELECT all * FROM ls_dr{data_release}.tractor WHERE " + f"dec<{(coordinates_transformed.dec + width).to_value(u.deg)} and " + f"ra>{coordinates_transformed.ra.to_value(u.deg) - width.to_value(u.deg) / np.cos(coordinates_transformed.dec)} and " + f"ra<{coordinates_transformed.ra.to_value(u.deg) + width.to_value(u.deg) / np.cos(coordinates_transformed.dec)}") + + tap_result = tap_service.run_sync(qstr) + tap_result = tap_result.to_table() + # filter out duplicated lines from the table + mask = tap_result['type'] != 'D' + filtered_table = tap_result[mask] + + return filtered_table + + def get_images(self, position, *, pixels=None, width=0.5*u.arcmin, data_release=9, show_progress=True, image_band='g'): + """ + Downloads the images for a certain region of interest. + + Parameters + ---------- + position : `astropy.coordinates`. + coordinates around which we define our region of interest. + width : `~astropy.units.Quantity`, optional + the width of our region of interest. + data_release : int, optional + the data release of the LegacySurvey to use. + show_progress : bool, optional + Whether to display a progress bar if the file is downloaded + from a remote server. Default is True. + image_band : str, optional + Default to 'g' + + Returns + ------- + list : A list of `~astropy.io.fits.HDUList` objects. + """ + + position_transformed = position.transform_to(coord.ICRS) + + image_size_arcsec = width.arcsec + pixsize = 2 * image_size_arcsec / pixels + + image_url = (f"{conf.legacysurvey_service_url}?" + f"ra={position_transformed.ra.deg}&" + f"dec={position_transformed.dec.deg}&" + f"size={pixels}&" + f"layer=ls-dr{data_release}&" + f"pixscale={pixsize}&" + f"bands={image_band}") + + file_container = commons.FileContainer(image_url, encoding='binary', show_progress=show_progress) + + try: + fits_file = file_container.get_fits() + except (requests.exceptions.HTTPError, urllib.error.HTTPError) as exp: + fits_file = None + warnings.warn(f"{str(exp)} - Problem retrieving the file at the url: {image_url}", NoResultsWarning) + + return [fits_file] + + +DESILegacySurvey = DESILegacySurveyClass() diff --git a/astroquery/desi/tests/__init__.py b/astroquery/desi/tests/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/astroquery/desi/tests/data/dummy_table.txt b/astroquery/desi/tests/data/dummy_table.txt new file mode 100644 index 0000000000..da5be1746a --- /dev/null +++ b/astroquery/desi/tests/data/dummy_table.txt @@ -0,0 +1,8 @@ +ra,dec,objid,type +166.10552527002142,38.20797162140221,877,PSF +166.10328347620825,38.211862682863625,855,PSF +166.1146193911762,38.20826292586543,991,PSF +166.1138080401007,38.20883307659884,3705,DUP +166.11382707824612,38.20885008952696,982,SER +166.11779248975387,38.211159276708706,1030,PSF +166.11865123008005,38.2147201669633,1039,PSF diff --git a/astroquery/desi/tests/data/dummy_tractor.fits b/astroquery/desi/tests/data/dummy_tractor.fits new file mode 100644 index 0000000000..254b4cb3ba Binary files /dev/null and b/astroquery/desi/tests/data/dummy_tractor.fits differ diff --git a/astroquery/desi/tests/data/hdu_list_images.fits b/astroquery/desi/tests/data/hdu_list_images.fits new file mode 100644 index 0000000000..8b77018ffb Binary files /dev/null and b/astroquery/desi/tests/data/hdu_list_images.fits differ diff --git a/astroquery/desi/tests/setup_package.py b/astroquery/desi/tests/setup_package.py new file mode 100644 index 0000000000..2e89d237d6 --- /dev/null +++ b/astroquery/desi/tests/setup_package.py @@ -0,0 +1,9 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +import os + + +def get_package_data(): + paths = [os.path.join('data', '*.txt'), + os.path.join('data', '*.fits'), + ] + return {'astroquery.desi.tests': paths} diff --git a/astroquery/desi/tests/test_desi.py b/astroquery/desi/tests/test_desi.py new file mode 100644 index 0000000000..0d9ae737ff --- /dev/null +++ b/astroquery/desi/tests/test_desi.py @@ -0,0 +1,115 @@ +import astropy.io.votable +import pytest +import os +import numpy as np +import pyvo as vo + +from astropy.table import Table +from astropy.io import fits +from astropy import coordinates as coord +from astroquery.utils.mocks import MockResponse +from astroquery.utils import commons +from astroquery import desi +from urllib import parse +from contextlib import contextmanager + + +DATA_FILES = { + 'dummy_tap_table': 'dummy_table.txt', + 'dummy_tractor_fits': 'dummy_tractor.fits', + 'dummy_hdu_list_fits': 'hdu_list_images.fits' +} + +coords = coord.SkyCoord('11h04m27s +38d12m32s', frame='icrs') +width = coord.Angle(0.5, unit='arcmin') +pixels = 60 +data_release = 9 +emispheres_list = ['north', 'south'] + + +@pytest.fixture +def patch_get(request): + mp = request.getfixturevalue("monkeypatch") + + mp.setattr(desi.DESILegacySurvey, '_request', get_mockreturn) + return mp + + +@pytest.fixture +def patch_get_readable_fileobj(request): + @contextmanager + def get_readable_fileobj_mockreturn(filename, **kwargs): + file_obj = data_path(DATA_FILES['dummy_hdu_list_fits']) # TODO: add images option + # f = None + with open(file_obj, 'rb') as f: + yield f + + mp = request.getfixturevalue("monkeypatch") + + mp.setattr(commons, 'get_readable_fileobj', + get_readable_fileobj_mockreturn) + return mp + + +@pytest.fixture +def patch_tap(request): + mp = request.getfixturevalue("monkeypatch") + + mp.setattr(vo.dal.TAPService, 'run_sync', tap_mockreturn) + return mp + + +def get_mockreturn(method, url, params=None, timeout=10, **kwargs): + parsed_url = parse.urlparse(url) + splitted_parsed_url = parsed_url.path.split('/') + url_filename = splitted_parsed_url[-1] + filename = None + content = None + if url_filename.startswith('tractor-'): + filename = data_path(DATA_FILES['dummy_tractor_fits']) + + if filename is not None: + with open(filename, 'rb') as f: + content = f.read() + return MockResponse(content) + + +def tap_mockreturn(url, params=None, timeout=10, **kwargs): + content_table = Table.read(data_path(DATA_FILES['dummy_tap_table']), + format='ascii.csv', comment='#') + votable_table = astropy.io.votable.from_table(content_table) + return vo.dal.TAPResults(votable_table) + + +def data_path(filename): + data_dir = os.path.join(os.path.dirname(__file__), 'data') + return os.path.join(data_dir, filename) + + +def compare_result_data(result, data): + for col in result.colnames: + if result[col].dtype.type is np.string_ or result[col].dtype.type is np.str_: + assert np.array_equal(result[col], data[col]) + else: + np.testing.assert_allclose(result[col], data[col]) + + +def image_tester(images, filetype): + assert type(images) == list + with fits.open(data_path(DATA_FILES[filetype])) as data: + assert images[0][0].header == data[0].header + assert np.array_equal(images[0][0].data, data[0].data) + + +def test_coords_query_region(patch_tap): + result = desi.DESILegacySurvey.query_region(coords, width=width) + data = Table.read(data_path(DATA_FILES['dummy_tap_table']), + format='ascii.csv', comment='#') + data['objid'] = data['objid'].astype(np.int64) + compare_result_data(result, data) + + +def test_coords_get_images(patch_get_readable_fileobj, dr=data_release): + images_list = desi.DESILegacySurvey.get_images(coords, data_release=dr, width=width, pixels=pixels) + + image_tester(images_list, 'dummy_hdu_list_fits') diff --git a/astroquery/desi/tests/test_desi_remote.py b/astroquery/desi/tests/test_desi_remote.py new file mode 100644 index 0000000000..be607955cd --- /dev/null +++ b/astroquery/desi/tests/test_desi_remote.py @@ -0,0 +1,41 @@ +import pytest + +from astroquery.desi import DESILegacySurvey +from astropy.io.fits import HDUList +from astropy.coordinates import SkyCoord +from astropy.coordinates import Angle +from astropy.table import Table +from astroquery.exceptions import NoResultsWarning + + +@pytest.mark.remote_data +class TestLegacySurveyClass: + + def test_query_region(self): + ra = 166.1125 + dec = 38.209 + coordinates = SkyCoord(ra, dec, unit='degree') + + width = Angle(15, unit='arcsec') + + query1 = DESILegacySurvey.query_region(coordinates, width=width, data_release=9) + + assert isinstance(query1, Table) + + @pytest.mark.parametrize(('ra', 'dec', 'width', 'pixels'), + ((166.1125, 38.209, 0.5, 60),)) + def test_get_images(self, ra, dec, width, pixels): + pos = SkyCoord(ra, dec, unit='degree') + width = Angle(width, unit='arcmin') + + query1 = DESILegacySurvey.get_images(pos, pixels=pixels, width=width, data_release=9) + assert isinstance(query1, list) + assert isinstance(query1[0], HDUList) + + def test_noresults_warning(self): + # Using position with no coverage + pos = SkyCoord(86.633212, 22.01446, unit='degree') + width = Angle(3, unit='arcmin') + + with pytest.warns(NoResultsWarning): + DESILegacySurvey.get_images(pos, width=width, pixels=100) diff --git a/astroquery/utils/commons.py b/astroquery/utils/commons.py index 3d149fa43a..0f99b3cf87 100644 --- a/astroquery/utils/commons.py +++ b/astroquery/utils/commons.py @@ -24,10 +24,8 @@ from ..exceptions import TimeoutError, InputWarning from .. import version - CoordClasses = (SkyCoord, BaseCoordinateFrame) - __all__ = ['send_request', 'parse_coordinates', 'TableList', diff --git a/docs/desi/desi.rst b/docs/desi/desi.rst new file mode 100644 index 0000000000..eb66967644 --- /dev/null +++ b/docs/desi/desi.rst @@ -0,0 +1,99 @@ +.. _astroquery.desi: + +********************************************* +DESI LegacySurvey Queries (`astroquery.desi`) +********************************************* + +Getting started +=============== + +This module provides a way to query the DesiLegacySurvey service. +Presented below are examples that illustrate the different types of queries +that can be formulated. + +Query a region +============== + +This example shows how to query a certain region with DesiLegacySurvey. +We'll use a set of coordinates that define the region of interest, +and search within a 15 arcsec box. + +.. doctest-remote-data:: + + >>> from astroquery.desi import DESILegacySurvey + >>> from astropy.coordinates import Angle, SkyCoord + >>> ra = Angle('11h04m27s', unit='hourangle').degree + >>> dec = Angle('+38d12m32s', unit='hourangle').degree + >>> coordinates = SkyCoord(ra, dec, unit='degree') + >>> width = Angle(15, unit='arcsec') + >>> table_out = DESILegacySurvey.query_region(coordinates, width=width, data_release=9) + >>> print(table_out[:5]) # doctest: +IGNORE_OUTPUT + ls_id dec ra ... type wise_coadd_id + ... + ---------------- ----------------- ------------------ ... ---- ------------- + 9907734382838387 38.12628495570797 166.0838654387131 ... PSF 1667p378 + 9907734382838488 38.12798336424771 166.0922968862182 ... PSF 1667p378 + 9907734382838554 38.12858283671958 166.0980673954384 ... REX 1667p378 + 9907734382838564 38.12840803351445 166.09921863682337 ... EXP 1667p378 + 9907734382838584 38.12836885301038 166.10070750146636 ... REX 1667p378 + +The result is an astropy.Table. + +Get images +========== + +To download images for a certain region of interest, +we can define our region in the same way used in the example above. + +.. doctest-remote-data:: + + >>> from astroquery.desi import DESILegacySurvey + >>> from astropy.coordinates import Angle, SkyCoord + >>> + >>> ra = Angle('11h04m27s', unit='hourangle').degree + >>> dec = Angle('+38d12m32s', unit='hourangle').degree + >>> radius_input = 0.5 + >>> pixels = 60 + >>> + >>> pos = SkyCoord(ra, dec, unit='degree') + >>> width = Angle(radius_input, unit='arcmin') + >>> im = DESILegacySurvey.get_images(pos, pixels=pixels, width=width, data_release=9) + +All the information we need can be found within the object "im". + +.. doctest-remote-data:: + + >>> hdul = im[0] + >>> hdul[0].header + SIMPLE = T / file does conform to FITS standard + BITPIX = -32 / number of bits per data pixel + NAXIS = 2 / number of data axes + NAXIS1 = 60 / length of data axis 1 + NAXIS2 = 60 / length of data axis 2 + EXTEND = T / FITS dataset may contain extensions + COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy + COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H + BANDS = 'g ' + BAND0 = 'g ' + CTYPE1 = 'RA---TAN' / TANgent plane + CTYPE2 = 'DEC--TAN' / TANgent plane + CRVAL1 = 166.1125 / Reference RA + CRVAL2 = 38.2088888888889 / Reference Dec + CRPIX1 = 30.5 / Reference x + CRPIX2 = 30.5 / Reference y + CD1_1 = -0.000277777777777778 / CD matrix + CD1_2 = 0. / CD matrix + CD2_1 = 0. / CD matrix + CD2_2 = 0.000277777777777778 / CD matrix + IMAGEW = 60. / Image width + IMAGEH = 60. / Image height + +The variable "im" is a list of `~astropy.io.fits.HDUList` objects, one entry for +each corresponding object. + + +Reference/API +============= + +.. automodapi:: astroquery.desi + :no-inheritance-diagram: diff --git a/docs/index.rst b/docs/index.rst index 8e5ab7e85a..5c7e75bc21 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -185,6 +185,7 @@ The following modules have been completed using a common API: cds/cds.rst linelists/cdms/cdms.rst dace/dace.rst + desi/desi.rst esa/hsa/hsa.rst esa/hubble/hubble.rst esa/iso/iso.rst