Improve spaxel unit handling for cubes

[rosteen]

#3221 broke MaNGA data cube loading, because we didn't account for spaxel being in the flux units anywhere in our code. I added that to the list of recognized solid angle units in check_if_unit_is_per_solid_angle, but this probably needs a little more thought about whether to expose spaxel in the unit conversion plugin along with pix^2, or whether to always convert spaxel to pix^2 on data load and just have pix^2 in the unit conversion plugin. @havok2063 @camipacifici any thoughts on those two options?

Note that I milestoned this as 4.1, loading still works in 4.0.x because #3221 was milestoned to 4.1.

Still investigating why the current MaNGA loading test didn't catch this.

[rosteen]

It seems like MaNGA files might have changed format somewhat in the last few years. The test file we have is from 2020 and loads fine on main, where spaxel doesn't show up in the units when hovering over the cube in Cubeviz. But the file I was testing with today that @camipacifici linked me to fails on main and shows spaxel in the unit when hovering. However, loading the old cube with this PR also shows spaxel in the units on hover, so I'm not sure where that difference is happening. Quite confused tbh, @havok2063 any insight?

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 88.80%. Comparing base (9a90da0) to head (8e50ce7).
Report is 1 commits behind head on main.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #3307   +/-   ##
=======================================
  Coverage   88.80%   88.80%           
=======================================
  Files         125      125           
  Lines       19118    19118           
=======================================
  Hits        16978    16978           
  Misses       2140     2140           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[rosteen]

Btw, I don't think we were really handling this correctly before #3221 either, that just exposed the problem in a way that errors out.

[havok2063]

@rosteen Regarding spaxel vs pix2, I'm in favor of spaxel personally. I do feel it's becoming a more common unit of expression, however it's not a formal astropy-supported unit. If it's easy to keep in, I'd say use it, but I also don't necessarily want Jdaviz to adopt a non-standard unit.

There haven't been any changes to the MaNGA files since their last data release in 2021, reduced in 2020. You can check the header keyword VERSDRP2. If it says v3_1_1 then it's the latest file. Otherwise it's an older version.

The flux units have always been set to the following in the headers, for both new and older files.

BSCALE  =              1.00000 / Intensity unit scaling
BZERO   =              0.00000 / Intensity zeropoint
BUNIT   = '1E-17 erg/s/cm^2/Angstrom/spaxel' / Specific intensity (per spaxel)

It's impossible to change the MaNGA data itself now so a fix has to be either in jdaviz or specutils. Do the manga data loaders for specutils need updating? Maybe it's an issue between some of the latest jdaviz and specutils changes on units, parsing, and loaders?

[rosteen]

It looks like we are using an older file version in our test:

VERSDRP2= 'trunk 9661'         / MaNGA DRP version (2d processing)
VERSDRP3= 'trunk 63586'        / MaNGA DRP Version (3d processing)
VERSPLDS= 'v2_52   '           / Platedesign Version
VERSFLAT= 'trunk 161238'       / Specflat Version
VERSCORE= 'trunk 9678'         / MaNGAcore Version
VERSUTIL= 'trunk 63620'        / Version of idlutils
[...]
BUNIT   = '1E-17 erg/s/cm^2/Angstrom/spaxel' / Flux units are per spaxel

In regards to spaxel being an astropy-supported unit, while u.spaxel doesn't work, u.Unit('spaxel') does (although it doesn't know the physical type).

I think the fix belongs here in Jdaviz, this PR already lets us load the data properly, it was just a matter of whether I add a little more spaxel support throughout the app for after it's loaded.

Edit: from the edit I made in the test, looks like the old file was version 2_1_2

[rosteen]

Gah, I see now what you were saying that spaxel isn't fully supported by astropy, @havok2063 . Making a Quantity out of a unit string with spaxel in it silently drops the spaxel, which messes up how our Unit Conversion plugin gets the angular unit. I think we may need to convert spaxel to pix^2 on load to get this to work correctly, unfortunately. I guess I could also open something upstream in astropy to get spaxel working correctly.

[havok2063]

Ahh then in that case, if it's easier to switch to pix2, then I think that's ok. Do we even need an angular unit? How did Jdaviz handle manga data before before, say last year or so? v2_1_2 is quite old indeed. Maybe if no one was using newer manga, then they never ran into this issue? Although I've been using it with later data.

[rosteen]

We've had a whole push to handle surface brightness vs flux better in the last year (amidst the bigger Unit Conversion work), which included adding /pix^2 to units of flux that were really representing a surface brightness.

[james-trayford]

Thanks for looking into this @rosteen - I tried your changes in a patch to let me read in MaNGA cubes - this works for me now (thanks!) though I came across a unit error that may be related when trying to fit a model to the spectrum (when I trigger ADD COMPNENT) I get:

UnitConversionError: 'erg / (Angstrom s cm2)' (power density/spectral flux density wav) and '' (dimensionless) are not convertible

Details

---------------------------------------------------------------------------
UnitConversionError                       Traceback (most recent call last)
File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/ipyvue/VueTemplateWidget.py:60, in Events._handle_event(self, _, content, buffers)
     58     getattr(self, "vue_" + event)(data, buffers)
     59 else:
---> 60     getattr(self, "vue_" + event)(data)

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/jdaviz/configs/default/plugins/model_fitting/model_fitting.py:763, in ModelFitting.vue_add_model(self, event)
    762 def vue_add_model(self, event):
--> 763     self.create_model_component()

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/jdaviz/configs/default/plugins/model_fitting/model_fitting.py:462, in ModelFitting.create_model_component(self, model_component, model_component_label, poly_order)
    459 if comp_label in [cm['id'] for cm in self.component_models]:
    460     raise ValueError(f"model component label '{comp_label}' already in use")
--> 462 new_model = self._initialize_model_component(model_comp, comp_label, poly_order=poly_order)
    463 self.component_models = self.component_models + [new_model]
    464 # update the default label (likely adding the suffix)

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/jdaviz/configs/default/plugins/model_fitting/model_fitting.py:541, in ModelFitting._initialize_model_component(self, model_comp, comp_label, poly_order)
    538 # equivs for spectral density and flux<>flux/pix2. revisit
    539 # when generalizing plugin UC equivs.
    540 equivs = _eqv_flux_to_sb_pixel() + u.spectral_density(init_x)
--> 541 init_y = init_y.to(self._units['y'], equivs)
    543 initialized_model = initialize(
    544     MODELS[model_comp](name=comp_label,
    545                        **initial_values,
    546                        **new_model.get("model_kwargs", {})),
    547     init_x, init_y)
    549 # need to loop over parameters again as the initializer may have overridden
    550 # the original default value. However, if we toggled cube_fit, we may need to override

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/astropy/units/quantity.py:943, in Quantity.to(self, unit, equivalencies, copy)
    939 unit = Unit(unit)
    940 if copy:
    941     # Avoid using to_value to ensure that we make a copy. We also
    942     # don't want to slow down this method (esp. the scalar case).
--> 943     value = self._to_value(unit, equivalencies)
    944 else:
    945     # to_value only copies if necessary
    946     value = self.to_value(unit, equivalencies)

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/astropy/units/quantity.py:896, in Quantity._to_value(self, unit, equivalencies)
    893     equivalencies = self._equivalencies
    894 if not self.dtype.names or isinstance(self.unit, StructuredUnit):
    895     # Standard path, let unit to do work.
--> 896     return self.unit.to(
    897         unit, self.view(np.ndarray), equivalencies=equivalencies
    898     )
    900 else:
    901     # The .to() method of a simple unit cannot convert a structured
    902     # dtype, so we work around it, by recursing.
    903     # TODO: deprecate this?
    904     # Convert simple to Structured on initialization?
    905     result = np.empty_like(self.view(np.ndarray))

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/astropy/units/core.py:1196, in UnitBase.to(self, other, value, equivalencies)
   1194     return UNITY
   1195 else:
-> 1196     return self._get_converter(Unit(other), equivalencies)(value)

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/astropy/units/core.py:1125, in UnitBase._get_converter(self, other, equivalencies)
   1122             else:
   1123                 return lambda v: b(converter(v))
-> 1125 raise exc

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/astropy/units/core.py:1108, in UnitBase._get_converter(self, other, equivalencies)
   1106 # if that doesn't work, maybe we can do it with equivalencies?
   1107 try:
-> 1108     return self._apply_equivalencies(
   1109         self, other, self._normalize_equivalencies(equivalencies)
   1110     )
   1111 except UnitsError as exc:
   1112     # Last hope: maybe other knows how to do it?
   1113     # We assume the equivalencies have the unit itself as first item.
   1114     # TODO: maybe better for other to have a `_back_converter` method?
   1115     if hasattr(other, "equivalencies"):

File ~/Documents/Code/strauss_dev/venv/lib/python3.11/site-packages/astropy/units/core.py:1086, in UnitBase._apply_equivalencies(self, unit, other, equivalencies)
   1083 unit_str = get_err_str(unit)
   1084 other_str = get_err_str(other)
-> 1086 raise UnitConversionError(f"{unit_str} and {other_str} are not convertible")

UnitConversionError: 'erg / (Angstrom s cm2)' (power density/spectral flux density wav) and '' (dimensionless) are not convertible

[rosteen]

Thanks for looking into this @rosteen - I tried your changes in a patch to let me read in MaNGA cubes - this works for me now (thanks!) though I came across a unit error that may be related when trying to fit a model to the spectrum (when I trigger ADD COMPNENT) I get:

UnitConversionError: 'erg / (Angstrom s cm2)' (power density/spectral flux density wav) and '' (dimensionless) are not convertible

Details

Thanks for testing this, I was aware of a couple things that weren't working but not this specifically - I'm working on refactoring a bit in a way that I think will fix this. I'll take this PR out of draft once I think I have it all working.

[rosteen]

Ok, this PR fixes the actual data loading of the MaNGA cubes, but exposed two related bugs/problems that I'm going to make a follow up for:

1. Spectral extraction isn't considering the mask cube if there was one loaded in with the flux cube. This can lead to incorrect results for the extracted spectrum, e.g., if the masked results were set to 0 and we do a mean extraction.

2. Model Fitting is failing for very small flux values, since we removed the factor of e.g. 1E-17 from the unit and put it on the flux values themselves rather than keeping it in the unit. I don't think the fitter will converge for such small values. So we probably need to take out that factor for model fitting and then re-apply it to the results.

[rosteen]

I think the docs build failure is related to Astropy release, not this PR. Not 100% sure though.

[kecnry]

should we just always replace infs (maybe with nans instead of zeros) regardless of the uncertainty type? Or should that be handled in individual plugins so that we do not modify original data?

[rosteen]

Not sure, let's discuss offline.

[pllim]

This is already in specutils parser (but not exposed). Should this be moved upstream eventually?

[rosteen]

I think so, for now I just made our spaxel definition match the specutils one so they would actually divide properly. And whoops, looks like I have a double spaxel = there.

[pllim]

This is probably gonna conflict with #3228 . Should @cshanahan1 review?

[pllim]

jdaviz.utils.flux_conversion:22: WARNING: py:obj reference target not found: astropy.units.equivalencies

Use intersphinx

:ref:`astropy:unit_equivalencies`

[rosteen]

• Spectral extraction isn't considering the mask cube if there was one loaded in with the flux cube. This can lead to incorrect results for the extracted spectrum, e.g., if the masked results were set to 0 and we do a mean extraction.

• Model Fitting is failing for very small flux values, since we removed the factor of e.g. 1E-17 from the unit and put it on the flux values themselves rather than keeping it in the unit. I don't think the fitter will converge for such small values. So we probably need to take out that factor for model fitting and then re-apply it to the results.

Upon further investigation, 2 is just a manifestation of 1, not due to the values being too small. So only one follow up issue!

[gibsongreen]

Everything looks good to me here!

I did notice one bug that I confirmed was not created by this PR. When a manga (or similarly counts) cube is loaded, if you translate, the the values do not update according to the scale factor. This PR made this more visible when changing flux units of the manga cube after translating (resulting in values not updating and zoom/other features no longer functioning until reverting to the original spectrum-y value).

[cshanahan1]

in the aperture photometry plugin, there is some logic specifically referring to PIX2 that should also probably be applied to spaxel. the pixel area is fixed to 1.0 for pix2 cubes, which I think should also be the case for spaxel? (the PIX2 logic in lines ~210-220 in aper_phot_simple.py are fixed in #3228).

There is also some logic for PIX2 done in coords_info when filtering out certain data types so their units aren't converted for mouseover, and spaxel should be added to that too, but I would do this off #3228 since there are changes to this part of the code there as well.

[cshanahan1]

we should probably define the PIX2 unit like this right? should i make a ticket?

[cshanahan1]

I'm not sure what the best order of operations here is since these changes will overlap with some of the changes made in #3228. It might make more sense for that to go in first? That way, you can add 'spaxel' to the parameterized tests on that branch that test unit conversion in plugins. Additionally, the same logic added to utils.flux_conversion here will need to be added to unit_conversion_utils.flux_conversion_general on my branch so spaxel will work in plugins with unit conversions. (I hope these two functions can be unified at some point, but for now utils.flux_conversion_general handles the spectral axis and unit_conversion_utils.flux_conversion_general handles all other flux unit conversions).

Another idea would be to merge this since it fixes loading for now, with the caveat that unit conversion is broken. Then merge my PR, then make a follow up ticket to add spaxel to the logic that my PR introduced.

[rosteen]

@cshanahan1 I'm happy to merge yours first and rebase this, that might be the easier way. I'll respond to the other comments shortly.

[rosteen]

@cshanahan1 just pointing out that this is where spaxel is getting converted to PIX2 in the parser - as discussed offline we shouldn't need spaxel handling elsewhere.