Rubin_sim calc_mag returnig NaN's?

Science Survey Strategy
#1

Hi,

I’m trying to take some reflectance spectra I have of solar system objects, turn them into SED’s and then calculate magnitudes for them in LSST filters (much like the example notebook provided in the repository for rubin_sim). Unfortunately, no matter what I seem to do regarding extending spectra to redder+bluer wavelengths allows for a calculation of a magnitude that isn’t nan. I have tried the exact code with an asteroid from the SMASSII dataset that is used in the example and I get accurate colour calculations, but not with my own spectra.

I have included a link to my github repository with my notebook if this helps. The spectra I am using is just a simple .dat file with two columns, one with reflectance and another wavelength (in um’s). Is there a major issue with what I am doing that could cause the supposed lack of filter-SED wavelength overlap? Thanks in advance for any help anyone can provide!

Joe

#2

Hi Joe -

I’m not particularly familiar with rubin_sim, but I do note that the docs at this link state that the wavelength/flux values need to be on the same grid as the bandpasses that are used. Maybe you have to do some resampling?

Otherwise, I’m hoping @yoachim or @ljones can help you with more detailed responses.

-Jeff

#3

Hi Jeff,

Thanks for the response. I thought that might be the case, and I see I have forgotten to add that line to my notebook, but I have tried resampling the SED using resample_sed(wavelen_match=lsst[‘g’].wavelen) (for example) to no avail. I have assumed this method also resamples the flux values as well as the wavelengths, but perhaps this is not the case and that is my issue? Hopefully someone knows, but regardless thanks anyways!

~Joe

#4

Hi Joe,

Sorry - I also see that a previous community post linked to a notebook illustrating this issue and what to do about it, but the notebook was renamed … here’s an updated link so hopefully this calculating magnitudes notebook will be helpful.

I could investigate further by running through your notebook (is at least one of the SEDs you’re using available in the repo you linked? Or a dummy replacement example) tomorrow but I’ll try to give some pointers you can try –

  • resampling onto the same wavelength grid will happen automatically when you calculate a magnitude … it’s just nice to know this happens “under the hood” (and moreover, it’s good to know that the resampling occurs onto the wavelength grid of the filter).
  • the wavelength range of the SED has to exceed the wavelength range of the filter … otherwise it’s incomplete information when calculating the magnitude (even if the throughput of the filter is 0 at those wavelengths … the code isn’t currently smart enough to say “oh yeah, there’s no throughput here, so I don’t need Sed information”).
  • if the wavelength of the filter (min and max, as defined in .wavelen) exceeds in either the min or the max the wavelength range of the sed (min and max, as defined in .wavelen), then the result will be a nan. This indicates there was “missing information”.

Workarounds:

  • extend your SED, as you did in your notebook. It’s possible that you did not extend to beyond the wavelength ranges you’re getting from the filters? (you could just look at min/max values for each of your objects to check).
  • cut down the filters. If you know you’re not interested and do not have relevant information beyond 900 nm, for example, you could cut off the filters at that point (apply " .resample_bandpass(wavelen_min=, wavelen_max=<max_value>)" ) and then proceed as you were.

Some examples of this are given in that notebook linked above.

#5

Hi Lynne,

Thanks for the reply. I had previously seen that notebook and was looking through the suggested fixes. Whilst it may be strictly true that the redder and bluer bands won’t contain any useful information, I don’t necessarily want to disregard them altogether in doing any calculations.

As for extending the SED wavelength range, I had previously extended them to just a nm above and below the filter bandpasses and still gotten the nan result. I’ve just tried extending them 100nm either end below and above the filters and same result.

I’ve uploaded the spectrum that I’ve been using to the same github repo (first column wavelength in um, second column reflectance spectrum, so will need Solar spectrum multiplied back in). My gut feeling is I’m not doing the resampling correctly somehow, or that perhaps the flux values are too sparse and need smoothing? Hopefully if you have a chance to check over it you’ll be able to see where I’m making a mistake…

Thanks again for your help in this!

#6

Hi Joe! I forked your repo and sent you back a PR with an update to the notebook. Update CenCols notebook by rhiannonlynne · Pull Request #1 · astjoephysics/CentaurColours · GitHub

I’m not sure what happened when you previously extended the sed a bit further than the filter coverage … that should have worked. It was not the case in the version of the notebook that you had in the repo, though, so I did basically just ensure that the sed covered the full wavelength range (including making sure that numpy.arange was including the end points as I wanted), and then calc_mags works.
I did wonder if perhaps the sun_limited section was causing you some trouble, or at least, maybe causing some unexpected behavior under the hood. You do want to resample that sun spectra onto the pholus spectra wavelength, if you’re going to multiply them together … so I did that instead, so it’s possible that I also missed some error that was arising from that.

1 Like
#7

Hi Lynne,

I’ve just gone over your pull request - thanks for taking the time to sift through my code! I thought I had updated the wavelength range in that notebook before uploading…apparently not! Sorry about that!

Regardless, it seems that wasn’t my only issue anyways, and that resampling seems to have also been affecting this. I had previously been following the solar_system_object_colors example notebook and replacing the SMASSII sample spectrum with my own. Is there a reason I wouldn’t have to limit the solar wavelength range for my data similar to how it was done in your example notebook?

Thanks for all your help on this, I really appreciate it!

#8

Ah, good question. I think you’re referring to cell 17 from solar_system_object_colors

# Limit the wavelength range for the sun, just to reduce warning messages
sun_limited = phot_utils.Sed()
idx = np.where((sun.wavelen > 299) & (sun.wavelen < 1105))
sun_limited.set_sed(wavelen=sun.wavelen[idx], flambda=sun.flambda[idx])
sun_limited.resample_sed(wavelen_match=lsst['g'].wavelen)
wavelen_step = np.unique(np.diff(lsst['g'].wavelen)).min()
obs_smass = {}
for k in smass:
    smass[k].resample_sed(wavelen_match=lsst['g'].wavelen)
    obs_smass[k] = smass[k].multiply_sed(sun_limited, wavelen_step=wavelen_step)

whereas in the PR I sent to you I just did this between cell 5 and 6 in CenColors

# read in solar spectrum and MATCH wavelengths for pholus sed
sun = phot_utils.Sed()
sun.read_sed_flambda(os.path.join(get_data_dir(), 'movingObjects', 'kurucz_sun.gz'))

# don't need to sub-select sun idxes, this will do fine. 
sun.resample_sed(wavelen_match=phol_sed.wavelen)
# turn reflectance into sed and plot
# Looks like we have a bug, so set wavelen_step to side-step that .. sorry
phol_obs_sed = phol_sed.multiply_sed(sun, wavelen_step=0.1)

So … to be honest, I see warning messages in my notebook solar_system_object_colors after that cell … which are new and I didn’t notice when I reran it recently. (these errors could in fact be giving ME nans except that in a later cell I shortened the filter wavelength range which covered this problem!) and are the result of updates in Bandpass that changed how we read things in by default (compared to when I wrote the notebook in the first place). And if I wrote it now, I don’t think I’d do the subselect at all … I’d do this:

sun.resample_sed(wavelen_match=lsst['g'].wavelen)
obs_smass = {}
for k in smass:
    smass[k].resample_sed(wavelen_match=lsst['g'].wavelen)
    obs_smass[k] = smass[k].multiply_sed(sun_limited, wavelen_step=0.1)

which looks pretty similar to what’s in CenColors, right? (I should go update that notebook).

There is a slight difference that’s worth noting in terms of what is being used to define the resampled wavelengths. The version in CenColors is likely easier to read and more logical (resample the sun to match the reflectance spectra, then multiply them). You can then multiply the resulting spectra by whatever bandpass objects you have, even if they don’t have exactly the same wavelength range (as long as your SED completely overlaps all of them … J band colors, for example, if your SED goes into IR?)

The version in solar_system_object colors resampled both the sun and the reflectance spectra onto a third wavelength grid, the Bandpass and is a bit more specialized in that
you have to know that your Bandpass objects (if you’re calculating in multiple filters) will cover the same wavelength range and have the same wavelength gridding, and you need to know that these bandpasses are well sampled at a tight wavelength range. However, if you do know those things, it will be faster for use with multiple reflectance spectra, because you only resample the solar spectra once, only resample the reflectance spectra once per object, and since they all match the bandpass grid, there is no further resampling for calculating magnitudes.