Filters example

Defining a wavelength array

We first define a wavelength range (in Angstrom) over which we want to define our filters. This wavelength array can be passed to the Filter at initialisation, and the Filter transmission curve will be defined on the passed wavelength array.

In practice, we would normally use the wavelength array from the Grid object to avoid conflicts between different wavelength arrays on different objects. Additionally, it is recommended to initialise a Filter / FilterCollection with a rest frame wavelength grid (i.e. the Grid wavelength array). Although nothing will break if you use an observer frame wavelength grid, the code expects rest frame wavelengths to ensure *.lam arrays agree across synthesizer objects. Any redshifting of Filter transmission curves is then handled automatically with minimal overhead.

One exception to this rule is when calculating photometry for a large number of objects all at a specific redshift. In that case, it is more efficient to use an observer frame wavelength grid on a Filter / FilterCollection, and simply be aware that Filter.lam differs from Grid.lam and SED.lam.

import numpy as np
from unyt import angstrom

from synthesizer.exceptions import (
    InconsistentWavelengths,
    WavelengthOutOfRange,
)
from synthesizer.filters import Filter, FilterCollection

lams = np.linspace(2000, 55000, 1000) * angstrom

Creating different types of Filter

To define a generic Filter we need to define some arbitray transmission curve that maps onto our wavelength array.

trans = np.zeros(lams.size)
trans[int(lams.size / 4) : int(lams.size / 2)] = 1
filt = Filter("generic/filter.1", transmission=trans, new_lam=lams)

To define a Top Hat Filter we define either the minimum and maximum wavelength of the transmission, or the effective wavelength of the transmission and its FWHM.

filt1 = Filter("top_hat/filter.1", lam_min=3000 * angstrom, lam_max=5500 * angstrom, new_lam=lams)
filt2 = Filter("top_hat/filter.2", lam_eff=7000 * angstrom, lam_fwhm=2000 * angstrom, new_lam=lams)

To create a Filter from the SVO database we simply pass the corresponding SVO filter identifier when intialising the Filter. You can get a summary of a Filter simply by printing it.

filt3 = Filter("JWST/NIRCam.F444W", new_lam=lams)
print(filt3)

Filter Code: JWST/NIRCam.F444W
Observatory: JWST
Instrument: NIRCam
Filter: F444W
Filter Type: SVO
SVO URL: http://svo2.cab.inta-csic.es/theory/fps/getdata.php?format=ascii&id=JWST/NIRCam.F444W
Wavelength Array: shape = (1000,), min = 2000.0 Å, max = 55000.0 Å
Frequency Array: shape = (1000,), min = 5.45e+13 Hz, max = 1.50e+15 Hz
Original Wavelength Array: shape = (1027,), min = 37131.8 Å, max = 50995.5 Å
Original Frequency Array: shape = (1027,), min = 5.88e+13 Hz, max = 8.07e+13 Hz
Transmission Curve: shape = (1000,), min = 0.0, max = 0.5146128183516849
Original Transmission Curve: shape = (1027,), min = 0.000105079, max = 0.514676

Creating a FilterCollection

We can intitialise a FilterCollection for each different type of Filter in a very similar way to what is shown above, only now the inputs are ensembles of filter properties, defined using a dictionary or a list.

You may have noticed in the above examples that each Filter was initialised with a string identifier, regardless of whether it was an SVO filter or not. These are the filter codes, used to identify each Filter object, and act as a key in dictionaries when required.

For a set of generic filters we provide a dictionary of transmission curves, where the key is the filter code for that Filter.

generics = {"generic_filter1": trans, "generic_filter2": trans}
generic_filters = FilterCollection(generic_dict=generics, new_lam=lams)

For a set of top hat filters we pass a dictionary containing the properties of each Filter, where the keys are the filter code of each Filter.

tophats = {
    "top_hat1": {"lam_eff": 25000 * angstrom, "lam_fwhm": 8000 * angstrom},
    "top_hat12": {"lam_min": 15000 * angstrom, "lam_max": 30000 * angstrom},
}
top_hat_filters = FilterCollection(tophat_dict=tophats, new_lam=lams)

For a set of SVO filters we pass a list of filter codes to be extracted from the database.

filter_codes = [f"JWST/NIRCam.{f}" for f in ["F070W", "F444W"]]
svo_filters = FilterCollection(filter_codes=filter_codes, new_lam=lams)

You can also pass a list of existing Filter objects.

filt_lst = [filt1, filt2, filt3]
filt_filters = FilterCollection(filters=filt_lst, new_lam=lams)

You are not tied to having a single type of Filter in a FilterCollection, they can be mixed and matched by passing in each one to the relevant argument. In the example below we combine SVO, tophat, generic and pre-defined filters in a single FilterCollection.

fs = [f"JWST/NIRCam.{f}" for f in ["F090W", "F250M"]]
tophats = {
    "U": {"lam_eff": 3650 * angstrom, "lam_fwhm": 660 * angstrom},
    "V": {"lam_eff": 5510 * angstrom, "lam_fwhm": 880 * angstrom},
    "J": {"lam_eff": 12200 * angstrom, "lam_fwhm": 2130 * angstrom},
}
generics = {"generic_filter3": trans}
filt_lst = [filt1, filt2]
combined_filters = FilterCollection(
    filter_codes=fs,
    tophat_dict=tophats,
    generic_dict=generics,
    filters=filt_lst,
    new_lam=lams,
)

One important thing to note when creating a FilterCollection is that all Filter objects in a FilterCollection will always be interpolated onto a common wavelength array. If a wavelength array has been passed to new_lam then all transmission curves will be interpolated on to it; otherwise, a common wavelength array is calculated using the wavelength arrays of all filters, with any gaps filled with values at the minimum mean resolution of all filters. Any overlapping areas will adopt the wavelength array of the higher resolution filter.

Working with a FilterCollection

A FilterCollection object behavea a bit like a Python list, and a bit like a Python dict; you can treat them as you would either of these data structures. Additionally, some helpful Python operators have been overloaded to make working with FilterCollection objects easier.

You can get their length, loop over them, compare them, extract specific keys, and add them; some of this functionality is demonstrated below.

print("We have %d filters" % len(combined_filters))

print("My Filters:")
for f in combined_filters:
    print(f.filter_code)

if combined_filters == combined_filters:
    print("This is the same filter collection!")
if combined_filters != generic_filters:
    print("These are not the same filter collection!")

filt = combined_filters["JWST/NIRCam.F090W"]
print(filt.filter_code)

new_filters = combined_filters + top_hat_filters
new_filters += generic_filters
print("My combined Filters:", new_filters.filter_codes)

We have 8 filters
My Filters:
JWST/NIRCam.F090W
JWST/NIRCam.F250M
U
V
J
generic_filter3
top_hat/filter.1
top_hat/filter.2
This is the same filter collection!
These are not the same filter collection!
JWST/NIRCam.F090W
My combined Filters: ['JWST/NIRCam.F090W', 'JWST/NIRCam.F250M', 'U', 'V', 'J', 'generic_filter3', 'top_hat/filter.1', 'top_hat/filter.2', 'top_hat1', 'top_hat12', 'generic_filter1', 'generic_filter2']

Note that adding individual Filter or FilterCollection objects to a FilterCollection will result in interpolation of the transmission curves onto the new FilterCollection objects wavelength array. In the event the wavelength arrays differ, the wavelengths of the FilterCollection being added to are maintained, and a warning is printed if this will truncate the transmission curve of the Filter being added. If the truncation results in 0 transmission then an error is raised.

filter_codes = [f"JWST/NIRCam.{f}" for f in ["F070W", "F444W"]]
new_filters2 = FilterCollection(filter_codes=filter_codes)
out_of_range_filt = Filter("JWST/MIRI.F2550W")

try:
    new_filters2 = new_filters2 + out_of_range_filt
except InconsistentWavelengths as e:
    print("InconsistentWavelengths:", e)

Calculated wavelength array:
min = 5.93e+03 Angstrom
max = 5.11e+04 Angstrom
FilterCollection.lam.size = 26891
InconsistentWavelengths: Interpolated transmission curve has no non-zero values. Consider removing this filter or extending the wavelength range.

Helper methods

We provide a number of helper methods and functions to make certain functionality quick and easy.

You can plot the transmission curves.

fig, ax = new_filters.plot_transmission_curves(show=True)

You provide a helper function to generate a UVJ top hat FilterCollection:

from synthesizer.filters import UVJ

lams = np.linspace(2500, 15000, 1000) * angstrom
fc = UVJ(new_lam=lams)
fig, ax = fc.plot_transmission_curves(show=True)

Each individual Filter has a number of methods to calculate useful filter properties.

print(f"The pivot wavelength of {filt.filter_code} is {filt.pivwv():.2e}.")
print(
    "The transmission at the pivot wavelength of "
    f"{filt.filter_code} is {filt.pivT():.2e}."
)
print(f"The mean wavelength of {filt.filter_code} is {filt.meanwv():.2e}.")
print(f"The bandwidth of {filt.filter_code} is {filt.bandw():.2e}.")
print(f"The FWHM of {filt.filter_code} is {filt.fwhm():.2e}.")
print(f"The transmission peak of {filt.filter_code} is {filt.Tpeak():.2e}.")
print(f"The rectangular width of {filt.filter_code} is {filt.rectw():.2e}.")
print(
    "The maximum wavelength with non-zero transmission of "
    f"{filt.filter_code} is {filt.max():.2e}."
)
print(
    "The minimum wavelength with non-zero transmission of "
    f"{filt.filter_code} is {filt.min():.2e}."
)
min_max = filt.mnmx()
print(
    "The minimum and maximum wavelength (regardless of transmission) of "
    f"{filt.filter_code} is {min_max[0]:.2e} and {min_max[1]}:.2e."
)

The pivot wavelength of JWST/NIRCam.F090W is 9.02e+03 Å.
The transmission at the pivot wavelength of JWST/NIRCam.F090W is 2.92e-01.
The mean wavelength of JWST/NIRCam.F090W is 8.98e+03 Å.
The bandwidth of JWST/NIRCam.F090W is 6.11e+02 Å.
The FWHM of JWST/NIRCam.F090W is 1.44e+03 Å.
The transmission peak of JWST/NIRCam.F090W is 3.19e-01.
The rectangular width of JWST/NIRCam.F090W is 1.94e+03.
The maximum wavelength with non-zero transmission of JWST/NIRCam.F090W is 1.02e+04 Å.
The minimum wavelength with non-zero transmission of JWST/NIRCam.F090W is 7.90e+03 Å.
The minimum and maximum wavelength (regardless of transmission) of JWST/NIRCam.F090W is 7.84e+03 Å and 10355.5 Å:.2e.

Saving and loading FilterCollections

You may encounter situations where you want to save a FilterCollection. This can be particularly useful where there is no network available (e.g. on certain isolated compute nodes), to avoid periods where the SVO database goes down, or when using a single FilterCollection across a number of MPI ranks where the multiple SVO database calls can cause failures.

To do this call the write_filters method, passing the path where you want to save the file. This will create a HDF5 file containing the attributes and datasets defining the FilterCollection.

# Write out the filter collection to a file in the synthesizer root directory
new_filters.write_filters(path="../../test_filter_collection.hdf5")

Note: this HDF5 file contains the version of synthesizer used to create it as an attribute (hdf.attrs["synthesizer_version"]). A warning will be printed if this differs from the version you are using, just in case future changes break the functionality of past FilterCollection objects.

To load a FilterCollection at a later date simply pass the path of the file into the FilterCollection at instantiation.

# Load a FilterCollection from a file
new_new_filters = FilterCollection(path="../../test_filter_collection.hdf5")

# Compare filter objects to show they are the same
new_filters == new_new_filters

True

Finding a Filter

If you want to know the Filter in a FilterCollection which will probe a particular wavelength, either in the rest–frame or observer–frame (given some redshift), you can use the find_filter method on a FilterCollection. This method will print the resulting filter code and return the Filter itself.

This function takes a method argument. There are 3 methods implemented for doing this: - “pivot”: Find the filter with the closest pivot wavelength to the desired wavelength (the default). - “mean”: Find the filter with the closest mean wavelength to the desired wavelength. - “transmission”: Find the filter with the maximum transmission at the desired wavelength.

# Define the wavelength we want to know the filter for (in the rest
# frame, with Angstrom units)
search_lam = 25000 * angstrom

# First lets search in the rest frame with the default pivot method
rest_frame_filter = new_filters.find_filter(search_lam)

# What if we used the other methods?
rest_frame_filter = new_filters.find_filter(search_lam, method="mean")
rest_frame_filter = new_filters.find_filter(search_lam, method="transmission")

# Now lets be more realistic, we have an observation redshift 7
# which filter can we use?
search_lam = 1000 * angstrom
obs_filter = new_filters.find_filter(search_lam, redshift=7, method="pivot")

Filter containing rest_frame_lam=2.50e+04 Angstrom: JWST/NIRCam.F250M
Filter containing rest_frame_lam=2.50e+04 Angstrom: JWST/NIRCam.F250M
Filter containing rest_frame_lam=2.50e+04 Angstrom: generic_filter3
Filter containing rest_frame_lam=1.00e+03 Angstrom (with observed wavelength=8.00e+03 Angstrom): top_hat/filter.2

Note that the method will produce an error if the transmission of all Filter objects in the FilterCollection is 0 at the desired wavelength. In certain circumstances "pivot" and "mean" can fail when the wavelength falls inside a wide band filter.

search_lam = 2000 * angstrom
try:
    obs_filter = new_filters.find_filter(
        search_lam, redshift=7, method="pivot"
    )
except WavelengthOutOfRange as e:
    print(e)

The wavelength (rest_frame_lam=2.00e+03 Å Angstrom, observed_lam=1.60e+04 Å Angstrom) has 0 transmission in the closest Filter (J). Try method='transmission'.

obs_filter = new_filters.find_filter(
    search_lam, redshift=7, method="transmission"
)

Filter containing rest_frame_lam=2.00e+03 Angstrom (with observed wavelength=1.60e+04 Angstrom): generic_filter3

In most cases, however, if a filter isn’t found the "transmission" method will also fail.

try:
    obs_filter = new_filters.find_filter(
        search_lam, redshift=15, method="transmission"
    )
except WavelengthOutOfRange as e:
    print(e)

The wavelength (rest_frame_lam=2.00e+03 Å Angstrom, observed_lam=3.20e+04 Å Angstrom) does not fall in any Filters.