radio-astro-tools / spectral-cube

Library for reading and analyzing astrophysical spectral data cubes
http://spectral-cube.rtfd.org
BSD 3-Clause "New" or "Revised" License
98 stars 65 forks source link

Masks for cubes that are created directly are not correctly handled #739

Open smaret opened 3 years ago

smaret commented 3 years ago

I think there is a problem with how masks are handled for cubes that are created directly. For example, this:

import numpy as np
from spectral_cube import SpectralCube
from spectral_cube import BooleanArrayMask
from astropy.wcs import WCS

data = np.zeros((51, 101, 101))
w = WCS(naxis=3)
w.wcs.crpix = [51, 51, 26]
w.wcs.crval = [0, 0, 0]
w.wcs.cdelt = [.1, .1, .5]
w.wcs.cunit = ['deg', 'deg', 'km/s']
w.wcs.ctype = ['RA---CAR', 'DEC--CAR', 'VELO-LSR']

cube = SpectralCube(data, w)
print(cube)

gives the following output, which looks fine:

SpectralCube with shape=(51, 101, 101):
 n_x:    101  type_x: RA---CAR  unit_x: deg    range:     5.000000 deg:  355.000000 deg
 n_y:    101  type_y: DEC--CAR  unit_y: deg    range:    -5.000000 deg:    5.000000 deg
 n_s:     51  type_s: VOPT      unit_s: m / s  range:   -12500.000 m / s:   12500.000 m / s

However accessing the cube values:

cube[:, 0, 0]

raises a TypeError:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "spectral_cube/spectral_cube.py", line 1219, in __getitem__
    mask=self.mask[view],
TypeError: 'NoneType' object is not subscriptable

Indeed, the cube does not have a mask:

print(cube.mask)
None

To solve this issue, one need to attach a mask to the datacube:

mask = BooleanArrayMask(np.ones_like(data, dtype=bool), wcs=w)
cube = SpectralCube(data, w, mask=mask)
cube[:, 0, 0]  # this works

I think that the mask should be created automatically when SpectralCube is instantiated with mask = None (which is the default).

keflavich commented 3 years ago

Actually, we're supposed to have all of the mask-using tools check before attempting to index. We do this now: https://github.com/radio-astro-tools/spectral-cube/blob/master/spectral_cube/spectral_cube.py#L1324

so it looks like you're on an old version of spectral-cube. I recommend updating to the latest version, and I think we need to issue a new release.

smaret commented 3 years ago

You're right: I'm using 0.4.5. I though it was the latest version because it is tagged as such on GitHub. I didn't notice that there was a more recent version on PyPi. Upgrading to 0.5.0 (from PyPI) solves this problem. Thanks!

smaret commented 3 years ago

I've found another error in 0.5.0 which seems related:

cube = SpectralCube(data, w)

# Attach a beam that correspond to the pixel area.
pixel_area = abs(cube.wcs.wcs.cdelt[0]) * abs(cube.wcs.wcs.cdelt[1]) * u.deg ** 2
beam_sigma = np.sqrt(pixel_area / (2 * np.pi))
beam_fwhm = beam_sigma * np.sqrt(8 * np.log(2))
cube = cube.with_beam(Beam(beam_fwhm))

# Convolve with the instrument beam
beam = Beam(major=1 * u.deg, minor=0.5 * u.deg, pa=38.7 * u.deg)
cube = cube.convolve_to(beam)
Traceback (most recent call last):
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/joblib/parallel.py", line 820, in dispatch_one_batch
    tasks = self._ready_batches.get(block=False)
  File "/nix/store/kz0lkr3nvqbzq7idiil10cak77prxah6-python3-3.8.9/lib/python3.8/queue.py", line 167, in get
    raise Empty
_queue.Empty

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "test.py", line 33, in <module>
    cube = cube.convolve_to(beam)
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/spectral_cube/utils.py", line 49, in wrapper
    return function(self, *args, **kwargs)
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/spectral_cube/spectral_cube.py", line 3130, in convolve_to
    newcube = self.apply_function_parallel_spatial(convfunc,
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/spectral_cube/spectral_cube.py", line 2846, in apply_function_parallel_spatial
    return self._apply_function_parallel_base(images, function,
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/spectral_cube/spectral_cube.py", line 2765, in _apply_function_parallel_base
    Parallel(n_jobs=num_cores,
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/joblib/parallel.py", line 1041, in __call__
    if self.dispatch_one_batch(iterator):
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/joblib/parallel.py", line 831, in dispatch_one_batch
    islice = list(itertools.islice(iterator, big_batch_size))
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/spectral_cube/spectral_cube.py", line 2767, in <genexpr>
    max_nbytes=None)(delayed(applicator)(arg, outcube,
  File "/Users/smaret/Documents/Software/sc-bug/.direnv/python-3.8.9/lib/python3.8/site-packages/spectral_cube/spectral_cube.py", line 2841, in <genexpr>
    self.mask.include(view=(ii, slice(None), slice(None))),
AttributeError: 'NoneType' object has no attribute 'include'

The same cube with a mask:

mask = BooleanArrayMask(np.ones_like(data, dtype=bool), wcs=w)
cube = SpectralCube(data, w, mask=mask)

gives no error.

keflavich commented 3 years ago

This has been sitting untested for a while. It looks like the problem is that cubes initialized with masks can still exist and have their mask attributes called. We need to develop a MWE test and a fix.