SpectralCube use of `.max` on cube with NaNs needs modification

[jmangum]

The following line:

max_map = peak_amplitude = brightest_cube.max(axis=0)

...produces an all-NaN slice warning under some circumstances. I think that this could be fixed with the use of a function like numpy.nanmax. How does one do this with SpectralCube?

[keflavich]

SpectralCube always uses nanmax. This is saying that you're trying to do np.nanmax([nan, nan, nan]), which will still return nan. You can silence the warning?

[jmangum]

Thanks. Probably should not silence as I don't think that brightest_cube.max(axis=0) should be an all-NaN slice. Looking for a bug here as brightest_cube in this case is a cutoutcube that is not the same as `cube' (change implemented last week).

[jmangum]

Hmmmm. Tested simple use case where cube and cutoutcube are the same and still get this warning. Still searching for source of all-NaN slice in cutoutcube.

[jmangum]

Now wondering if there is a bug in the max_map calculation. The pertinent section of code (with some debugging print statements thrown-in) is:

    # 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)
    print('vz: ',vz)
    print('velocity_half_range: ',velocity_half_range)
    print('brightest_cube spectral axis: ',brightest_cube.spectral_axis)
    print('brightest_cube unmasked values: ',brightest_cube.unmasked_data[0,:,:])

    # compute various moments & statistics along the spcetral dimension
    peak_velocity = brightest_cube.spectral_axis[brightest_cube.argmax(axis=0)]
    print('Peak Velocity: ',peak_velocity)
    print('All NaN Slice Here?')
    max_map = peak_amplitude = brightest_cube.max(axis=0) # This is sometimes as all-NaN slice
    print('max_map: ',max_map)
    print('...maybe not...')
    width_map = brightest_cube.linewidth_sigma() # or vcube.moment2(axis=0)**0.5
    centroid_map = brightest_cube.moment1(axis=0)

This produces the following output:

vz:  236.0 km / s
velocity_half_range:  400.0 km / s
brightest_cube spectral axis:  [ 638.82525781  628.79156544  618.75854325  608.72619118  598.69450915
  588.66349709  578.63315496  568.60348266  558.57448015  548.54614734
  538.51848418  528.49149059  518.46516652  508.43951189  498.41452663
  488.39021068  478.36656398  468.34358644  458.32127802  448.29963863
  438.27866822  428.25836672  418.23873405  408.21977016  398.20147497
  388.18384843  378.16689045  368.15060098  358.13497994  348.12002728
  338.10574292  328.0921268   318.07917885  308.066899    298.05528719
  288.04434335  278.03406741  268.02445931  258.01551898  248.00724634
  237.99964135  227.99270392  217.986434    207.98083151  197.97589639
  187.97162857  177.96802798  167.96509456  157.96282825  147.96122897
  137.96029665  127.96003124  117.96043267  107.96150086   97.96323575
   87.96563727   77.96870537   67.97243996   57.97684099   47.98190839
   37.98764209   27.99404202   18.00110812    8.00884032   -1.98276144
  -11.97369723  -21.96396713  -31.95357119  -41.94250948  -51.93078206
  -61.91838902  -71.9053304   -81.89160628  -91.87721673 -101.86216181
 -111.84644158 -121.83005612 -131.81300549 -141.79528975 -151.77690898
 -161.75786324] km / s
brightest_cube unmasked values:  [[ 3.96811520e-05 -8.41866713e-05 -7.88558391e-05 ... -3.87204345e-04
  -5.62617672e-04 -5.18065877e-04]
 [ 4.96062567e-05 -4.27419436e-05 -1.57168834e-05 ... -5.08279598e-04
  -6.71427115e-04 -7.07802130e-04]
 [ 5.01758186e-05  7.55995279e-06 -6.29120623e-05 ... -5.45957300e-04
  -7.76681176e-04 -7.66127894e-04]
 ...
 [ 4.02163831e-04  3.15971964e-04  1.90121267e-04 ...  3.12264514e-04
   4.26833227e-04  6.56572112e-04]
 [ 4.44979232e-04  2.97025108e-04  1.76950940e-04 ...  1.40180746e-05
   2.56236293e-04  4.44111705e-04]
 [ 2.50054116e-04  1.91611151e-04 -1.24303915e-05 ...  8.32177466e-06
   1.18798533e-04  3.42549465e-04]] Jy / beam
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.argmax at 0x19e905a60>) requires loading the entire cube into memory and may therefore be slow. [spectral_cube.utils]
Peak Velocity:  [[638.82525781 638.82525781 638.82525781 ... 638.82525781 638.82525781
  638.82525781]
 [638.82525781 638.82525781 638.82525781 ... 638.82525781 638.82525781
  638.82525781]
 [638.82525781 638.82525781 638.82525781 ... 638.82525781 638.82525781
  638.82525781]
 ...
 [638.82525781 638.82525781 638.82525781 ... 638.82525781 638.82525781
  638.82525781]
 [638.82525781 638.82525781 638.82525781 ... 638.82525781 638.82525781
  638.82525781]
 [638.82525781 638.82525781 638.82525781 ... 638.82525781 638.82525781
  638.82525781]] km / s
All NaN Slice Here?
WARNING: PossiblySlowWarning: This function (<function BaseSpectralCube.max at 0x19e905700>) 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,
max_map:  [[nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 ...
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]
 [nan nan nan ... nan nan nan]] Jy / beam
...maybe not...

Not at all clear to me why the max_map calculation is producing an all-NaN array when the brightest_cube appears to have reasonable values.

[jmangum]

I think that I understand now (though there still might be a bug). The following:

    max_map = peak_amplitude = brightest_cube.max(axis=0) # This is sometimes as all-NaN slice
    print('shape(max_map): ',max_map.shape)
    print('type(max_map): ',type(max_map))
    print('max_map.max (all dims): ',max_map.max())
    print('np.nanmax(max_map) (all dims): ',np.nanmax(max_map))

...produces the following output:

/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,
shape(max_map):  (334, 400)
type(max_map):  <class 'spectral_cube.lower_dimensional_structures.Projection'>
max_map.max (all dims):  nan Jy / beam
np.nanmax(max_map) (all dims):  0.16510973870754242 Jy / beam

This makes me suspect that two things are (possibly) happening: (1) The calculation of max_map, which produces the maximum spectral value for each pixel position, is correct. Some pixel positions are completely blanked, which means that all spectral values are NaN, while others have a mixture of real and NaN values. The warning message says that there was an All-NaN slice encountered, meaning that at least one spectrum toward a pixel was all-NaN, which is correct. (2) @keflavich Why isn't max_map.max() the same as np.nanmax(max_map)?

Found an additional issue that I will open a separate issue for... max_map is derived from a masked version of cutoutcube, which is the so-called "brightest line cube". max_map is used later in cubelinemoment_multiline to calculate peak_sn = max_map / noisemap. noisemap is derived from cube (the target line cube), not from cutoutcube. This looks inconsistent...

[keflavich]

(2) @keflavich Why isn't max_map.max() the same as np.nanmax(max_map)?

max_map.max() is np.max(array) because max_map is not a SpectralCube, it's a Projection

[jmangum]

ok. Can close this issue.