Setting `signal_mask_limit` and `spatial_mask_limit` to `None` Does Not Result in an Unmasked Image

[jmangum]

Using Bayesian inference chemical models with integrated intensities has led us to believe that integrated intensities with no signal mask limitations should be used. Expected that setting signal_mask_limit and spatial_mask_limit to None should have resulted in an integrated intensity image with only velocity range masking. Instead I get an image that looks very much like an integrated intensity image with 3-sigma masking. Tried setting the mask_negatives boolean in the cubelinement_setup call to False, but that only partially solved the problem. Still get masked pixels. Searched through script to see if I could spot obvious error but found none. @keflavich designed the masking scheme so perhaps this is not a bug but a feature? Should be able to just integrate over a velocity range and get an unmasked image, though.

[keflavich]

I'll have a look.

[keflavich]

https://github.com/keflavich/cube-line-extractor/blob/33a7688aa72a94a3572573f078899c6a9935b63c/CubeLineMoment.py#L176-L179 and https://github.com/keflavich/cube-line-extractor/blob/33a7688aa72a94a3572573f078899c6a9935b63c/CubeLineMoment.py#L259-L262 are the only places where any comparison to the flux is done.

How sure are you that extra masking is being done? Best thing to do is test on a single sample pixel

[jmangum]

Thank you @keflavich. The moment0 image contains NaNs in the regions where it would normally mask noise, so something is masking when I try to not signal mask. Should not be velocity range masking either as single spectral channels should still have noise over entire image.

I am not sure how I would test on a single sample pixel.

[keflavich]

nans are always masked out, they can't be un-masked. Are there any valid pixels in the velocity range you're looking at?

And can you clarify the specific symptom you've identified: is is that there are NaNs in the final output moment map in spatial pixels you do not expect them? Or is it that the flux values in some pixels in the final moment map are higher than you expect?

[jmangum]

There are no NaNs in the input image cube (within the primary beam corrected region of the image). The symptom is the first thing you mention. It is masking (specifically, setting pixel values to NaN) in the final output moment map in spatial pixels that I would expect to see non-NaN values.

In case you want to try this test yourself, I have posted the necessary files to astrocloud. Will send link separately. You will need to edit the YAML file to update cube and cutoutcube appropriately.

[keflavich]

ok I've run it locally. Do you expect to see no nan pixels?

[keflavich]

OK, this is a feature: the problem is, moment2 returns negative values, so going from moment2 -> width is sqrt(negative). That's why the moment maps have NaNs in them.

[jmangum]

Where is it calculating width from sqrt(moment2)?

[keflavich]

https://github.com/keflavich/cube-line-extractor/blob/33a7688aa72a94a3572573f078899c6a9935b63c/CubeLineMoment.py#L207

[jmangum]

Thanks. Shouldn't there be a way to simply integrate along the spectral axis without applying any intensity-level masking?

[keflavich]

This isn't an intensity-level masking. Remember the way we define the spectral region to integrate over is to model the emission in each pixel with a simple Gaussian and perform a threshold on that Gaussian to determine what spectral pixels to include. If we can't measure a width, we can't define that Gaussian.

An option would be to specify some minimal width to plug in if no width can be measured from the data. Then, however, we run into the next problem: What do we adopt for moment1 if it is badly-behaved? Do we just set it to the mean velocity, say? How do we test for bad behavior? I'll propose some answers to these questions.

[jmangum]

Mean velocity seems reasonable for moment1 default. Could also use velocity_half_range for moment2 default if needed.

[keflavich]

I pushed a change in commit 89a98b9. Try adding parameter use_default_width: True and see if that gets what you want.

I still see NaNs in the final moment image. Some of that appears to be because there are some tiny widths, some is... not obvious yet.

[jmangum]

Thanks. Added use_default_width: True to YAML file. Script crashed. Did you mean to add the new parameter at the call?

[jmangum]

In [9]: run ~/Python/mangum_galaxies-master/CubeLineMoment.py HCN10.yaml
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x117f221f0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.argmax at 0x117f228b0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x117f22550>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/spectral_cube/spectral_cube.py:441: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x117f221f0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
noisemapbright peak = 0.0007277460535988212 Jy / beam
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x117f221f0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/Python/mangum_galaxies-master/CubeLineMoment.py:913, in <module>
    909     return locals()
    912 if __name__ == "__main__":
--> 913     new_locals = main()
    914     locals().update(new_locals)

File ~/Python/mangum_galaxies-master/CubeLineMoment.py:863, in main()
    859 params.pop('cube')
    860 params.pop('sample_pixel')
--> 863 cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
    864                          peak_velocity=peak_velocity,
    865                          centroid_map=centroid_map, max_map=max_map,
    866                          noisemap=noisemap, noisemapbright=noisemapbright,
    867                          width_map=width_map, sample_pixel=sample_pixel,
    868                          fit=False, **params)
    870 # params.pop('signal_mask_limit')
    871 # cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
    872 #                          peak_velocity=peak_velocity,
   (...)
    896
    897 # Clean up open figures
    898 pl.close('all')

File ~/Python/mangum_galaxies-master/CubeLineMoment.py:430, in cubelinemoment_multiline(cube, peak_velocity, centroid_map, max_map, noisemap, noisemapbright, signal_mask_limit, spatial_mask_limit, brightest_line_name, brightest_line_frequency, my_line_list, my_line_widths, my_line_names, target, spatial_mask, width_map, sample_pixel, width_map_scaling, width_cut_scaling, use_default_width, fit, apply_width_mask, **kwargs)
    428     bad_widths = np.isnan(width_map)
    429 if np.any(width_map <= 0):
--> 430     raise ValueError("Negative or zero width found in width map")
    432 # Now loop over EACH line, extracting moments etc. from the appropriate region:
    433 # we'll also apply a transition-dependent width (my_line_widths) here because
    434 # these fainter lines do not have peaks as far out as the bright line.
    436 for line_name,line_freq,line_width in zip(my_line_names,my_line_list,my_line_widths):

ValueError: Negative or zero width found in width map

[keflavich]

well that shouldn't happen

[keflavich]

you should specify these arguments in the yaml:

use_default_width: True
min_width: 10

min_width: 10 is ~10 km/s, just a smidge bigger than the channel spacing

[keflavich]

I pushed a new commit (c9829a6); try again

[keflavich]

What was the crash for use_default_width: True, though? That also shouldn't crash.

[jmangum]

Pull latest. Ran fine but still gives me a masked moment0 image.

[keflavich]

Oooooh, it's because we impose a S/N mask on the Gaussian width masking

[jmangum]

Right! I think that it is the following (and note the comment...):

            # NOTE: Following line sometimes produces image with NaNs at some positions.  Should try to fix...
            gauss_mask_cube = np.exp(-(np.array(centroid_map)[None,:,:] -
                                       np.array(subcube.spectral_axis)[:,None,None])**2 /
                                     (2*np.array(width_map*width_map_scaling)[None,:,:]**2))

[keflavich]

try bb44e64

the problem was that the peak velocity (velocity of peak intensity) was being used as the reference but not as the center of the Gaussian, and in low S/N cases, the selected regions by the Gaussian and the simple v-cut methods did not agree

[jmangum]

Thanks. Still masking. Looks very similar to before.

[keflavich]

I now get the above

params are:

cube: ngc253.88632MHz.HCN.1-0.regrid.12mC12mE.image.pbcor.fits
cuberegion: CubeLineMoment.reg
cutoutcube: ngc253.88632MHz.HCN.1-0.regrid.12mC12mE.image.pbcor.fits
cutoutcuberegion: CubeLineMoment.reg
vz: 236
target: NGC253
brightest_line_name: HCN_10
brightest_line_frequency: 88.6316022
width_line_frequency: 88.6316022
velocity_half_range: 400
noisemapbright_baseline: [[0,24],[77,100]]
noisemap_baseline: [[0,24],[77,100]]
my_line_list: 88.6316022
my_line_widths: 150
my_line_names: HCN_10
signal_mask_limit: None
spatial_mask_limit: None
mask_negatives: False
sample_pixel: Region6.reg
use_default_width: True
min_width: 10
min_gauss_threshold: 0.1

[jmangum]

Cannot reproduce. Did another git pull and re-ran with update YAML file (added min_gauss_threshold: 0.1. Now get a completely masked image.

[jmangum]

YAML file:

cube: /Users/jmangum/Science/ALCHEMI/ScienceProjects/HCNHNC/createCubes/ngc253.88632MHz.HCN.1-0.regrid.12mC12mE.image.pbcor.fits
cuberegion: CubeLineMoment.reg
cutoutcube: /Users/jmangum/Science/ALCHEMI/ScienceProjects/HCNHNC/createCubes/ngc253.88632MHz.HCN.1-0.regrid.12mC12mE.image.pbcor.fits
cutoutcuberegion: CubeLineMoment.reg
vz: 236
target: NGC253
brightest_line_name: HCN_10
brightest_line_frequency: 88.6316022
width_line_frequency: 88.6316022
velocity_half_range: 400
noisemapbright_baseline: [[0,24],[77,100]]
noisemap_baseline: [[0,24],[77,100]]
my_line_list: 88.6316022
my_line_widths: 150
my_line_names: HCN_10
signal_mask_limit: None
spatial_mask_limit: None
sample_pixel: Region6.reg
use_default_width: True
min_width: 10
min_gauss_threshold: 0.1

[keflavich]

ok that's not great, I don't see any reason we should suddenly have diverged

[jmangum]

Found the problem. Had not added mask_negatives: False to YAML. Now get the following:

[jmangum]

I think that this issue is resolved. I will open an issue as a reminder to update the documentation to describe the new input parameters and assign to me. Will try to do soon but may be after ALMA proposal deadline. Thanks!

[jmangum]

Unfortunately, further testing suggests that this issue is not yet resolved. Tried to create integrated intensity images with another image cube (HCN 2-1) and found that I get an all-blanked moment0 image when I try to integrated with no signal limits. Also tried to reproduce the moment0 image produced before using a signal_mask_limit and spatial_mask_limit of 3, but that also yielded an all-blanked moment0 image. Have posted image cube and other CubeLineMoment files used to do this HCN 2-1 test to astrocloud.

[keflavich]

I'm traveling today so I might not get around to helping with this.

[jmangum]

For reference, run with HCN 2-1 image cube (posted to astrocloud) produced the following messages:

In [14]: run ~/Python/mangum_galaxies-master/CubeLineMoment.py ../HCN21.yaml
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.argmax at 0x11726baf0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x11726b790>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/spectral_cube/spectral_cube.py:441: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
noisemapbright peak = 0.017179029062390327 Jy / beam
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py:455: UserWarning: Found minimum width in width_map 0.0 km / s.
  warnings.warn(f"Found minimum width in width_map {np.nanmin(width_map)}.")
INFO: Line: HCN_21, 177.2611112 GHz, 150.0 km / s [__main__]
/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py:496: RuntimeWarning: divide by zero encountered in divide
  gauss_mask_cube = np.exp(-(np.array(velocity_map)[None,:,:] -
Peak S/N: 325.4417419433594
Highest Threshold: 0.7893054485321045
Lowest Positive Threshold: 0.0030727465637028217
Lowest Threshold: nan
Sample Pixel =  (180, 178, 'Region 6')
SP Threshold: 0.003982150927186012
SP S, N, S/N: 0.5061072111129761 Jy / beam, 0.0020153953228145838 Jy / beam, 251.1205596923828
Number of values above threshold: 0
Number of spatial pixels excluded: 133600
Max value in the mask cube: 0.3126231329466763
shapes: mask cube=(1, 334, 400)  threshold: (334, 400)
INFO: Auto-setting vmin to  0.000e+00 [aplpy.core]
INFO: Auto-setting vmax to  0.000e+00 [aplpy.core]
Moment 0 for sample pixel Region 6 is nan Jy km / (beam s)
INFO: Auto-setting vmin to  0.000e+00 [aplpy.core]
INFO: Auto-setting vmax to  0.000e+00 [aplpy.core]
Moment 1 for sample pixel Region 6 is nan km / s
INFO: Auto-setting vmin to  0.000e+00 [aplpy.core]
INFO: Auto-setting vmax to  0.000e+00 [aplpy.core]
Moment 2 for sample pixel Region 6 is nan km / s

[jmangum]

Decided to try some more tests thinking that they might help isolate the source of the failures. Since we know that our HCN 1-0 test was successful and that HCN 2-1 failed, did some more HCN/HNC testing where I already have 3-sigma-clipped integrated intensity images from a previous analysis (so I know that all of these species/transitions should produce usable integrated intensities). Here is what I found:

• HCN 1-0: Good
• HCN 2-1: Fail
• HCN 3-2: Fail
• HCN 4-3: Script crash...see below
• HNC 1-0: Good
• HNC 2-1: Fail
• HNC 3-2: Fail
• HNC 4-3: Fail

Note that the "Fail" were all the same. Produced an all-blanked moment0. The HCN 4-3 script crash was different:

In [20]: run ~/Python/mangum_galaxies-master/CubeLineMoment.py HCN43.yaml
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.argmax at 0x11726baf0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x11726b790>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/spectral_cube/spectral_cube.py:441: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
noisemapbright peak = 0.01124945655465126 Jy / beam
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
INFO: Line: HCN_43, 354.5054773 GHz, 130.0 km / s [__main__]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/astropy/units/quantity.py:611: RuntimeWarning: divide by zero encountered in divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
Peak S/N: inf
Highest Threshold: 3.26678466796875
Lowest Positive Threshold: 0.10000000149011612
Lowest Threshold: nan
Sample Pixel =  (180, 178, 'Region 6')
SP Threshold: 0.10000000149011612
SP S, N, S/N: 1.3806281089782715 Jy / beam, 0.0005810492439195514 Jy / beam, 2376.0947265625
Number of values above threshold: 585
Number of spatial pixels excluded: 133015
Max value in the mask cube: 0.9472874424110594
shapes: mask cube=(1, 334, 400)  threshold: (334, 400)
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
File ~/Python/mangum_galaxies-master/CubeLineMoment.py:959, in <module>
    955     return locals()
    958 if __name__ == "__main__":
--> 959     new_locals = main()
    960     locals().update(new_locals)

File ~/Python/mangum_galaxies-master/CubeLineMoment.py:909, in main()
    905 params.pop('cube')
    906 params.pop('sample_pixel')
--> 909 cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
    910                          peak_velocity=peak_velocity,
    911                          centroid_map=centroid_map, max_map=max_map,
    912                          noisemap=noisemap, noisemapbright=noisemapbright,
    913                          width_map=width_map, sample_pixel=sample_pixel,
    914                          fit=False, **params)
    916 # params.pop('signal_mask_limit')
    917 # cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
    918 #                          peak_velocity=peak_velocity,
   (...)
    942
    943 # Clean up open figures
    944 pl.close('all')

File ~/Python/mangum_galaxies-master/CubeLineMoment.py:553, in cubelinemoment_multiline(cube, peak_velocity, centroid_map, max_map, noisemap, noisemapbright, signal_mask_limit, spatial_mask_limit, brightest_line_name, brightest_line_frequency, my_line_list, my_line_widths, my_line_names, target, spatial_mask, width_map, sample_pixel, width_map_scaling, width_cut_scaling, use_default_width, fit, apply_width_mask, min_width, min_gauss_threshold, use_peak_for_velcut, **kwargs)
    551     spatially_masked_pixels = (width_mask_cube.max(axis=0) == 0)
    552     spatially_masked_pixels2 = msubcube.mask.include().max(axis=0) == 0
--> 553     assert np.all(spatially_masked_pixels == spatially_masked_pixels2)
    554     assert spatially_masked_pixels.sum() == spatially_masked_pixels2.sum()
    556 # this part makes a cube of velocities

AssertionError:

Unfortunately I don't see anything inherently different between the HCN 4-3 and the other image cubes.

[keflavich]

problem was NaNs in the velocity map. fixed in 3979277

[jmangum]

Thanks @keflavich. Tested update with all eight HCN/HNC combinations. Six produce reasonable results (still a few blanks but probably usable). Two, HNC21 and HNC43, failed:

HNC21:

In [29]: run ~/Python/mangum_galaxies-master/CubeLineMoment.py HNC21.yaml
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.argmax at 0x11726baf0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x11726b790>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/spectral_cube/spectral_cube.py:441: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
noisemapbright peak = 0.05820803344249725 Jy / beam
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
INFO: Line: HNC_21, 181.324758 GHz, 130.0 km / s [__main__]
Peak S/N: 58.462646484375
Highest Threshold: 1.1932094097137451
Lowest Positive Threshold: 0.10000000149011612
Lowest Threshold: 0.10000000149011612
Sample Pixel =  (180, 178, 'Region 6')
SP Threshold: 0.10000000149011612
SP S, N, S/N: 0.35602399706840515 Jy / beam, 0.009534697979688644 Jy / beam, 37.33982849121094
Number of values above threshold: 5960018
Number of spatial pixels excluded: 5282
Max value in the mask cube: 0.9999999999999913
shapes: mask cube=(101, 334, 400)  threshold: (334, 400)
Unexpected exception formatting exception. Falling back to standard exception
Traceback (most recent call last):
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/interactiveshell.py", line 2738, in safe_execfile
    py3compat.execfile(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/utils/py3compat.py", line 55, in execfile
    exec(compiler(f.read(), fname, "exec"), glob, loc)
  File "/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py", line 980, in <module>
    new_locals = main()
  File "/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py", line 930, in main
    cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
  File "/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py", line 574, in cubelinemoment_multiline
    assert np.all(spatially_masked_pixels == spatially_masked_pixels2)
AssertionError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/interactiveshell.py", line 1993, in showtraceback
    stb = self.InteractiveTB.structured_traceback(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 1118, in structured_traceback
    return FormattedTB.structured_traceback(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 1012, in structured_traceback
    return VerboseTB.structured_traceback(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 865, in structured_traceback
    formatted_exception = self.format_exception_as_a_whole(etype, evalue, etb, number_of_lines_of_context,
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 818, in format_exception_as_a_whole
    frames.append(self.format_record(r))
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 736, in format_record
    result += ''.join(_format_traceback_lines(frame_info.lines, Colors, self.has_colors, lvals))
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/utils.py", line 145, in cached_property_wrapper
    value = obj.__dict__[self.func.__name__] = self.func(obj)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/core.py", line 698, in lines
    pieces = self.included_pieces
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/utils.py", line 145, in cached_property_wrapper
    value = obj.__dict__[self.func.__name__] = self.func(obj)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/core.py", line 649, in included_pieces
    pos = scope_pieces.index(self.executing_piece)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/utils.py", line 145, in cached_property_wrapper
    value = obj.__dict__[self.func.__name__] = self.func(obj)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/core.py", line 628, in executing_piece
    return only(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/executing/executing.py", line 164, in only
    raise NotOneValueFound('Expected one value, found 0')
executing.executing.NotOneValueFound: Expected one value, found 0

HNC43:

In [31]: run ~/Python/mangum_galaxies-master/CubeLineMoment.py HNC43.yaml
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: StokesWarning: Cube is a Stokes cube, returning spectral cube for I component [spectral_cube.io.core]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.argmax at 0x11726baf0>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x11726b790>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/spectral_cube/spectral_cube.py:441: RuntimeWarning: All-NaN slice encountered
  out = function(self._get_filled_data(fill=fill,
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
noisemapbright peak = 0.017722057178616524 Jy / beam
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x11726b430>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1878: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,
INFO: Line: HNC_43, 362.630303 GHz, 100.0 km / s [__main__]
Peak S/N: 620.5407104492188
Highest Threshold: 3.4910507202148438
Lowest Positive Threshold: 0.10000000149011612
Lowest Threshold: 0.10000000149011612
Sample Pixel =  (180, 178, 'Region 6')
SP Threshold: 0.10000000149011612
SP S, N, S/N: 0.8129488825798035 Jy / beam, 0.0025105164386332035 Jy / beam, 323.8173828125
Number of values above threshold: 2877240
Number of spatial pixels excluded: 36656
Max value in the mask cube: 0.9999999999997937
shapes: mask cube=(71, 334, 400)  threshold: (334, 400)
Unexpected exception formatting exception. Falling back to standard exception
Traceback (most recent call last):
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/interactiveshell.py", line 2738, in safe_execfile
    py3compat.execfile(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/utils/py3compat.py", line 55, in execfile
    exec(compiler(f.read(), fname, "exec"), glob, loc)
  File "/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py", line 980, in <module>
    new_locals = main()
  File "/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py", line 930, in main
    cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
  File "/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py", line 574, in cubelinemoment_multiline
    assert np.all(spatially_masked_pixels == spatially_masked_pixels2)
AssertionError

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/interactiveshell.py", line 1993, in showtraceback
    stb = self.InteractiveTB.structured_traceback(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 1118, in structured_traceback
    return FormattedTB.structured_traceback(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 1012, in structured_traceback
    return VerboseTB.structured_traceback(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 865, in structured_traceback
    formatted_exception = self.format_exception_as_a_whole(etype, evalue, etb, number_of_lines_of_context,
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 818, in format_exception_as_a_whole
    frames.append(self.format_record(r))
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/IPython/core/ultratb.py", line 736, in format_record
    result += ''.join(_format_traceback_lines(frame_info.lines, Colors, self.has_colors, lvals))
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/utils.py", line 145, in cached_property_wrapper
    value = obj.__dict__[self.func.__name__] = self.func(obj)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/core.py", line 698, in lines
    pieces = self.included_pieces
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/utils.py", line 145, in cached_property_wrapper
    value = obj.__dict__[self.func.__name__] = self.func(obj)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/core.py", line 649, in included_pieces
    pos = scope_pieces.index(self.executing_piece)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/utils.py", line 145, in cached_property_wrapper
    value = obj.__dict__[self.func.__name__] = self.func(obj)
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/stack_data/core.py", line 628, in executing_piece
    return only(
  File "/Users/jmangum/anaconda3/envs/python39/lib/python3.9/site-packages/executing/executing.py", line 164, in only
    raise NotOneValueFound('Expected one value, found 0')
executing.executing.NotOneValueFound: Expected one value, found 0

Posted HNC21 and HNC43 files to astrocloud.

[jmangum]

I have been experimenting with the debug option and have yet to find an instance when the following if-statement is true (i.e. the block inside the if-statement is never executed):

https://github.com/keflavich/cube-line-extractor/blob/397927791f254e88ac21a96044d0bef13dbab908/CubeLineMoment.py#L227-L231

Not sure this has anything to do with the remaining issue(s), but seemed odd as I thought that the previous problem was related to bad velocity widths toward noisy measurements.

[keflavich]

Could you put up the different .yml files you used in this post? https://github.com/keflavich/cube-line-extractor/issues/43#issuecomment-1480417238

[keflavich]

nvm, I see them all

[keflavich]

HNCs look fine now

[jmangum]

@keflavich : Have encountered another instance of blanking when told not to. I am trying to generate moment images with no clipping (as before). As noted above we got this working for some ALCHEMI HCN, HNC, and isotopomer data. Trying it now on other ALCHEMI image cubes (specifically, H3Op and SO). After producing a yaml file with the appropriate new keywords, I get a moment0 image with blanks, which should not be there. Uploaded the yaml file and associated image cubes to astrocloud (same folder as before, so you will see the old HCN files there also). It is H3Op 1(11)-2(10) that I am trying to produce no-clip moment images for.

[keflavich]

I get the above. Is this consistent with what you're getting? I think these are the same sorts of images we got last time around; there are a few NaNs where there are no data at all, I suspect.

[jmangum]

Thank you @keflavich. Yes, this is what I get. I don't see how, though, there could ever be a blanked pixel when the signal-to-noise cutoff is set to 0. The current images may be exhibiting the same behaviour as the moment images with just a few blanks encountered earlier, but I now don't see how any moment image could have any blanked pixels when we have asked for no signal blanking. Shouldn't the default assumed minimum line width, along with the signal level of a given pixel, ensure that we always get moment values at all pixels?

[keflavich]

This is the behavior I would expect. You are using H3Op as the reference "brightest" line, but it contains no signal in those pixels:

If you want to force extraction at this position, I've added a new parameter max_gauss_threshold in f53e7f4. Set that to, say, 0.9 or something to force inclusion of some spectral pixels along all spatial pixels.

In the example case, the selected pixel had peak S/N = 0.6, which resulted in a threshold = 1.6. Because our Gaussians are peak-normalized, 1 < 1.6, so all pixels were excluded.

[jmangum]

Thank you @keflavich. I understand the logic, which is a bit different than the way I was thinking about how the brightest_line is used in CubeLineMoment. I have always thought of it as a "guide" to where emission should be, while before adding max_gauss_threshold the brightest_line was being treated more like a rule. This was all good when we also wanted to apply a signal-to-noise cutoff. Now that we are using Bayesian inference modeling, our thinking is that applying a signal-to-noise cutoff no longer makes sense, since the model result diagnostic process will apply low-weight to any solutions that derive from noise in an integrated intensity image. I think that the new max_gauss_threshold allows us to create moment images which allow us to completely ignore any signal thresholds and just integrate over a spectral line. Initial testing with H3Op 1(11)-2(10) look great. Will test some more. Thanks!

[jmangum]

Should have posted H3Op 1(11)-2(10) integrated intensity using max_gauss_threshold = 0.9. Attached. . Just a few blanked pixels.

[keflavich]

yeah I don't know what those last few outliers are, but can you safely ignore them?

[jmangum]

Let's close this issue as its original intent was to resolve some blanking issues when attempting to use the no-clip options within CubeLineMoment. I believe that those issues are largely resolved. Note, though, that @keflavich and I have concluded that using these no-clip options to produce an unclipped moment0 image might not be the best approach. Decided to leave the no-clip approach as an option, but to document best-practices. Will do this as part of documentation upgrade.