Adding DESILegacySurvey module

Author: volodymyrss

astroquery/desi/__init__.py

@@ -0,0 +1,41 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+"""
+DESI LegacySurvery
+
+https://www.legacysurvey.org/
+-------------------------
+
+:author: Gabriele Barni ([email protected])
+"""
+
+# Make the URL of the server, timeout and other items configurable
+# See <http://docs.astropy.org/en/latest/config/index.html#developer-usage>
+# 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`.
+    """
+    server = _config.ConfigItem(
+        ['https://portal.nersc.gov/cfs/cosmo/data/legacysurvey/',
+         ],
+        'base url')
+
+    timeout = _config.ConfigItem(
+        30,
+        'Time limit for connecting to template_module server.')
+
+
+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',
+           ]

astroquery/desi/core.py

@@ -0,0 +1,77 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+import urllib.error
+import requests
+
+import pyvo as vo
+import numpy as np
+
+from astroquery.exceptions import NoResultsWarning
+from astroquery.query import BaseQuery
+from astroquery.utils import commons, async_to_sync
+
+from . import conf
+
+
+__all__ = ['DESILegacySurvey', 'DESILegacySurveyClass']
+
+
+@async_to_sync
+class DESILegacySurveyClass(BaseQuery):
+    URL = conf.server
+    TIMEOUT = conf.timeout
+
+    def query_region(self, coordinates, radius, data_release=9):
+        """
+        Queries a region around the specified coordinates.
+
+        Parameters
+        ----------
+        coordinates : str or `astropy.coordinates`.
+            coordinates around which to query
+        radius : str or `astropy.units.Quantity`.
+            the radius of the cone search
+
+        Returns
+        -------
+        response : `astropy.table.Table`
+        """
+
+        url = 'https://datalab.noirlab.edu/tap'
+        tap_service = vo.dal.TAPService(url)
+        qstr = "SELECT all * FROM ls_dr" + str(data_release) + ".tractor WHERE dec>" + str(coordinates.dec.deg - radius.deg) + " and dec<" + str(
+            coordinates.dec.deg + radius.deg) + " and ra>" + str(coordinates.ra.deg - radius.deg / np.cos(coordinates.dec.deg * np.pi / 180.)) + " and ra<" + str(
+            coordinates.ra.deg + radius.deg / np.cos(coordinates.dec.deg * np.pi / 180))
+
+        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, data_release=9, pixels=None, radius=None, show_progress=True, image_band='g'):
+        """
+        Returns
+        -------
+        A list of `astropy.io.fits.HDUList` objects.
+        """
+
+        image_size_arcsec = radius.arcsec
+        pixsize = 2 * image_size_arcsec / pixels
+
+        image_url = 'https://www.legacysurvey.org/viewer/fits-cutout?ra=' + str(position.ra.deg) + '&dec=' + str(position.dec.deg) + '&size=' + str(
+            pixels) + '&layer=ls-dr' + str(data_release) + '&pixscale=' + str(pixsize) + '&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 e:
+            # TODO not sure this is the most suitable exception
+            raise NoResultsWarning(f"{str(e)} - Problem retrieving the file at the url: {str(image_url)}")
+
+        return [fits_file]
+
+
+DESILegacySurvey = DESILegacySurveyClass()

astroquery/desi/tests/__init__.py

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

astroquery/desi/tests/data/dummy_table.txt

@@ -0,0 +1,10 @@
+# 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'),
+             os.path.join('data', '*.fits.gz'),
+             ]
+    return {'astroquery.desi.tests': paths}

astroquery/desi/tests/data/dummy_tractor.fits

@@ -0,0 +1,132 @@
+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 pyvo.dal import TAPResults
+from urllib import parse
+from contextlib import contextmanager
+
+from ... import desi
+from ...utils.mocks import MockResponse
+from ...utils import commons
+
+DATA_FILES = {
+    'dummy_tap_table': 'dummy_table.txt',
+    'dummy_tractor_fits': 'dummy_tractor.fits',
+    'dummy_hdu_list_fits': 'hdu_list_images.fits'
+}
+
+coords = commons.ICRSCoord('11h04m27s +38d12m32s')
+radius = commons.ArcminRadiusGenerator(0.5)
+pixels = 60
+data_release = 9
+emispheres_list = ['north', 'south']
+
+
+@pytest.fixture
+def patch_get(request):
+    try:
+        mp = request.getfixturevalue("monkeypatch")
+    except AttributeError:  # pytest < 3
+        mp = request.getfuncargvalue("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
+        encoding = kwargs.get('encoding', None)
+        f = None
+        try:
+            if encoding == 'binary':
+                f = open(file_obj, 'rb')
+                yield f
+            else:
+                f = open(file_obj, 'rb')
+                yield f
+        finally:
+            if f is not None:
+                f.close()
+
+    try:
+        mp = request.getfixturevalue("monkeypatch")
+    except AttributeError:  # pytest < 3
+        mp = request.getfuncargvalue("monkeypatch")
+
+    mp.setattr(commons, 'get_readable_fileobj',
+               get_readable_fileobj_mockreturn)
+    return mp
+
+
+@pytest.fixture
+def patch_tap(request):
+    try:
+        mp = request.getfixturevalue("monkeypatch")
+    except AttributeError:  # pytest < 3
+        mp = request.getfuncargvalue("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:
+        content = open(filename, 'rb').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 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
+    data = fits.open(data_path(DATA_FILES[filetype]))
+    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, coords=coords, radius=radius):
+    result = desi.DESILegacySurvey.query_region(coords, radius)
+    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, radius=radius, pixels=pixels)
+
+    image_tester(images_list, 'dummy_hdu_list_fits')

astroquery/desi/tests/data/hdu_list_images.fits

@@ -0,0 +1,53 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+
+import pytest
+from astropy.io.fits import HDUList
+from astroquery.exceptions import NoResultsWarning
+
+
+@pytest.mark.remote_data
+class TestLegacySurveyClass:
+
+    def test_query_region(self):
+        import astroquery.desi
+        from astropy.coordinates import SkyCoord
+        from astropy.coordinates import Angle
+        from astropy.table import Table
+
+        ra = Angle('11h04m27s', unit='hourangle').degree
+        dec = Angle('+38d12m32s', unit='hourangle').degree
+        coordinates = SkyCoord(ra, dec, unit='degree')
+
+        radius = Angle(5, unit='arcmin')
+
+        query1 = astroquery.desi.DESILegacySurvey.query_region(coordinates, radius=radius, data_release=9)
+
+        assert isinstance(query1, Table)
+
+    @pytest.mark.parametrize("valid_inputs", [True, False])
+    def test_get_images(self, valid_inputs):
+        import astroquery.desi
+        from astropy.coordinates import SkyCoord
+        from astropy.coordinates import Angle
+
+        if valid_inputs:
+            ra = Angle('11h04m27s', unit='hourangle').degree
+            dec = Angle('+38d12m32s', unit='hourangle').degree
+            radius_input = 0.5  # arcmin
+            pixels = 60
+        else:
+            ra = Angle('86.633212', unit='degree').degree
+            dec = Angle('22.01446', unit='degree').degree
+            radius_input = 3  # arcmin
+            pixels = 1296000
+
+        pos = SkyCoord(ra, dec, unit='degree')
+        radius = Angle(radius_input, unit='arcmin')
+
+        if valid_inputs:
+            query1 = astroquery.desi.DESILegacySurvey.get_images(pos, data_release=9, radius=radius, pixels=pixels)
+            assert isinstance(query1, list)
+            assert isinstance(query1[0], HDUList)
+        else:
+            with pytest.raises(NoResultsWarning):
+                astroquery.desi.DESILegacySurvey.get_images(pos, data_release=9, radius=radius, pixels=pixels)

astroquery/desi/tests/setup_package.py

@@ -42,10 +42,13 @@ def FK4CoordGenerator(*args, **kwargs):
     return coord.SkyCoord(*args, frame='fk4', **kwargs)
 
 
+def ArcminRadiusGenerator(*args, **kwargs):
+    return coord.Angle(*args, unit='arcmin', **kwargs)
+
+
 ICRSCoord = coord.SkyCoord
 CoordClasses = (coord.SkyCoord, BaseCoordinateFrame)
 
-
 __all__ = ['send_request',
            'parse_coordinates',
            'TableList',