Hi,
I’ve been trying to use rubin_sim.photUtils.Sed.calcMag on real-life spectra and I discovered behaviour that seems strange to me. If I try to calculate mags on a SED that does not completely overlap in wavelengths with the LSST wavelength range, I get a warning rubin_sim/photUtils/Sed.py:870: UserWarning: There is an area of non-overlap between desired wavelength range (which is reasonable), and then the produced magnitudes are all NaN. The easiest way to reproduce this is to delete the first rows (up to the #nm=301.393) in the ‘kurucz_sun’ SED, included in the tutorial.
Is there a bug or it is intended behaviour? Or I’m simply doing something wrong?
Hi @RShA , that does sound like intended behaviour to me, but I’ll call @ljones into this discussion.
Thanks for the tag Melissa.
Yes - this is intended behavior. The options, from our design point of view if the SED does not overlap the full wavelength range, were:
- fill in zeros
- extend the last value
- add Nans for the wavelengths which don’t match
Out of these, the first two will silently change the resulting magnitude and someone may not notice the error. The last option (which is what we do) results in a nan for the magnitude, but is clear.
I added some examples of actions a user can take to get around the problem to the end of the demo notebook - Calculating Magnitudes in GitHub lsst/rubin_sim_notebooks.
(notebook was moved – now at https://github.com/lsst/rubin_sim_notebooks/blob/main/photometry/calculating_magnitudes.ipynb)
@ljones hi, thanks for the answer!
I understand the logic behind this behaviour, but I am still not sure that the implementation is entirely correct. The thing is, in the given example with short_sun SED not only z_mag will be nan, but all other magnitudes as well. This doesn’t seem to be reasonable behaviour: if I lack SED data on the edge of the LSST range, I believe, I still should be able to generate mags in all the bands for which SED is present.
If I understand correctly, the current behaviour is caused by two lines in calcFlux function:
dlambda = wavelen[1] - wavelen[0]
flux = (fnu * bandpass.phi).sum() * dlambda
first, dlambda here is going to be nan, if SED is incomplete in u band, and second, sum() for the array with nans will produce nan. Since it sums fnu as a whole, not only the part of fnu which belongs to the given band, any nans in the fnu will result in nans for all mags.
For my purposes I modified the code, changing the abovementioned lines to
dlambda = bandpass.wavelen_step
flux = numpy.nansum(fnu * bandpass.phi) * dlambda
but the right move here, probably, would be writing a solution that would produce nans for the bands for which SED is absent and correct mags for the bands that have full overlap with SED.
So here is the thing: the filter bandpasses DO cover the wavelength range 300-1150 nm, in all bandpasses. Yes, the sensitivity is not high throughout most of this range, but the classes do not assume that. (which is why saying “I should be able to generate mags in all the bands for which the SED is present” – is true. But you are assuming that the wavelength range for each bandpass is different, which is not the case, for reasons which are also spelled out in the notebook in the previous ‘faster mag calculation’ stages).
Using np.nansum is equivalent to assuming fnu = 0 in the non-overlap region, which is an option we chose not to take. You could also test that phi is 0 in the range of non-overlap, which we also do not do (but as you point out, we could – it’s additional overhead that didn’t really match the original purpose of the code, which was primarily for calculating magnitudes in all bandpasses, for mock catalog generation).
As in the examples in the notebook, if you know your bandpass sensitivity is zero outside the wavelength range of your input SED, you can set the wavelength range to this abbreviated version for the filter / Bandpass object. Then the NaN will not occur.
I suppose a bigger point here is that the use of the photUtils package is maybe primarily NOT the original use case. Would an acceptable/useful alternative from your point of view be to simply add a method to the Bandpass class called something like “trim_bandpass” which would reduce the wavelength range for the bandpass class itself to the minimum range required? Then, after reading in the bandpasses, this could optionally be called to make the bandpass wavelength range smaller and produce the desired effect of magnitudes being reported in bands which overlap the actual transmission range, and NaNs when there is relevant missing information.
It would have the effect of slowing down “faster magnitude calculation” but should have no impact on the ‘standard’ use. I would make the default behavior be to trim, but allow the user to turn this off (thus setting up things for simpler use as standard, but leaving the option for advanced/faster use).
Hi again,
ok, thanks for your explanation, originally I missed the point about why we need to set all bandpasses wavelength range to the same values.
I think creating a method for trimming the bandpass is a very good idea. It’s definitely going to be useful, considering that there are many source SEDs that can be used for generating mock catalogues and not all of them have the same coverage as LSST filters. I am not even sure even if it has to be the default behaviour. Maybe it’s easier and better to leave the current version as the default one but add the suggestion to use ‘trim_bandpass’ to the warning about non-overlapping areas of the wavelength scale. I am honestly not sure which one is going to be the more common case: the need to generate mags for large catalogues with good wavelength range or the need to generate mags for incomplete SEDs.