Add a medium parameter to WaveCoord.coord

CHANGELOG

@@ -3,6 +3,9 @@
 
 Development version, not released yet.
 
+- The `WaveCoord.coord` now accepts a medium to force getting the air or vacuum
+  wavelength coordinates.
+
 
 3.1 (03/04/2019)
 ----------------

lib/mpdaf/obj/coords.py

@@ -1640,7 +1640,7 @@ class WaveCoord:
                             atol=1E-2, rtol=0) and
                 self.wcs.wcs.ctype[0] == other.wcs.wcs.ctype[0])
 
-    def coord(self, pixel=None, unit=None):
+    def coord(self, pixel=None, unit=None, medium=None):
         """Return the coordinate corresponding to pixel.
 
         If pixel is None (default value), the full coordinate array is
@@ -1649,9 +1649,13 @@ class WaveCoord:
         Parameters
         ----------
         pixel : int, array or None.
-            pixel value.
+            Pixel value.
         unit : `astropy.units.Unit`
-            type of the wavelength coordinates
+            Unit of the wavelength coordinates
+        medium : str or None
+            Medium in which the wavelengths are returned: 'air' or 'vacuum'.
+            If None (default), the wavelength corresponding to the spectrum
+            CTYPE are returned.
 
         Returns
         -------
@@ -1669,7 +1673,25 @@ class WaveCoord:
         res = self.wcs.wcs_pix2world(pixelarr, 0)[0]
         if unit is not None:
             res = (res * self.unit).to(unit).value
-        return res[0] if np.isscalar(pixel) else res
+
+        result = res[0] if np.isscalar(pixel) else res
+
+        if medium is not None:
+            if medium not in ['air', 'vacuum']:
+                raise ValueError("Unknown 'medium' parameter value.")
+            ctype = self.wcs.wcs.ctype[0]  # Wave type
+            if ctype[:4] not in ['WAVE', 'AWAV']:
+                raise ValueError("No method to convert from %s to %s." %
+                                 (ctype, medium))
+
+            from .spectrum import airtovac, vactoair  # To avoid circ. import
+
+            if medium == "air" and ctype.startswith("WAVE"):
+                result = vactoair(result)
+            elif medium == "vacuum" and ctype.startswith("AWAV"):
+                result = airtovac(result)
+
+        return result
 
     def pixel(self, lbda, nearest=False, unit=None):
         """Return the decimal pixel corresponding to the wavelength lbda.

lib/mpdaf/obj/tests/test_coords.py

@@ -36,7 +36,9 @@ import numpy as np
 from astropy import wcs as pywcs
 from astropy.io import fits
 from mpdaf.obj import WCS, WaveCoord, deg2sexa, sexa2deg, determine_refframe
+from mpdaf.obj.spectrum import airtovac, vactoair
 from numpy.testing import assert_allclose, assert_array_equal
+import pytest
 
 from mpdaf.tests.utils import get_data_file
 
@@ -160,6 +162,29 @@ class TestWaveCoord:
         wave2 = wave.copy()
         assert wave.isEqual(wave2)
 
+    def test_coord(self):
+        """WaveCoord class: testing getting the coordinates"""
+        wave = WaveCoord(crval=0, cunit=u.nm, shape=10)
+        # By default the CTYPE is LINEAR and can't be converted to air or
+        # vacuum.
+        with pytest.raises(ValueError):
+            wave.coord(medium='air')
+        with pytest.raises(ValueError):
+            wave.coord(medium='vacuum')
+        wave.wcs.wcs.ctype = ['WAVE']
+        with pytest.raises(ValueError):
+            # Unknown parameter value
+            wave.coord(medium='vacuu')
+        np.testing.assert_array_equal(
+            wave.coord(medium='vacuum'), np.arange(10))
+        np.testing.assert_array_equal(
+            wave.coord(medium='air'), vactoair(np.arange(10)))
+        wave.wcs.wcs.ctype = ['AWAV']
+        np.testing.assert_array_equal(
+            wave.coord(medium='air'), np.arange(10))
+        np.testing.assert_array_equal(
+            wave.coord(medium='vacuum'), airtovac(np.arange(10)))
+
     def test_coord_transform(self):
         """WaveCoord class: testing coordinates transformations"""
         wave = WaveCoord(crval=0, cunit=u.nm, shape=10)