[BUG] Cubeviz: Moment 2 fails for NIRSpec IFU cube processed with JWST pipeline

[PatrickOgle]

Describe the bug
When attempting to run the Moment Maps plugin on NIRSpec IFU data, the map is not created for Moment = 2, but the plugin works for Moment = 1, 0

To Reproduce
Steps to reproduce the behavior:

1. Load a NIRSpec IFU cube into Cubeviz
2. Run moment map plugin on the [SCI] extension
3. Enter '2' for Moment (With or without Spectral Region selection).
4. Click 'calcuclate'
5. See error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/miniconda3/envs/jdaviz2.0.0/lib/python3.8/site-packages/ipyvue/VueTemplateWidget.py in _handle_event(self, _, content, buffers)
     58                 getattr(self, "vue_" + event)(data, buffers)
     59             else:
---> 60                 getattr(self, "vue_" + event)(data)
     61 
     62 

~/Desktop/jdat/jdaviz/jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py in vue_calculate_moment(self, *args)
    134         except ValueError:
    135             raise ValueError("Moment must be a positive integer")
--> 136         self.moment = analysis.moment(slab, order=n_moment)
    137 
    138         moment_ccd = CCDData(self.moment, unit=self.moment.unit)

~/miniconda3/envs/jdaviz2.0.0/lib/python3.8/site-packages/specutils/analysis/moment.py in moment(spectrum, regions, order, axis)
     40 
     41     """
---> 42     return computation_wrapper(_compute_moment, spectrum, regions,
     43                                order=order, axis=axis)
     44 

~/miniconda3/envs/jdaviz2.0.0/lib/python3.8/site-packages/specutils/analysis/utils.py in computation_wrapper(func, spectrum, region, **kwargs)
     17     # No region, therefore whole spectrum.
     18     if region is None:
---> 19         return func(spectrum, **kwargs)
     20 
     21     # Single region

~/miniconda3/envs/jdaviz2.0.0/lib/python3.8/site-packages/specutils/analysis/moment.py in _compute_moment(spectrum, regions, order, axis)
     80             m1 = np.tile(m1, _shape).T
     81 
---> 82         return np.sum(flux * (dispersion - m1) ** order, axis=axis) / m0

~/miniconda3/envs/jdaviz2.0.0/lib/python3.8/site-packages/astropy/coordinates/spectral_quantity.py in __array_ufunc__(self, function, method, *inputs, **kwargs)
     86     def __array_ufunc__(self, function, method, *inputs, **kwargs):
     87         # We always return Quantity except in a few specific cases
---> 88         result = super().__array_ufunc__(function, method, *inputs, **kwargs)
     89         if ((function is np.multiply
     90             or function is np.true_divide and inputs[0] is self)

~/miniconda3/envs/jdaviz2.0.0/lib/python3.8/site-packages/astropy/units/quantity.py in __array_ufunc__(self, function, method, *inputs, **kwargs)
    484 
    485         # Call our superclass's __array_ufunc__
--> 486         result = super().__array_ufunc__(function, method, *arrays, **kwargs)
    487         # If unit is None, a plain array is expected (e.g., comparisons), which
    488         # means we're done.

ValueError: operands could not be broadcast together with shapes (33,39,2059) (39,33,2059) 

Expected behavior
Create Moment 2 map with no error.

Screenshots
If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

• OS: OSX 10.15.7
• Browser: firefox
• Jupyter: [run "jupyter --version"]

jupyter core : 4.6.3
jupyter-notebook : 6.0.3
qtconsole : 4.7.4
ipython : 7.14.0
ipykernel : 5.3.0
jupyter client : 6.1.11
jupyter lab : not installed
nbconvert : 6.0.7
ipywidgets : 7.5.1
nbformat : 5.0.6
traitlets : 5.0.5

Package versions (please complete the following information):

macOS-10.15.7-x86_64-i386-64bit
Python 3.8.2 | packaged by conda-forge | (default, Apr 24 2020, 07:56:27)
[Clang 9.0.1 ]
Numpy 1.18.5
astropy 4.3.1
specutils 1.4.1
spectral-cube 0.5.0
pyyaml 5.4.1
click 7.1.2
asteval 0.9.23
idna 2.10
traitlets 5.0.5
bqplot 0.12.30
bqplot-image-gl 1.4.4
glue-core 1.2.2
glue-jupyter 0.10.1
glue-astronomy 0.3.2
echo 0.5
ipyvue 1.5.0
ipyvuetify 1.8.1
ipysplitpanes 0.1.0
ipygoldenlayout 0.3.0
voila 0.2.6
vispy 0.6.6
Jdaviz 2.0.0

Additional context (e.g. data files)

Run on this PR: #547

🐱

[pllim]

Notes for dev who will work on this:

np.sum(flux * (dispersion - m1) ** order, axis=axis) / m0

What in that equation has shape (33,39,2059) and what has (39,33,2059)? Which shape is the correct one? Why the other one is incorrect?

[pllim]

I am not sure exactly which file Patrick was using. I get a different traceback from a jw00636-o003_t002_nirspec_prism-clear_s3d.fits I have in my stash but I don't remember how I got it. So this might or might not be the same problem...

.../specutils/manipulation/extract_spectral_region.py in _to_edge_pixel(subregion, spectrum)
     58                 if hasattr(spectrum.wcs, "spectral"):
---> 59                     left_index = spectrum.wcs.spectral.world_to_pixel(left_reg_in_spec_unit)
     60                 else:

.../astropy/wcs/wcsapi/high_level_api.py in world_to_pixel(self, *world_objects)
    308 
--> 309         world_values = high_level_objects_to_values(*world_objects, low_level_wcs=self.low_level_wcs)
    310 

.../astropy/wcs/wcsapi/high_level_api.py in high_level_objects_to_values(low_level_wcs, *world_objects)
    232         if callable(attr):
--> 233             world.append(attr(objects[key]))
    234         else:

.../astropy/wcs/wcsapi/fitswcs.py in value_from_spectralcoord(spectralcoord)
    562                     else:
--> 563                         return spectralcoord.with_observer_stationary_relative_to(observer).to_value(**kwargs)
    564 

.../astropy/units/decorators.py in wrapper(*func_args, **func_kwargs)
    303             with add_enabled_equivalencies(self.equivalencies):
--> 304                 return_ = wrapped_function(*func_args, **func_kwargs)
    305 

.../astropy/coordinates/spectral_coordinate.py in with_observer_stationary_relative_to(self, frame, velocity, preserve_observer_frame)
    582         if self.observer is None or self.target is None:
--> 583             raise ValueError("This method can only be used if both observer "
    584                              "and target are defined on the SpectralCoord.")

ValueError: This method can only be used if both observer and target are defined on the SpectralCoord.

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
.../ipyvue/VueTemplateWidget.py in _handle_event(self, _, content, buffers)
     58                 getattr(self, "vue_" + event)(data, buffers)
     59             else:
---> 60                 getattr(self, "vue_" + event)(data)
     61 
     62 

.../jdaviz/configs/cubeviz/plugins/moment_maps/moment_maps.py in vue_calculate_moment(self, *args)
    130         spec_min = float(self.spectral_min) * u.Unit(self.spectral_unit)
    131         spec_max = float(self.spectral_max) * u.Unit(self.spectral_unit)
--> 132         slab = manipulation.spectral_slab(cube, spec_min, spec_max)
    133 
    134         # Calculate the moment and convert to CCDData to add to the viewers

.../specutils/manipulation/extract_spectral_region.py in spectral_slab(spectrum, lower, upper)
    241     region = SpectralRegion(lower, upper)
    242 
--> 243     return extract_region(spectrum, region)
    244 
    245 

.../specutils/manipulation/extract_spectral_region.py in extract_region(spectrum, region, return_single_spectrum)
    145     extracted_spectrum = []
    146     for subregion in region._subregions:
--> 147         left_index, right_index = _to_edge_pixel(subregion, spectrum)
    148 
    149         # If both indices are out of bounds then return an empty spectrum

.../specutils/manipulation/extract_spectral_region.py in _to_edge_pixel(subregion, spectrum)
     62                 left_index = int(np.ceil(left_index))
     63             except Exception as e:
---> 64                 raise ValueError(
     65                     "Lower value, {}, could not convert using spectrum's WCS "
     66                     "{}. Exception: {}".format(

ValueError: Lower value, 5.77931382762185e-07 m, could not convert using spectrum's WCS WCS Keywords

Number of WCS axes: 3
CTYPE : 'WAVE'  'DEC--TAN'  'RA---TAN'  
CRVAL : 5.77931382762185e-07  -2.2157499999999e-05  -3.9430000012079e-06  
CRPIX : 1.0  26.0  25.0  
PC1_1 PC1_2 PC1_3  : 1.0  0.0  0.0  
PC2_1 PC2_2 PC2_3  : 0.0  1.0  0.0  
PC3_1 PC3_2 PC3_3  : 0.0  0.0  1.0  
CDELT : 4.999999888241291e-09  2.77777781916989e-05  -2.77777781916989e-05  
NAXIS : 950  51  49.

Exception: This method can only be used if both observer and target are defined on the SpectralCoord.

[pllim]

To get moment 2 to work, I needed three different patches with the data that I have at hand. I can see if any of the patches is actually unnecessary if someone can give me a properly sanctioned NIRSpec IFU cube file.

1. EXP: Ignore obsgeo in JWST WCS for moment 2 calculation #1011
2. Allow descending spectral axes and more flexible SpectralRegion bounds astropy/specutils#911
3. BUG: Transpose m1 before transposing again astropy/specutils#912

[pllim]

Update: @kecnry pointed me to https://www.stsci.edu/jwst/science-planning/proposal-planning-toolbox/simulated-data (thanks, Kyle!). Using NRS00001-S1S2-F290LP-G395H-04_1_491_SE_2021-10-29T12h15m00_noflat_s3d.fits from there, I only need this patch for Moment 2 calculation to work:

• BUG: Transpose m1 before transposing again astropy/specutils#912

[pllim]

Fixed upstream.

[pllim]

The quick fix was reverted in favor for a more proper fix after specutils release.