desihub / desisim

DESI simulations
BSD 3-Clause "New" or "Revised" License
16 stars 22 forks source link

new problem with how the spectral templates are being normalized #535

Closed moustakas closed 4 years ago

moustakas commented 4 years ago

Clearly we need more unit tests (either here or in speclite).

In hunting down #534 I'm finding that speclite (specifically, speclite.filters._get_table) is not correctly normalizing the templates generated by the various desisim.templates methods. I think the issue is the newer versions astropy, which used to work on input spectra that did not have astropy.units attached to them. However, now they return errors along the lines of

*** ValueError: Values units None not convertible to erg / (Angstrom cm2 s).

Subsequently, normmaggies (e.g., here) returns all zeros (since all the table rows returned by speclite are masked).

This issue may be related to #507, too, although I'm not certain yet.

moustakas commented 4 years ago

MWE--

from desisim.templates import ELG
ELG().make_templates()
INFO:io.py:959:read_basis_templates: Reading /Users/ioannis/work/desi/spectro/templates/basis_templates/v3.2/elg_templates_v2.2.fits
/Users/ioannis/repos/desihub/desisim/py/desisim/templates.py:847: RuntimeWarning: divide by zero encountered in true_divide
  magnorm = 10**(-0.4*mag[ii]) / normmaggies
moustakas commented 4 years ago

Setting mask_invalid=False in the call to the get_ab_maggies returns the traceback below. Previously, mask_invalid=True was used to capture the case where the filter response function extended outside the wavelength range of the spectrum, but that's not the case here.

Oddly, using the example in the speclite documentation without specifying the units does not raise the same error, so I'm somewhat confused-- https://github.com/dkirkby/speclite/blob/master/speclite/filters.py#L7-L46

from desisim.templates import ELG
ELG().make_templates()
[snip]
~/repos/desihub/speclite/speclite/filters.py in get_ab_maggies(self, spectrum, wavelength, axis)
   1083         else:
   1084             # Allow interpolation since this is a convenience method.
-> 1085             convolution = self.convolve_with_array(
   1086                 wavelength, spectrum, photon_weighted=True,
   1087                 interpolate=True, axis=axis, units=default_flux_unit)

~/repos/desihub/speclite/speclite/filters.py in convolve_with_array(self, wavelength, values, photon_weighted, interpolate, axis, units, method)
   1039         convolution = FilterConvolution(
   1040             self, wavelength, photon_weighted, interpolate, units)
-> 1041         return convolution(values, axis, method)
   1042
   1043

~/repos/desihub/speclite/speclite/filters.py in __call__(self, values, axis, method, plot)
   1452                 'Must specify expected units for values with units.')
   1453         except astropy.units.UnitConversionError:
-> 1454             raise ValueError(
   1455                 'Values units {0} not convertible to {1}.'
   1456                 .format(values.unit, self.input_units))

ValueError: Values units None not convertible to erg / (Angstrom cm2 s).
moustakas commented 4 years ago

Fixed in #536.

Fortunately the fix was much easier than originally thought. Briefly, a flux array was inadvertently being passed to speclite as an astropy.table.Column, which astropy (internally to speclite) interpreted as something with inconsistent units (definitely something that used to work with earlier versions of astropy). So, closing this ticket for now--until the next problem!