Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

update phot_utils to be able to read throughput files #321

Merged
merged 6 commits into from
May 16, 2023

Conversation

yoachim
Copy link
Member

@yoachim yoachim commented May 10, 2023

The Bandpass class was setting wavelength parameters on instantiation that then prevented it from being able to read/set filters outside that range (e.g., 2MASS filters). Refactored so Bandpass only uses a wavelen attribute and no longer uses self.wavelen_min etc.

@rhiannonlynne
Copy link
Member

This isn't a bad idea, but I'd like to make sure it actually works in the context of a different-than-expected but not-uncommon use case where the input bandpass file only has a few data points and they're unevenly spread over wavelength space. The Sed would now be resampled only that uneven wavelength grid, rather than onto the currently enforced relatively tight grid (0.1nm for the wavelength spacing was chosen on purpose as a compromise between photometric accuracy and speed of calculation of magnitudes for our SED library). Then the photometry calculated from that would be inaccurate and potentially miss features in the spectrum. In the same case, the calculation of magnitudes is pretty basic and unsophisticated and assumes that the wavelength grid is uniform - without this in place, the resulting magnitudes will not be accurate.
(yes, I agree these are places that the unit tests are missing)

Did you consider just make the wavelength range (min/max) determined from the input data file instead to solve the problem of "input data files don't match the 300 to 1100 nm range expected from our standard throughput files"?

@yoachim
Copy link
Member Author

yoachim commented May 11, 2023

Maybe have Sed throw a warning if it is resampling below some resolution or to an uneven grid? Then Bandpass is more flexible, but users still know if they are entering the danger zone.

@rhiannonlynne
Copy link
Member

rhiannonlynne commented May 11, 2023 via email

@rhiannonlynne
Copy link
Member

rhiannonlynne commented May 11, 2023

Ok, great - this addresses the "sampling is too wide" but the calculation of magnitudes (in Sed) currently relies on bandpass being regularly gridded in wavelength (see "calc_flux" and how the total flux is calculated).
Can you update this to be a proper integral instead (to account for non-gridded wavelengths) or add a method to Sed that then imposes a regular gridding on bandpass? If Sed changes Bandpass, it may be wise to make a copy of the Bandpass object, since it presumably shouldn't (unexpectedly) change another object.

(one of the things that causes some overhead in Sed and Bandpass is that each time you call a magnitude or change something, the question arises "should I irretrievably change this or another item, or just make a temporary copy and act on the temporary copy" which is why you see "update_self" in various places. If you decide this shouldn't be a question any longer, it would be good to update this documentation as well and consider the implication of irretrievably changing the wavelength ranges of each of these when you may be calculating magnitudes for objects in different bandpasses that don't overlap. The use-cases for Sed that require this are laid out in the docstrings at the start of the class .. it may be fine to decide not to support some of those use-cases any longer).

@rhiannonlynne
Copy link
Member

OK, cool. So, can you please check that the notebooks to calculate m5 values in the syseng_throughputs repo give the same results as with the previous classes, that would be good. (we should see if we can get some CI going that includes those).
Clearly there are a lot of unit tests missing here that would have caught the problems that would come up with non-gridded wavelength files (you should also update the docstring in the classes that refers to uniformly spaced wavelength grids).

@rhiannonlynne
Copy link
Member

I'm still somewhat concerned that the warning isn't very strong and that people will get values that look reasonable but are reasonably inaccurate (see screenshot). And that telling them to just call resample_bandpass will result in them ending back up on 300-1150 nm wavelength range anyway (and maybe without realizing that's what's happening in that case either).
But it does look like it all hangs together, although I note that we get the "resample bandpass" warning about several throughput curves we have in our throughputs repository as examples (i.e. the SDSS bandpasses) now.

Screenshot 2023-05-15 at 6 35 17 PM

@yoachim
Copy link
Member Author

yoachim commented May 16, 2023

Just needed to update Bandpass.sb_tophi to do a proper integral. Now mags from arbitrary wavelength sampling should be fine.

@yoachim yoachim merged commit e995f5b into main May 16, 2023
5 checks passed
@rhiannonlynne rhiannonlynne deleted the tickets/OPSIM-1012 branch May 17, 2023 20:02
@rmjarvis
Copy link

rmjarvis commented Aug 28, 2023

Hi Peter, Lynne,

We just ran into the new wavelength sampling warning in our use of this in imSim.

/Users/Shared/cvmfs/sw.lsst.eu/darwin-x86_64/lsst_distrib/w_2023_34/conda/envs/lsst-scipipe-7.0.1-
ext/lib/python3.11/site-packages/rubin_sim/phot_utils/bandpass.py:67: UserWarning: Wavelength 
sampling of 16.1 nm is > 0.2 nm, this may not work well with a Sed object. Consider resampling with 
resample_bandpass method.

Does the last comment here:

Just needed to update Bandpass.sb_tophi to do a proper integral. Now mags from arbitrary wavelength sampling should be fine.

imply that it's actually fine that our bandpasses have fairly wide sampling? Or is this really a useful warning that we are probably doing something bad?

@cwwalter
Copy link
Member

Hi Peter and Lynne,

I have confirmed in our case these messages are being triggered by calling the skybrightness model:

image

@rhiannonlynne
Copy link
Member

rhiannonlynne commented Aug 29, 2023

So there are a couple of things -- yes, doing the 'proper integral' helps. But also, no because the sed object is resampled to whatever the bandpass object's wavelength array looks like. So if you have features in the Sed that are in a region where the Bandpass does not have good sampling (e.g. some of the older filter throughput curves just have points at the edges of the places where the filter actually was spec'd for throughput, as well as the 300/1200 nm), then you would lose that information.
(see

wavelen, fnu = self.resample_sed(wavelen, fnu, wavelen_match=bandpass.wavelen)
)
The Sed and Bandpass class have been heavily modified over time from where they started, which is I think generally the problem here .. although you do have a problem that you could miss features in one direction or the other so automatically just using the Sed sampling instead isn't quite the right answer either).

Thanks for the heads up about this issue.
The warning here is coming from reading the Canon filter curves from disk, these are only used in the calculation of skybrightness during twilight. If you're not predicting skybrightness during twilight time, it's completely irrelevant. If you are using twilight .. it's not entirely obvious to me if it's a problem or not, but we can fix it in a couple of different ways.

@cwwalter
Copy link
Member

Hi Lynne,

OK thanks! With your information, I was able to see what was going on. The twilight component of the sky model (which we initialize) is reading the three canon .csv files in:

skybrightness/Canon/

in the downloaded data directory. These are falling foul of your condition check and so is giving us the three warning errors. I'm not sure if we will ever simulate this situation in imSim for our tests or whether we will assume we have a calibrated system but I know @astrostubbs is interested in doing observations/flats during this period.

A few suggestions/requests:

  • It was quite hard to find these so putting the name of the throughput/filter file in the warning message would have been quite helpful and we could have more immediately had an idea of what was going on.

  • It would indeed be nice to fix this in someway (if is is actually a problem get a better file?). Otherwise, we are going to see these errors every time we start imSim and we would like to get rid of that.

@cwwalter
Copy link
Member

One other thought:

I assume for the sky we don't need the fine grained SED which is how it is being used here. This error is basically saying "Hey! don't use this filter curve you just read in with objects!" Maybe the place that the check was done could be moved so only if you try to use the curve in that way would the warning be triggered.

@rhiannonlynne
Copy link
Member

rhiannonlynne commented Aug 29, 2023

That's not a bad idea for moving the warning. You would, of course, get that warning still in this case .. because the Canon filter is used to calculate a magnitude with the solar spectrum in order to get twilight sky brightness estimates.
For now, I'm we're just going to quiet the warning by resampling the bandpass. This restores the older behavior, so should keep things consistent between this and older usage.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

4 participants