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

Add a "spectral_only=True" option for minimal_subcube() ? #824

Open apmtowner opened 2 years ago

apmtowner commented 2 years ago

Short version: For https://spectral-cube.readthedocs.io/en/latest/api/spectral_cube.SpectralCube.html#spectral_cube.SpectralCube.minimal_subcube, there is a spatial_only=True boolean option. It would be extremely useful to my (admittedly weird) use case to have a spectral_only=True option as well, i.e. find the minimal subcube in the spectral dimension only and ignore the two spatial dimensions. This is not currently an option but would be very nice to have.

Longer context: I'm using SpectralCube masking to create subcubes where the frequencies that are cut out vary depending on the pixel. In my case, I want to cut out the central 5-10 km/s of the cube, but the vlsr varies by up to 5 km/s across the field of view. This is annoying for my science. I'm using masking to tie the central velocity of each pixel to the actual vlsr of that pixel instead of using the same (mostly wrong) vlsr for the whole field. The basic code looks like this:

mycube = SpectralCube.read('mycube.fits') #the data cube I care about. In my case this is SiO but could be anything.
myreference = SpectralCube.read('myreference.fits') #a different data cube (a good dense gas tracer, in my case) that tells me how the VLSR varies across the field. 

reference_mom1 = myreference.moment(order=1) #make a mom1 map of the reference molecule, to get the vlsr for each pixel

lowerlimitmask = mycube.spectral_axis[:,None,None] > reference_mom1 - mylinewidth #where 'mylinewidth' is 5 km/s, or 10 km/s, or whatever I need cut out
upperlimitmask = mycube.spectral_axis[:,None,None] < reference_mom1 + mylinewidth
combinedmask = lowerlimitmask & upperlimitmask

mycutout = mycube.with_mask(combinedmask)

I had previously been using mycutout = mycube.with_mask(combinedmask).minimal_subcube() instead, because hey, why keep 400 channels in memory if you only need 80? But unfortunately minimal_subcube() was also cutting off exactly one column in my n_x axis as well. I don't know why and I can't make it stop. One good (but currently nonexistent) workaround that doesn't torture the RAM on my machine would be to let minimal_subcube have a spectral_only=True option just as it does for spatial_only.

I freely admit I might be tormenting SpectralCube with purposes it wasn't designed for, but I do think having a spectral_only option would be very helpful in general.

keflavich commented 2 years ago

This should be pretty easy to implement; we can add it.

The way to do this (what we'll implement under the hood, I expect) is:

cube = <your cube>
spectral_inclusion = np.any(cube.mask.include(), axis=(1,2))
start_zpix = np.argmax(spectral_inclusion)
end_zpix = spectral_inclusion.size - np.argmax(spectral_inclusion[::-1])
cube = cube[start_zpix:end_zpix,:,:]
apmtowner commented 2 years ago

Thanks!