SpectralCube use of `.max` on cube with NaNs needs modification
jmangum opened this issue · 7 comments
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?
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?
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).
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.
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.
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...
(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
ok. Can close this issue.