mfouesneau / pyphot

suite to deal with passband photometry
https://mfouesneau.github.io/pyphot/
MIT License
57 stars 18 forks source link

Astropy units? #18

Closed christopherlovell closed 4 years ago

christopherlovell commented 4 years ago

[an enhancement request] Are there plans to move to astropy units in the future?

mfouesneau commented 4 years ago

yes but man power limited, so no timescale

christopherlovell commented 4 years ago

I'd be happy to take a look if someone can provide some direction.

From a quick glance I can see two ezunit functions, unit and hasUnit, that require replacing in the following files:

There's also _drop_units in phot.py that needs redefining.

mfouesneau commented 4 years ago

There is more because of the internal conversion of the units. Some pieces assume a ezunit object with a given API.

christopherlovell commented 4 years ago

Astropy has methods for unit conversion, so as long as we can isolate any calls to the ezunit API we should be able to switch to the astropy API for this, right?

mfouesneau commented 4 years ago

Yes of course and adapt to the Astropy API. But that's time consuming

mfouesneau commented 4 years ago

Have a look into the astropy folder. I just made a first version (incomplete and not fully tested)

christopherlovell commented 4 years ago

OK great, I'll take a look. Is it worth putting all this on an 'astropy' development branch for the moment?

mfouesneau commented 4 years ago

Not sure that a branch is what needs to be done. I need to keep the old units too so both should co-exist in some ways.

christopherlovell commented 4 years ago

sure, I agree that keeping the old units is worthwhile, but any breaking changes will make it in to the public release if they're on the main branch is all. A separate branch keeps the astropy changes isolated until we're happy to merge

mfouesneau commented 4 years ago

The astropy code (in the astropy folder) seems to work. I added the notebooks in the examples folder.

Lick API not completely independent from the previous version yet, but very close.

mfouesneau commented 4 years ago

from pyphot import astropy as pyphot from pyphot.astropy import Vega, Sun ... should now be equivalent to the other version using Astropy units

christopherlovell commented 4 years ago

Testing this now.

If you wish to define a filter in fnu rather than flam, how do you access the @set_method_default_units decorator arguments at runtime?

For example, for the get_flux method, defining flux_unit='fnu' as an argument in either the UnitFilter instantiation, or directly in the call to get_flux, doesn't seem to be recognised.

mfouesneau commented 4 years ago

I've never done that sort of definition. 'fnu' is not a unit that can be easily converted. Maybe making a decorator that does fnu to flambda conversion?

If you were to give an example maybe I could be more constructive, sorry

christopherlovell commented 4 years ago
import numpy as np
from astropy import units as u
from pyphot.astropy import UnitFilter

wl = np.arange(845,855,0.1) * u.micron
flux = np.random.rand(len(wl)) * (u.erg * u.s**-1 * u.Hz**-1 * u.cm**-2)

filt_wl = [845,846,850,854,855] * u.micron
filt_trans = [0.,1.,1.,1.,0.]
tophat = UnitFilter(filt_wl, filt_trans, name='tophat', dtype='energy', unit='micron')
flux_tophat = tophat.get_flux(wl, flux)

So the above example fails with the following:

astropy.units.core.UnitConversionError: 'erg / (cm2 Hz s)' (spectral flux density) and 'flam' (spectral flux density wav) are not convertible

I thought that the previous version of pyphot handled this conversion for you, i.e. if you define the filter in wavelength it will convert to frequency space if the flux is defined in f_nu. But perhaps not?

christopherlovell commented 4 years ago

In fact, even using consistent frequency units for the filter and the spectrum still fails, as 'flam' and the output units is hard coded in the decorator

@set_method_default_units('AA', 'flam',
                          output_unit='erg*s**-1*cm**-2*AA**-1')
def get_flux(self, slamb, sflux, axis=-1):
....
mfouesneau commented 4 years ago

Wait that's normal because Astropy does not know what to do here: pivot wavelength 𝜆𝑝 is necessary to convert between 𝑓𝜈 and 𝑓𝜆 𝑓𝜈 = (𝜆/c)^2𝑓𝜆

I don't think this is a code issue, but a conceptual conversion. I do not have a piece of code converting these quantities tho. Do you agree?

christopherlovell commented 4 years ago

yep, you're right, my misunderstanding.

I'll close this issue for now as it all seems to be working

mfouesneau commented 4 years ago

A conversion method could be implemented to the Filter class if that helps?

christopherlovell commented 4 years ago

Something like that would be really nice, but it's not urgent. It would also require updating the documentation a bit to make it clear what it does, and avoid misuse.