ccdproc fails with error "operand units do not match"
[Archive.zip](https://github.com/astropy/ccdproc/files/321651/Archive.zip)

Hi ccdproc gurus,

I work on Mac with latest version of OS and Python 2.7 ressources linked to ccdproc. I wrote a short script to deflat a few images, and I got:

Traceback (most recent call last):
  File "deflat.py", line 34, in <module>
    data_exposure = img_expos)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/ccdproc/log_meta.py", line 68, in wrapper
    result = func(_args, *_kwd)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/ccdproc/core.py", line 208, in ccd_process
    nccd = nccd.subtract(master_bias)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/ccdproc/ccddata.py", line 444, in subtract
    result = super(CCDData, self).subtract(other)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/astropy/nddata/mixins/ndarithmetic.py", line 215, in subtract
    operand, propagate_uncertainties, "subtraction", np.subtract)
  File "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/astropy/nddata/mixins/ndarithmetic.py", line 87, in _arithmetic
    raise ValueError("operand units do not match")
ValueError: operand units do not match

After digging into the code, I strongly suspect the value of the gain (2.35 electrons/adu) not be to properly applied. In line 349 of ccddata.py, the unit to be used (result_unit) is 'electron'. Then, in function subtract, line 433, other and self have different units. I attach here my script (dir_wrk to be modified) and relevant images.

diff --git a/CHANGES.rst b/CHANGES.rst
index 8c49004..ebbfeff 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -10,6 +10,9 @@ Other Changes and Additions
 
 - Auto-identify files with extension ``fts`` as FITS files. [#355, #364]
 
+- Raise more explicit exception if unit of uncalibrated image and master do
+  not match in ``subtract_bias`` or ``subtract_dark``. [#361, #366]
+
 Bug Fixes
 ^^^^^^^^^
 
diff --git a/ccdproc/core.py b/ccdproc/core.py
index 34fac91..8391a0b 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -13,6 +13,7 @@ from astropy.modeling import fitting
 from astropy import stats
 from astropy.nddata import StdDevUncertainty
 from astropy.wcs.utils import proj_plane_pixel_area
+import astropy  # To get the version.
 
 from scipy import ndimage
 
@@ -90,15 +91,21 @@ def ccd_process(ccd, oscan=None, trim=None, error=False, master_bias=None,
         Default is ``False``.
 
     master_bias : `~ccdproc.CCDData` or None, optional
-        A master bias frame to be subtracted from ccd.
+        A master bias frame to be subtracted from ccd. The unit of the
+        master bias frame should match the unit of the image **after
+        gain correction**.
         Default is ``None``.
 
     dark_frame : `~ccdproc.CCDData` or None, optional
-        A dark frame to be subtracted from the ccd.
+        A dark frame to be subtracted from the ccd. The unit of the
+        master dark frame should match the unit of the image **after
+        gain correction**.
         Default is ``None``.
 
     master_flat : `~ccdproc.CCDData` or None, optional
-        A master flat frame to be divided into ccd.
+        A master flat frame to be divided into ccd. The unit of the
+        master flat frame should match the unit of the image **after
+        gain correction**.
         Default is ``None``.
 
     bad_pixel_mask : `numpy.ndarray` or None, optional
@@ -541,7 +548,17 @@ def subtract_bias(ccd, master):
     result : `~ccdproc.CCDData`
         CCDData object with bias subtracted.
     """
-    result = ccd.subtract(master)
+
+    try:
+        result = ccd.subtract(master)
+    except ValueError as e:
+        if 'operand units' in str(e):
+            raise u.UnitsError("Unit '{}' of the uncalibrated image does not "
+                               "match unit '{}' of the calibration "
+                               "image".format(ccd.unit, master.unit))
+        else:
+            raise e
+
     result.meta = ccd.meta.copy()
     return result
 
@@ -623,14 +640,29 @@ def subtract_dark(ccd, master, dark_exposure=None, data_exposure=None,
             raise TypeError("exposure times must be astropy.units.Quantity "
                             "objects.")
 
-    if scale:
-        master_scaled = master.copy()
-        # data_exposure and dark_exposure are both quantities,
-        # so we can just have subtract do the scaling
-        master_scaled = master_scaled.multiply(data_exposure / dark_exposure)
-        result = ccd.subtract(master_scaled)
-    else:
-        result = ccd.subtract(master)
+    try:
+        if scale:
+            master_scaled = master.copy()
+            # data_exposure and dark_exposure are both quantities,
+            # so we can just have subtract do the scaling
+            master_scaled = master_scaled.multiply(data_exposure /
+                                                   dark_exposure)
+            result = ccd.subtract(master_scaled)
+        else:
+            result = ccd.subtract(master)
+    except (u.UnitsError, u.UnitConversionError, ValueError) as e:
+        # Astropy LTS (v1) returns a ValueError, not a UnitsError, so catch
+        # that if it appears to really be a UnitsError.
+        if (isinstance(e, ValueError) and
+                'operand units' not in str(e) and
+                astropy.__version__.startswith('1.0')):
+            raise e
+
+        # Make the error message a little more explicit than what is returned
+        # by default.
+        raise u.UnitsError("Unit '{}' of the uncalibrated image does not "
+                           "match unit '{}' of the calibration "
+                           "image".format(ccd.unit, master.unit))
 
     result.meta = ccd.meta.copy()
     return result
diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 3cdb2f9..0235234 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -10,6 +10,7 @@ import astropy.units as u
 from astropy.wcs import WCS
 
 from astropy.nddata import StdDevUncertainty
+import astropy
 
 from numpy.testing import assert_array_equal
 from astropy.tests.helper import pytest
@@ -299,7 +300,7 @@ def test_subtract_bias_fails(ccd_data):
         subtract_bias(ccd_data, bias)
     # should fail because units don't match
     bias = CCDData(np.zeros_like(ccd_data), unit=u.meter)
-    with pytest.raises(ValueError):
+    with pytest.raises(u.UnitsError):
         subtract_bias(ccd_data, bias)
 
 
@@ -353,28 +354,68 @@ def test_subtract_dark_fails(ccd_data):
     # can be anything.
     ccd_data.header['exptime'] = 30.0
     master = ccd_data.copy()
+
     # Do we fail if we give one of dark_exposure, data_exposure but not both?
     with pytest.raises(TypeError):
         subtract_dark(ccd_data, master, dark_exposure=30 * u.second)
     with pytest.raises(TypeError):
         subtract_dark(ccd_data, master, data_exposure=30 * u.second)
+
     # Do we fail if we supply dark_exposure and data_exposure and exposure_time
     with pytest.raises(TypeError):
         subtract_dark(ccd_data, master, dark_exposure=10 * u.second,
                       data_exposure=10 * u.second,
                       exposure_time='exptime')
+
     # Fail if we supply none of the exposure-related arguments?
     with pytest.raises(TypeError):
         subtract_dark(ccd_data, master)
+
     # Fail if we supply exposure time but not a unit?
     with pytest.raises(TypeError):
         subtract_dark(ccd_data, master, exposure_time='exptime')
+
     # Fail if ccd_data or master are not CCDData objects?
     with pytest.raises(TypeError):
         subtract_dark(ccd_data.data, master, exposure_time='exptime')
     with pytest.raises(TypeError):
         subtract_dark(ccd_data, master.data, exposure_time='exptime')
 
+    # Fail if units do not match...
+
+    # ...when there is no scaling?
+    master = CCDData(ccd_data)
+    master.unit = u.meter
+
+    with pytest.raises(u.UnitsError) as e:
+        subtract_dark(ccd_data, master, exposure_time='exptime',
+                      exposure_unit=u.second)
+    assert "uncalibrated image" in str(e.value)
+
+
+def test_unit_mismatch_behaves_as_expected(ccd_data):
+    """
+    Test to alert us to any changes in how errors are raised in astropy when units
+    do not match.
+    """
+    bad_unit = ccd_data.copy()
+    bad_unit.unit = u.meter
+
+    if astropy.__version__.startswith('1.0'):
+        expected_error = ValueError
+        expected_message = 'operand units'
+    else:
+        expected_error = u.UnitConversionError
+        # Make this an empty string, which always matches. In this case
+        # we are really only checking by the type of error raised.
+        expected_message = ''
+
+    # Did we raise the right error?
+    with pytest.raises(expected_error) as e:
+        ccd_data.subtract(bad_unit)
+
+    # Was the error message as expected?
+    assert expected_message in str(e)
 
 # test for flat correction
 @pytest.mark.data_scale(10)
diff --git a/docs/ccdproc/reduction_toolbox.rst b/docs/ccdproc/reduction_toolbox.rst
index f3a70df..5029ccd 100644
--- a/docs/ccdproc/reduction_toolbox.rst
+++ b/docs/ccdproc/reduction_toolbox.rst
@@ -185,7 +185,7 @@ the overscan axis ends up being the *second* axis, which is numbered 1 in
 python zero-based numbering.
 
 With the arguments in this example the overscan is averaged over the overscan
-columns (i.e. 2000 through 2031) and then subtracted row-by-row from the
+columns (i.e. 200 through 231) and then subtracted row-by-row from the
 image. The ``median`` argument can be used to median combine instead.
 
 This example is not very realistic: typically one wants to fit a low-order
@@ -205,8 +205,8 @@ Trim an image
 +++++++++++++
 
 The overscan-subtracted image constructed above still contains the overscan
-portion. We are assuming came from a FITS file in which ``NAXIS1=2032`` and
-``NAXIS2=1000``, in which the last 32 columns along ``NAXIS1`` are overscan.
+portion. We are assuming came from a FITS file in which ``NAXIS1=232`` and
+``NAXIS2=100``, in which the last 32 columns along ``NAXIS1`` are overscan.
 
 Trim it using `~ccdproc.trim_image`,shown below in both python-
 style and FITS-style indexing:
@@ -299,15 +299,25 @@ Basic Processing
 ----------------
 
 All of the basic processing steps can be accomplished in a single step using
-`~ccdproc.ccd_process`.   This step will call overscan correct, trim, gain
-correct, add a bad pixel mask, create an uncertainty frame, subtract the master
-bias, and flat-field the image.  These can be run together as:
+`~ccdproc.ccd_process`. This step will call overscan correct, trim, gain
+correct, add a bad pixel mask, create an uncertainty frame, subtract the
+master bias, and flat-field the image. The unit of the master calibration
+frames must match that of the image *after* the gain, if any, is applied. In
+the example below, ``img`` has unit ``adu``, but the master frames have unit
+``electron``. These can be run together as:
 
      >>> ccd = ccdproc.CCDData(img, unit=u.adu)
-     >>> nccd = ccdproc.ccd_process(ccd, oscan='[1:10,1:100]', 
-     ...                            trim='[10:100, 1:100]',
-     ...                            error=True, gain=2.0*u.electron/u.adu,
-     ...                            readnoise = 5*u.electron) 
+     >>> ccd.header['exposure'] = 30.0  # for dark subtraction
+     >>> nccd = ccdproc.ccd_process(ccd, oscan='[201:232,1:100]',
+     ...                            trim='[1:200, 1:100]',
+     ...                            error=True,
+     ...                            gain=2.0*u.electron/u.adu,
+     ...                            readnoise=5*u.electron,
+     ...                            dark_frame=master_dark,
+     ...                            exposure_key='exposure',
+     ...                            exposure_unit=u.second,
+     ...                            dark_scale=True,
+     ...                            master_flat=master_flat)
 
 
 Reprojecting onto a different image footprint
['ccdproc/core.py', 'docs/ccdproc/reduction_toolbox.rst', 'ccdproc/tests/test_ccdproc.py', 'CHANGES.rst']