ValueError: Wavelengths do not cover filter response 3330.0-10980.0 Angstrom.

[moustakas]

One should be able to synthesize photometry in a filter that is covered by the wavelengths of the spectrum, even if some filters are not covered. However, even in the latter case it would be nice if the "Wavelengths do not cover..." error was non-catastrophic and something functional (like np.nan) is returned in the bands that are either too blue or too red.

from astropy.io import fits
from speclite import filters
from desisim.templates import ELG

elg = ELG()
flux, wave, meta = elg.make_templates(nmodel=10)
wave.min(), wave.max()
  (3600.0, 10000.0)
decam = filters.load_filters('decam2014-r')
decam.get_ab_maggies(flux,wave)
[snip]
/usr/local/lib/python2.7/site-packages/speclite-0.4.dev338-py2.7.egg/speclite/filters.pyc in __init__(self, response, wavelength, photon_weighted, interpolate, units)
   1069                 .format(self.response._wavelength[0],
   1070                         self.response._wavelength[-1],
-> 1071                         default_wavelength_unit))
   1072 
   1073         # Find the smallest slice that covers the non-zero range of the

ValueError: Wavelengths do not cover filter response 3330.0-10990.0 Angstrom.

[dkirkby]

This was a bit trickier to deal with but I agree its a useful feature. The default is still to raise an error, but I added an new mask_invalid boolean option to get_ab_maggies() and get_ab_magnitudes() that will silently handle errors and mask the corresponding columns in the returned table:

sdss = speclite.filters.load_filters('sdss2010-*')
wlen = np.arange(5000, 10000) * u.Angstrom
flux = np.ones((3, 5000)) * u.erg / (u.cm**2 * u.s * u.Angstrom)
mags = sdss.get_ab_magnitudes(flux, wlen, mask_invalid=True)

This results in a table with 3 masked columns:

print(mags)
sdss2010-u sdss2010-g   sdss2010-r     sdss2010-i   sdss2010-z
---------- ---------- -------------- -------------- ----------
        --         -- -21.3589222106 -21.7817639913         --
        --         -- -21.3589222106 -21.7817639913         --
        --         -- -21.3589222106 -21.7817639913         --

Masked table values are a bit more flexible than using np.nan, but you can always convert them if you prefer:

print(mags.filled(np.nan))
sdss2010-u sdss2010-g   sdss2010-r     sdss2010-i   sdss2010-z
---------- ---------- -------------- -------------- ----------
       nan        nan -21.3589222106 -21.7817639913        nan
       nan        nan -21.3589222106 -21.7817639913        nan
       nan        nan -21.3589222106 -21.7817639913        nan

See here for details on masking in astropy tables.

I just pushed the changes to the master branch.

[moustakas]

This is a perfect solution, but is there any way to suppress the predictable warning message? For example the standard stars will never go red enough to cover all four WISE bandpasses:

/usr/local/lib/python2.7/site-packages/numpy/ma/core.py:4089: UserWarning: Warning: converting a masked element to nan.
  warnings.warn("Warning: converting a masked element to nan.")

[moustakas]

Unfortunately this is still causing some problems with the QSOs, which only cover the 3600-10000 A wavelength range. Normalizing the spectrum by the synthesized r-band magnitude throws the error

ValueError: Wavelengths do not cover decam2014-r response 3330.0-10990.0 Angstrom.

as before, but I can't use mask_invalid=True here because I need the normalization factor. The r-band filter, meanwhile, only really covers the range 5500-7200 A (and indeed even the g-band filter curve should be usable): http://speclite.readthedocs.org/en/latest/filters.html#decam-filters

Perhaps the wavelengths at which the filter profile response drops below some threshold (like 1%) can safely be ignored?

[dkirkby]

Regarding the numpy warning messages, the python warnings module lets you silence specific warnings but I'd like to understand why you are getting this warning in the first place. Can you try the following in a new python session:

>>> import numpy as np
>>> import astropy.units as u
>>> import speclite
>>> sdss = speclite.filters.load_filters('sdss2010-*')
>>> wlen = np.arange(5000, 10000) * u.Angstrom
>>> flux = np.ones((3, 5000)) * u.erg / (u.cm**2 * u.s * u.Angstrom)
>>> sdss.get_ab_magnitudes(flux, wlen, mask_invalid=True)
<Table masked=True length=3>
sdss2010-u sdss2010-g   sdss2010-r     sdss2010-i   sdss2010-z
 float64    float64      float64        float64      float64
---------- ---------- -------------- -------------- ----------
        --         -- -21.3589222106 -21.7817639913         --
        --         -- -21.3589222106 -21.7817639913         --
        --         -- -21.3589222106 -21.7817639913         --

I don't get any warnings when I do this. Does this give you a warning, or what are you doing differently in your code?

You can also exclude any filters that you never need, e.g. if you never use W3 and W4:

filters = speclite.filters.load_filters('decam2014-*', 'wise2010-W1', 'wise2010-W2')

[dkirkby]

For cases like your 3500-10000A QSO spectrum that does not cover decam2014-r, there are two complementary approaches:

• Provide a new set of decam filter files that are truncated to eliminate long tails with small throughputs.
• Extrapolate your QSO spectrum to cover the range of the filter.

The first approach is probably the most convenient for users, but a bit dangerous because it makes implicit assumptions about the spectrum (e.g., that no bright line falls into a truncated region) and involves using something other than the official DECam filter curves.

The second approach puts the burden on the user (you) but at least forces you to make an explicit choice about how to deal with this issue. The numpy pad function makes it easy to implement some obvious extrapolation schemes. For example, to extrapolate from 5000-10000A to 2000-12000A using the median value:

padded_wlen = np.arange(2000, 12000) * u.Angstrom
padded_flux = np.pad(flux, ((5000-2000, 12000-10000)), 'median')
sdss.get_ab_magnitudes(padded_flux, padded_wlen)

This should really be wrapped in something like:

def extrapolate(wlen, flux, new_min_wlen, new_max_wlen, method='median'):
    ...

which I could add to speclite.

[moustakas]

Can you add this to speclite without too much hassle?

[dkirkby]

I will implement something simple tomorrow so we can tag once you have tested it with desisim.

[dkirkby]

@moustakas Does this look like it does everything you need? http://speclite.readthedocs.org/en/latest/api/speclite.filters.FilterResponse.html#speclite.filters.FilterResponse.pad_spectrum

In the case of a filter group, this method would be called for each filter in turn (but I haven't implemented this yet).

[dkirkby]

@moustakas I just pushed a release candidate to the github master branch. The basic usage in your case is:

filters = speclite.filters.load_filters(...)
padded_flux, padded_wave = filters.pad_spectrum(flux, wave)
maggies = filters.get_ab_maggies(padded_flux, padded_wave)

Of course, this is generally a dangerous thing to do (and I tried to be sufficiently cautionary in the docs) but you know what you are doing :-)

I will release this version as 0.4 once you have had a chance to stress test the changes (assuming you don't break anything).