Upgrade to python 3.9.13 breaks peak_sn calculation

[jmangum]

Upgraded from python 3.7.13 to 3.9.13 with commensurate update of all conda packages. Appears now that quantity is less-tolerant of calculating the ratio of two images when they are not the same dimension:

In [4]: 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]
Sample Pixel =  (180, 178)
 Sample Pixel Type =  <class 'tuple'>
{'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_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': <Quantity [88.6316022] GHz>, 'my_line_widths': <Quantity [150.] km / s>, 'my_line_names': ['HCN_10'], 'signal_mask_limit': 3, 'spatial_mask_limit': 3, 'sample_pixel': (180, 178)}
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 0x119c1b160>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.argmax at 0x119c1b820>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x119c1b4c0>) 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: VarianceWarning: Note that the second moment returned will be a variance map. To get a linewidth map, use the SpectralCube.linewidth_fwhm() or SpectralCube.linewidth_sigma() methods instead. [spectral_cube.spectral_cube]
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x119c1b160>) 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.0012017893604934216 Jy / beam
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.std at 0x119c1b160>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
Sample Pixel =  (180, 178)
 Sample Pixel Type =  <class 'tuple'>
INFO: Line: HCN_10, 88.6316022 GHz, 150.0 km / s [__main__]
/Users/jmangum/Python/mangum_galaxies-master/CubeLineMoment.py:399: RuntimeWarning: divide by zero encountered in divide
  gauss_mask_cube = np.exp(-(np.array(centroid_map)[None,:,:] -
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
File ~/Python/mangum_galaxies-master/CubeLineMoment.py:817, in <module>
    813     return locals()
    816 if __name__ == "__main__":
--> 817     new_locals = main()
    818     locals().update(new_locals)

File ~/Python/mangum_galaxies-master/CubeLineMoment.py:767, in main()
    762 (cube, spatialmaskcube, spatial_mask, noisemap, noisemapbright,
    763  centroid_map, width_map, max_map, peak_velocity) = cubelinemoment_setup(**params)
    765 params.pop('cube')
--> 767 cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
    768                          peak_velocity=peak_velocity,
    769                          centroid_map=centroid_map, max_map=max_map,
    770                          noisemap=noisemap, width_map=width_map,
    771                          regionlabel=regionlabel,
    772                          fit=False, **params)
    774 # params.pop('signal_mask_limit')
    775 # cubelinemoment_multiline(cube=cube, spatial_mask=spatial_mask,
    776 #                          peak_velocity=peak_velocity,
   (...)
    800
    801 # Clean up open figures
    802 pl.close('all')

File ~/Python/mangum_galaxies-master/CubeLineMoment.py:402, in cubelinemoment_multiline(cube, peak_velocity, centroid_map, max_map, noisemap, signal_mask_limit, my_line_list, my_line_widths, my_line_names, target, spatial_mask, width_map, regionlabel, width_map_scaling, width_cut_scaling, fit, apply_width_mask, sample_pixel, **kwargs)
    398 assert centroid_map.unit.is_equivalent(u.km/u.s)
    399 gauss_mask_cube = np.exp(-(np.array(centroid_map)[None,:,:] -
    400                            np.array(subcube.spectral_axis)[:,None,None])**2 /
    401                          (2*np.array(width_map*width_map_scaling)[None,:,:]**2))
--> 402 peak_sn = max_map / noisemap
    404 print("Peak S/N: {0}".format(np.nanmax(peak_sn)))
    406 # threshold at the fraction of the Gaussian corresponding to our peak s/n.
    407 # i.e., if the S/N=6, then the threshold will be 6-sigma
    408 # (this can be modified as you see fit)

File ~/anaconda3/envs/python39/lib/python3.9/site-packages/astropy/units/quantity.py:1148, in Quantity.__truediv__(self, other)
   1145     except UnitsError:  # let other try to deal with it
   1146         return NotImplemented
-> 1148 return super().__truediv__(other)

File ~/anaconda3/envs/python39/lib/python3.9/site-packages/astropy/units/quantity.py:611, in Quantity.__array_ufunc__(self, function, method, *inputs, **kwargs)
    608     arrays.append(converter(input_) if converter else input_)
    610 # Call our superclass's __array_ufunc__
--> 611 result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    612 # If unit is None, a plain array is expected (e.g., comparisons), which
    613 # means we're done.
    614 # We're also done if the result was None (for method 'at') or
    615 # NotImplemented, which can happen if other inputs/outputs override
    616 # __array_ufunc__; hopefully, they can then deal with us.
    617 if unit is None or result is None or result is NotImplemented:

ValueError: operands could not be broadcast together with shapes (800,800) (334,400)

If this is due to a change in how astropy.units.quantity works, perhaps a simple modification of the

peak_sn = max_map / noisemap

would fix this bug?

[keflavich]

No, you could never subtract two arrays with different shapes. I'm not sure what's causing this.

[keflavich]

the problem is that the max_map is 800x800 and the noise_map is 334x400. Why are they different?

[jmangum]

Don't know. Will look closer.

[jmangum]

Zeroing-in on why max_map has dimension 800x800 when it should be the same as noise_map (334x400). By switching back-and-forth between my python 3.7 and python 3.9 environments, have found that the issue appears to be how spectral_cube treats NaN slices.

For python 3.7, which gets the dimensioning right, I see just before the max_map calculation:

/Users/jmangum/anaconda3/lib/python3.7/site-packages/spectral_cube/spectral_cube.py:443: RuntimeWarning: All-NaN slice encountered
  **kwargs)
brightest_cube.shape =  (81, 334, 400)
max_map.shape =  (334, 400)

...while for python 3.9 I see the following at the same place in the script:

/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,
brightest_cube.shape =  (81, 800, 800)
max_map.shape =  (800, 800)

The max_map calculation is the following for both:

    # Create a copy of the cutoutcube with velocity units
    cutoutVcube = cutoutcube.with_spectral_unit(u.km/u.s,
                                                   rest_value=brightest_line_frequency,
                                                   velocity_convention='optical')

    # Use the brightest line to identify the appropriate peak velocities, but ONLY
    # from a slab including +/- width:
    brightest_cube = cutoutVcube.spectral_slab(vz-velocity_half_range,
                                               vz+velocity_half_range)

    # compute various moments & statistics along the spcetral dimension
    peak_velocity = brightest_cube.spectral_axis[brightest_cube.argmax(axis=0)]
    max_map = peak_amplitude = brightest_cube.max(axis=0)

Guessing that it is the cutoutVcube.spectral_slab calculation that is not handling NaN well. This does not appear to be quite right, though, as I am using the same version of spectral-cube in both cases (version 0.6.0).

[jmangum]

Found problem. Typo in recently-added try-except change to address deprecation of regions.read_ds9. Runs now but with another error. Closing this issue.