Add a spectral cutout example

Also include support for numpy arrays for band and channel, and autocorrect the band and channel values to be in ascending order

Author: jd-au

astroquery/casda/core.py

@@ -14,6 +14,7 @@
 from astropy.table import Table
 from astropy.io.votable import parse
 from astroquery import log
+import numpy as np
 
 # 3. local imports - use relative imports
 # commonly required local imports shown below as example
@@ -125,7 +126,7 @@ def _args_to_payload(self, radius=1*u.arcmin, **kwargs):
             if channel is not None:
                 raise ValueError("Either 'channel' or 'band' values may be provided but not both.")
 
-            if (not isinstance(band, (list, tuple))) or len(band) != 2 or \
+            if (not isinstance(band, (list, tuple, np.ndarray))) or len(band) != 2 or \
                     (band[0] is not None and not isinstance(band[0], u.Quantity)) or \
                     (band[1] is not None and not isinstance(band[1], u.Quantity)):
                 raise ValueError("The 'band' value must be a list of 2 wavelength or frequency values.")
@@ -137,22 +138,31 @@ def _args_to_payload(self, radius=1*u.arcmin, **kwargs):
             if bandBoundedLow or bandBoundedHigh:
                 unit = band[0].unit if bandBoundedLow else band[1].unit
                 if unit.physical_type == 'length':
-                    min_band = '-Inf' if not bandBoundedLow else str(band[0].to(u.m).value)
-                    max_band = '+Inf' if not bandBoundedHigh else str(band[1].to(u.m).value)
+                    min_band = '-Inf' if not bandBoundedLow else band[0].to(u.m).value
+                    max_band = '+Inf' if not bandBoundedHigh else band[1].to(u.m).value
                 elif unit.physical_type == 'frequency':
                     # Swap the order when changing frequency to wavelength
-                    min_band = '-Inf' if not bandBoundedHigh else str(band[1].to(u.m, equivalencies=u.spectral()).value)
-                    max_band = '+Inf' if not bandBoundedLow else str(band[0].to(u.m, equivalencies=u.spectral()).value)
+                    min_band = '-Inf' if not bandBoundedHigh else band[1].to(u.m, equivalencies=u.spectral()).value
+                    max_band = '+Inf' if not bandBoundedLow else band[0].to(u.m, equivalencies=u.spectral()).value
                 else:
                     raise ValueError("The 'band' values must be wavelengths or frequencies.")
+                # If values were provided in the wrong order, swap them
+                if bandBoundedLow and bandBoundedHigh and min_band > max_band:
+                    temp_val = min_band
+                    min_band = max_band
+                    max_band = temp_val
 
                 request_payload['BAND'] = f'{min_band} {max_band}'
 
         if channel is not None:
-            if not isinstance(channel, (list, tuple)) or len(channel) != 2 or \
-                    not isinstance(channel[0], int) or not isinstance(channel[1], int):
+            if not isinstance(channel, (list, tuple, np.ndarray)) or len(channel) != 2 or \
+                    not isinstance(channel[0], (int, np.integer)) or not isinstance(channel[1], (int, np.integer)):
                 raise ValueError("The 'channel' value must be a list of 2 integer values.")
-            request_payload['CHANNEL'] = f'{channel[0]} {channel[1]}'
+            if channel[0] <= channel[1]:
+                request_payload['CHANNEL'] = f'{channel[0]} {channel[1]}'
+            else:
+                # If values were provided in the wrong order, swap them
+                request_payload['CHANNEL'] = f'{channel[1]} {channel[0]}'
 
         return request_payload
 
@@ -310,7 +320,8 @@ def cutout(self, table, *, coordinates=None, radius=None, height=None,
 
         job_url = self._create_job(table, 'cutout_service', verbose)
 
-        cutout_spec = self._args_to_payload(radius=radius, **kwargs)
+        cutout_spec = self._args_to_payload(radius=radius, coordinates=coordinates, height=height, width=width, 
+               band=band, channel=channel, verbose=verbose)
 
         if not cutout_spec:
             raise ValueError("Please provide cutout parameters such as coordinates, band or channel.")

astroquery/casda/tests/test_casda.py

@@ -12,6 +12,7 @@
 from astropy.table import Table, Column
 from astropy.io.votable import parse
 from astroquery import log
+import numpy as np
 
 from astroquery.casda import Casda
 
@@ -348,6 +349,10 @@ def test_args_to_payload_band():
     assert payload['BAND'] == '0.195 0.215'
     assert list(payload.keys()) == ['BAND']
 
+    payload = casda._args_to_payload(band=(0.215*u.m, 0.195*u.m))
+    assert payload['BAND'] == '0.195 0.215'
+    assert list(payload.keys()) == ['BAND']
+
     payload = casda._args_to_payload(band=(0.195*u.m, 21.5*u.cm))
     assert payload['BAND'] == '0.195 0.215'
     assert list(payload.keys()) == ['BAND']
@@ -364,6 +369,10 @@ def test_args_to_payload_band():
     assert payload['BAND'] == '0.19986163866666667 0.21112144929577467'
     assert list(payload.keys()) == ['BAND']
 
+    payload = casda._args_to_payload(band=np.array([1.5, 1.42])*u.GHz)
+    assert payload['BAND'] == '0.19986163866666667 0.21112144929577467'
+    assert list(payload.keys()) == ['BAND']
+
     payload = casda._args_to_payload(band=(None, 1.5*u.GHz))
     assert payload['BAND'] == '0.19986163866666667 +Inf'
     assert list(payload.keys()) == ['BAND']
@@ -393,7 +402,7 @@ def test_args_to_payload_band_invalid():
     assert "The 'band' values must have the same kind of units." in str(excinfo.value)
 
     with pytest.raises(ValueError) as excinfo:
-        casda._args_to_payload(band=(1.42*u.radian, 21*u.deg))
+        casda._args_to_payload(band=[1.42*u.radian, 21*u.deg])
     assert "The 'band' values must be wavelengths or frequencies." in str(excinfo.value)
 
     with pytest.raises(ValueError) as excinfo:
@@ -407,7 +416,11 @@ def test_args_to_payload_channel():
     assert payload['CHANNEL'] == '0 30'
     assert list(payload.keys()) == ['CHANNEL']
 
-    payload = casda._args_to_payload(channel=(17, 23))
+    payload = casda._args_to_payload(channel=np.array([17, 23]))
+    assert payload['CHANNEL'] == '17 23'
+    assert list(payload.keys()) == ['CHANNEL']
+
+    payload = casda._args_to_payload(channel=(23, 17))
     assert payload['CHANNEL'] == '17 23'
     assert list(payload.keys()) == ['CHANNEL']
 

docs/casda/casda.rst

@@ -131,7 +131,7 @@ Spatial and spectral parameters can be combined to produce sub-cubes.
 
 Once completed, the cutouts can be downloaded as described in the section above.
 
-An example script to download a cutout from the Rapid ASKAP Continuum Survey (RACS) at a specified position is shown
+An example script to download a 2D cutout from the Rapid ASKAP Continuum Survey (RACS) at a specified position is shown
 below:
 
 .. doctest-skip::
@@ -152,6 +152,26 @@ below:
     >>> url_list = casda.cutout(subset[:1], coordinates=centre, radius=14*u.arcmin)
     >>> filelist = casda.download_files(url_list, savedir='/tmp')
 
+An example script to download a 3D cutout from the WALLABY Pre-Pilot Eridanus cube at a specified position and velocity
+is shown below:
+
+.. doctest-skip::
+
+    >>> from astropy import coordinates, units as u, wcs
+    >>> from astroquery.casda import Casda
+    >>> import getpass
+    >>> centre = coordinates.SkyCoord.from_name('NGC 1371')
+    >>> username = '[email protected]'
+    >>> password = getpass.getpass(str("Enter your OPAL password: "))
+    >>> casda = Casda(username, password)
+    >>> result = Casda.query_region(centre, radius=30*u.arcmin)
+    >>> public_data = Casda.filter_out_unreleased(result)
+    >>> eridanus_cube = public_data[public_data['filename'] == 'Eridanus_full_image_V3.fits']
+    >>> vel = np.array([1250, 1600])*u.km/u.s
+    >>> freq = vel.to(u.Hz, equivalencies=u.doppler_radio(1.420405751786*u.GHz))
+    >>> url_list = casda.cutout(eridanus_cube, coordinates=centre, radius=9*u.arcmin, band=freq)
+    >>> filelist = casda.download_files(url_list, savedir='/tmp')
+
 
 Reference/API
 =============