Merge 41ecabe into 5154874

.travis.yml

@@ -35,6 +35,10 @@ matrix:
           env: TOXENV='py36-test-casa'
           name: "Python 3.6 with minimal dependencies and CASA"
 
+        - python: 3.6
+          env: TOXENV='py36-test-casa-dev'
+          name: "Python 3.6, CASA, and dev versions of key dependencies"
+
         - python: 3.8
           env: TOXENV='py38-test-all-dev'
           name: "Python 3.8, all dependencies, and dev versions of key dependencies"

spectral_cube/base_class.py

@@ -577,14 +577,14 @@ def average_beams(self, threshold, mask='compute', warn=False):
             A new radio beam object that is the average of the unmasked beams
         """
         if mask == 'compute':
-            beam_mask = np.any(np.bitwise_and(self.mask.include(),
+            beam_mask = np.any(np.logical_and(self.mask.include(),
                                               self.goodbeams_mask[:,None,None]),
                                axis=(1,2))
         else:
             if mask.ndim > 1:
-                beam_mask = np.bitwise_and(mask, self.goodbeams_mask[:,None,None])
+                beam_mask = np.logical_and(mask, self.goodbeams_mask[:,None,None])
             else:
-                beam_mask = np.bitwise_and(mask, self.goodbeams_mask)
+                beam_mask = np.logical_and(mask, self.goodbeams_mask)
 
         # use private _beams here because the public one excludes the bad beams
         # by default

spectral_cube/io/casa_image.py

@@ -3,10 +3,12 @@
 import six
 import warnings
 import tempfile
+import shutil
 from astropy.io import fits
 from astropy.wcs import WCS
 from astropy import units as u
 from astropy.wcs.wcsapi.sliced_low_level_wcs import sanitize_slices
+from astropy import log
 import numpy as np
 from radio_beam import Beam, Beams
 
@@ -15,6 +17,7 @@
 from .. import SpectralCube, StokesSpectralCube, BooleanArrayMask, LazyMask, VaryingResolutionSpectralCube
 from .. import cube_utils
 from .. utils import BeamWarning, cached
+from .. import wcs_utils
 
 # Read and write from a CASA image. This has a few
 # complications. First, by default CASA does not return the
@@ -42,14 +45,20 @@ def wcs_casa2astropy(ia, coordsys):
     # to CASA by getting it to write out a FITS file and reading back in
     # using WCS
 
-    from casatasks import exportfits
-
     tmpimagefile = tempfile.mktemp() + '.image'
     tmpfitsfile = tempfile.mktemp() + '.fits'
     ia.newimagefromarray(outfile=tmpimagefile,
                          pixels=np.ones([1] * coordsys.naxes()),
                          csys=coordsys.torecord(), log=False)
-    exportfits(tmpimagefile, tmpfitsfile, stokeslast=False)
+    ia.done()
+    ia.open(tmpimagefile)
+    ia.tofits(tmpfitsfile, stokeslast=False)
+    ia.done()
+
+
+    # need to explicitly delete the tempfile to _force_ casa to close the
+    # handle
+    shutil.rmtree(tmpimagefile)
 
     return WCS(tmpfitsfile)
 
@@ -60,60 +69,112 @@ def __init__(self, filename, ia_kwargs={}):
 
         try:
             import casatools
-            self.ia = casatools.image()
+            self.iatool = casatools.image
+            tb = casatools.table()
         except ImportError:
             try:
-                from taskinit import iatool
-                self.ia = iatool()
+                from taskinit import iatool, tbtool
+                self.iatool = iatool
+                tb = tbtool()
             except ImportError:
                 raise ImportError("Could not import CASA (casac) and therefore cannot read CASA .image files")
 
-        # use the ia tool to get the file contents
-        try:
-            self.ia.open(filename)
-        except AssertionError as ex:
-            if 'must be of cReqPath type' in str(ex):
-                raise IOError("File {0} not found.  Error was: {1}"
-                              .format(filename, str(ex)))
-            else:
-                raise ex
 
         self.ia_kwargs = ia_kwargs
 
+        self.filename = filename
+
         self._cache = {}
 
-    @property
-    @cached
-    def shape(self):
-        return tuple(self.ia.shape()[::-1])
+        log.debug("Creating ArrayLikeCasa object")
+
+        # try to trick CASA into destroying the ia object
+        def getshape():
+            ia = self.iatool()
+            # use the ia tool to get the file contents
+            try:
+                ia.open(self.filename, cache=False)
+            except AssertionError as ex:
+                if 'must be of cReqPath type' in str(ex):
+                    raise IOError("File {0} not found.  Error was: {1}"
+                                  .format(self.filename, str(ex)))
+                else:
+                    raise ex
+
+            self.shape = tuple(ia.shape()[::-1])
+            self.dtype = np.dtype(ia.pixeltype())
+
+            ia.done()
+            ia.close()
+
+        getshape()
+
+        self.ndim = len(self.shape)
+
+        tb.open(self.filename)
+        dminfo = tb.getdminfo()
+        tb.done()
 
-    @property
-    @cached
-    def ndim(self):
-        return len(self.shape)
+        # unclear if this is always the correct callspec!!!
+        # (transpose requires this be backwards)
+        self.chunksize = dminfo['*1']['SPEC']['DEFAULTTILESHAPE'][::-1]
+
+
+        log.debug("Finished with initialization of ArrayLikeCasa object")
 
-    @property
-    @cached
-    def dtype(self):
-        return np.dtype(self.ia.pixeltype())
 
 
     def __getitem__(self, value):
 
+
+        log.debug(f"Retrieving slice {value} from {self}")
+
         value = sanitize_slices(value[::-1], self.ndim)
 
-        blc = [(slc.start or -1) if hasattr(slc, 'start') else slc for slc in value]
-        trc = [(slc.stop-1 if slc.stop is not None else -1)
+        # several cases:
+        # integer: just use an integer
+        # slice starting w/number: start number
+        # slice starting w/None: use -1
+        blc = [(-1 if slc.start is None else slc.start)
+               if hasattr(slc, 'start') else slc
+               for slc in value]
+        # slice ending w/number >= 1: end number -1 (CASA is end-inclusive)
+        # slice ending w/zero: use zero, not -1.
+        # slice ending w/negative: use it, but ?????
+        # slice ending w/None: use -1
+        trc = [((slc.stop-1 if slc.stop >= 1 else slc.stop)
+                if slc.stop is not None else -1)
                if hasattr(slc, 'stop') else slc for slc in value]
         inc = [(slc.step or 1) if hasattr(slc, 'step') else 1 for slc in value]
 
-        data = self.ia.getchunk(blc=blc, trc=trc, inc=inc, **self.ia_kwargs)
+
+        ia = self.iatool()
+        # use the ia tool to get the file contents
+        try:
+            ia.open(self.filename, cache=False)
+        except AssertionError as ex:
+            if 'must be of cReqPath type' in str(ex):
+                raise IOError("File {0} not found.  Error was: {1}"
+                              .format(self.filename, str(ex)))
+            else:
+                raise ex
+
+        log.debug(f'blc={blc}, trc={trc}, inc={inc}, kwargs={self.ia_kwargs}')
+        data = ia.getchunk(blc=blc, trc=trc, inc=inc, **self.ia_kwargs)
+        ia.done()
+        ia.close()
+
+        log.debug(f"Done retrieving slice {value} from {self}")
 
         # keep all sliced axes but drop all integer axes
         new_view = [slice(None) if isinstance(slc, slice) else 0
                     for slc in value]
 
-        return data[tuple(new_view)].transpose()
+        transposed_data = data[tuple(new_view)].transpose()
+
+        log.debug(f"Done transposing data with view {new_view}")
+
+        return transposed_data
 
 
 def load_casa_image(filename, skipdata=False,
@@ -127,17 +188,18 @@ def load_casa_image(filename, skipdata=False,
 
     try:
         import casatools
-        ia = casatools.image()
+        iatool = casatools.image
     except ImportError:
         try:
             from taskinit import iatool
-            ia = iatool()
         except ImportError:
             raise ImportError("Could not import CASA (casac) and therefore cannot read CASA .image files")
 
+    ia = iatool()
+
     # use the ia tool to get the file contents
     try:
-        ia.open(filename)
+        ia.open(filename, cache=False)
     except AssertionError as ex:
         if 'must be of cReqPath type' in str(ex):
             raise IOError("File {0} not found.  Error was: {1}"
@@ -147,25 +209,38 @@ def load_casa_image(filename, skipdata=False,
 
     # read in the data
     if not skipdata:
+        arrdata = ArraylikeCasaData(filename)
         # CASA data are apparently transposed.
-        data = dask.array.from_array(ArraylikeCasaData(filename))
+        data = dask.array.from_array(arrdata,
+                                     chunks=arrdata.chunksize,
+                                     name=filename
+                                    )
 
     # CASA stores validity of data as a mask
     if not skipvalid:
-        valid = dask.array.from_array(ArraylikeCasaData(filename,
-                                                        ia_kwargs={'getmask':
-                                                                   True}))
+        boolarr = ArraylikeCasaData(filename, ia_kwargs={'getmask': True})
+        valid = dask.array.from_array(boolarr, chunks=boolarr.chunksize,
+                                      name=filename+".mask"
+                                     )
 
     # transpose is dealt with within the cube object
 
     # read in coordinate system object
     casa_cs = ia.coordsys()
 
-    wcs = wcs_casa2astropy(ia, casa_cs)
-
     unit = ia.brightnessunit()
 
     beam_ = ia.restoringbeam()
+
+    ia.done()
+    ia.close()
+
+    wcs = wcs_casa2astropy(ia, casa_cs)
+
+    del casa_cs
+    del ia
+
+
     if 'major' in beam_:
         beam = Beam(major=u.Quantity(beam_['major']['value'], unit=beam_['major']['unit']),
                     minor=u.Quantity(beam_['minor']['value'], unit=beam_['minor']['unit']),
@@ -218,31 +293,36 @@ def load_casa_image(filename, skipdata=False,
         #    print new_order
         #    self.casa_cs.reorder(new_order)
 
-    # close the ia tool
-    ia.close()
 
     meta = {'filename': filename,
             'BUNIT': unit}
 
 
     if wcs.naxis == 3:
-        mask = BooleanArrayMask(valid, wcs)
+        data, wcs_slice = cube_utils._orient(data, wcs)
+        valid, _ = cube_utils._orient(valid, wcs)
+
+        mask = BooleanArrayMask(valid, wcs_slice)
         if 'beam' in locals():
-            cube = SpectralCube(data, wcs, mask, meta=meta, beam=beam)
+            cube = SpectralCube(data, wcs_slice, mask, meta=meta, beam=beam)
         elif 'beams' in locals():
-            cube = VaryingResolutionSpectralCube(data, wcs, mask, meta=meta, beams=beams)
+            cube = VaryingResolutionSpectralCube(data, wcs_slice, mask, meta=meta, beams=beams)
         else:
-            cube = SpectralCube(data, wcs, mask, meta=meta)
+            cube = SpectralCube(data, wcs_slice, mask, meta=meta)
+        # with #592, this is no longer true
         # we've already loaded the cube into memory because of CASA
         # limitations, so there's no reason to disallow operations
-        cube.allow_huge_operations = True
+        # cube.allow_huge_operations = True
+        assert cube.mask.shape == cube.shape
 
     elif wcs.naxis == 4:
+        valid, _ = cube_utils._split_stokes(valid, wcs)
         data, wcs = cube_utils._split_stokes(data, wcs)
         mask = {}
         for component in data:
             data_, wcs_slice = cube_utils._orient(data[component], wcs)
-            mask[component] = BooleanArrayMask(valid, wcs)
+            valid_, _ = cube_utils._orient(valid[component], wcs)
+            mask[component] = BooleanArrayMask(valid_, wcs_slice)
 
             if 'beam' in locals():
                 data[component] = SpectralCube(data_, wcs_slice, mask[component],
@@ -261,8 +341,11 @@ def load_casa_image(filename, skipdata=False,
 
 
         cube = StokesSpectralCube(stokes_data=data)
+        assert cube.I.mask.shape == cube.shape
+        assert wcs_utils.check_equality(cube.I.mask._wcs, cube.wcs)
     else:
         raise ValueError("CASA image has {0} dimensions, and therefore "
                          "is not readable by spectral-cube.".format(wcs.naxis))
 
+
     return cube

spectral_cube/io/casa_masks.py

@@ -3,7 +3,6 @@
 import numpy as np
 from astropy.io import fits
 import tempfile
-import warnings
 import os
 
 from ..wcs_utils import add_stokes_axis_to_wcs
@@ -88,6 +87,7 @@ def make_casa_mask(SpecCube, outname, append_to_image=True,
 
     cs = ia.coordsys()
 
+    ia.done()
     ia.close()
 
     temp2.close()
@@ -103,22 +103,26 @@ def make_casa_mask(SpecCube, outname, append_to_image=True,
     ia.newimagefromarray(outfile=maskpath,
                          pixels=mask_arr.astype('int16'),
                          overwrite=overwrite)
+    ia.done()
     ia.close()
 
-    ia.open(maskpath)
+    ia.open(maskpath, cache=False)
     ia.setcoordsys(cs.torecord())
 
+    ia.done()
     ia.close()
 
     if append_to_image:
         if img is None:
             raise TypeError("img argument must be specified to append the mask.")
 
-        ia.open(maskpath)
+        ia.open(maskpath, cache=False)
         ia.calcmask(maskname+">0.5")
+        ia.done()
         ia.close()
 
-        ia.open(img)
+        ia.open(img, cache=False)
         ia.maskhandler('copy', [maskpath+":mask0", maskname])
         ia.maskhandler('set', maskname)
+        ia.done()
         ia.close()

spectral_cube/io/fits.py

@@ -111,12 +111,8 @@ def read_data_fits(input, hdu=None, mode='denywrite', **kwargs):
 
     else:
 
-        hdulist = fits_open(input, mode=mode, **kwargs)
-
-        try:
+        with fits_open(input, mode=mode, **kwargs) as hdulist:
             return read_data_fits(hdulist, hdu=hdu)
-        finally:
-            hdulist.close()
 
     return array_hdu.data, array_hdu.header, beam_table