Images From Galaxy Particle Distributions¶
Synthesizer can create several different types of images and maps from particle data.
In the example below we demonstrate the main high-level imaging workflow for a galaxy stellar distribution taken from a CAMELS simulation galaxy. We will generate photometry, construct a PhotometricImager, and then call the galaxy image getters. Later sections briefly show how to work directly with the returned ImageCollection objects for PSF and noise manipulation.
[1]:
import time
import matplotlib.colors as cm
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
from astropy.cosmology import Planck18 as cosmo
from scipy import signal
from unyt import Hz, angstrom, arcsecond, erg, kpc, s
from synthesizer import TEST_DATA_DIR
from synthesizer.emission_models import BimodalPacmanEmission
from synthesizer.emission_models.attenuation import PowerLaw
from synthesizer.filters import FilterCollection as Filters
from synthesizer.grid import Grid
from synthesizer.kernel_functions import Kernel
from synthesizer.load_data.load_camels import load_CAMELS_IllustrisTNG
# Define the grid
grid_name = "test_grid"
grid = Grid(grid_name, new_lam=np.logspace(2, 5, 600) * angstrom)
# Create galaxy object
gal = load_CAMELS_IllustrisTNG(
TEST_DATA_DIR,
snap_name="camels_snap.hdf5",
group_name="camels_subhalo.hdf5",
physical=True,
)[1]
Getting the photometry¶
To make an image we need to generate spectra for each individual star particle. To do this we use the galaxy’s in built get_spectra method and pass a per_particle model. This will generate a spectrum for each particle in the galaxy.
We can then calculate the photometry on these spectra by defining a filter collection, and calculating the luminosities.
[2]:
# Get the stellar pacman model
model = BimodalPacmanEmission(
grid=grid,
tau_v_ism=1.0,
tau_v_birth=0.7,
dust_curve_ism=PowerLaw(slope=-1.3),
dust_curve_birth=PowerLaw(slope=-0.7),
fesc=0.1,
fesc_ly_alpha=0.9,
label="total",
per_particle=True,
)
print(model)
# And use it to generate the spectra
sed = gal.stars.get_spectra(model)
# Get the observed spectra
gal.get_observed_spectra(cosmo)
# Set up the filter collection for imaging
filter_codes = [
"JWST/NIRCam.F090W",
"JWST/NIRCam.F150W",
"JWST/NIRCam.F200W",
]
filters = Filters(filter_codes, new_lam=grid.lam)
gal.get_photo_lnu(filters)
|============================================ EmissionModel: total ===========================================|
|-------------------------------------------------------------------------------------------------------------|
| TOTAL (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_emergent, old_emergent |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_EMERGENT (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_attenuated_ism, young_escaped |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_ATTENUATED_ISM (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.dust_attenuation.PowerLaw'> |
| Apply to: young_attenuated_nebular |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - tau_v: 1.0 |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_ATTENUATED_NEBULAR (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.dust_attenuation.PowerLaw'> |
| Apply to: young_reprocessed |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - tau_v: 0.7 |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_REPROCESSED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_nebular, young_transmitted |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_NEBULAR (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_nebular_line, young_nebular_continuum |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_NEBULAR_LINE (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.escape_fraction.EscapedFraction'> |
| Apply to: _young_nebular_line_no_fesc |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - fesc_ly_alpha: 0.9 |
| Masks: |
| - log10ages < 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| _YOUNG_NEBULAR_LINE_NO_FESC (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: linecont |
| Use velocity shift: False |
| Save emission: False |
| Per particle emission: True |
| Masks: |
| - log10ages < 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_NEBULAR_CONTINUUM (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: nebular_continuum |
| Use velocity shift: False |
| Save emission: True |
| Per particle emission: True |
| Masks: |
| - log10ages < 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_TRANSMITTED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.escape_fraction.ProcessedFraction'> |
| Apply to: full_young_transmitted |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - fesc: 0.1 |
| Masks: |
| - log10ages < 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| FULL_YOUNG_TRANSMITTED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: transmitted |
| Use velocity shift: False |
| Save emission: True |
| Per particle emission: True |
| Masks: |
| - log10ages < 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_ESCAPED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.escape_fraction.EscapedFraction'> |
| Apply to: young_incident |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - fesc: 0.1 |
| Masks: |
| - log10ages < 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_INCIDENT (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: incident |
| Use velocity shift: False |
| Save emission: True |
| Per particle emission: True |
| Masks: |
| - log10ages < 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| OLD_EMERGENT (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: old_attenuated, old_escaped |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| OLD_ESCAPED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.escape_fraction.EscapedFraction'> |
| Apply to: old_incident |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - fesc: 0.1 |
| Masks: |
| - log10ages >= 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| OLD_INCIDENT (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: incident |
| Use velocity shift: False |
| Save emission: True |
| Per particle emission: True |
| Masks: |
| - log10ages >= 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| OLD_ATTENUATED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.dust_attenuation.PowerLaw'> |
| Apply to: old_reprocessed |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - tau_v: 1.0 |
|-------------------------------------------------------------------------------------------------------------|
| OLD_REPROCESSED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: old_nebular, old_transmitted |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| OLD_NEBULAR (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: old_nebular_line, old_nebular_continuum |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| OLD_NEBULAR_CONTINUUM (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: nebular_continuum |
| Use velocity shift: False |
| Save emission: True |
| Per particle emission: True |
| Masks: |
| - log10ages >= 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| OLD_NEBULAR_LINE (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.escape_fraction.EscapedFraction'> |
| Apply to: _old_nebular_line_no_fesc |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - fesc_ly_alpha: 0.9 |
| Masks: |
| - log10ages >= 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| _OLD_NEBULAR_LINE_NO_FESC (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: linecont |
| Use velocity shift: False |
| Save emission: False |
| Per particle emission: True |
| Masks: |
| - log10ages >= 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| OLD_TRANSMITTED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Transformer model: |
| Transformer: <class 'synthesizer.emission_models.transformers.escape_fraction.ProcessedFraction'> |
| Apply to: full_old_transmitted |
| Save emission: True |
| Per particle emission: True |
| Fixed parameters: |
| - fesc: 0.1 |
| Masks: |
| - log10ages >= 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| FULL_OLD_TRANSMITTED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Extraction model: |
| Grid: test_grid |
| Extract key: transmitted |
| Use velocity shift: False |
| Save emission: True |
| Per particle emission: True |
| Masks: |
| - log10ages >= 7 dimensionless |
|-------------------------------------------------------------------------------------------------------------|
| INTRINSIC (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_intrinsic, old_intrinsic |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| YOUNG_INTRINSIC (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_reprocessed, young_escaped |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| OLD_INTRINSIC (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: old_reprocessed, old_escaped |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| ATTENUATED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_attenuated_ism, old_attenuated |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| INCIDENT (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_incident, old_incident |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| TRANSMITTED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_transmitted, old_transmitted |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| NEBULAR (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_nebular, old_nebular |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| REPROCESSED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_reprocessed, old_reprocessed |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| EMERGENT (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_emergent, old_emergent |
| Save emission: True |
| Per particle emission: True |
|-------------------------------------------------------------------------------------------------------------|
| ESCAPED (stellar) |
|-------------------------------------------------------------------------------------------------------------|
|Combination model: |
| Combine models: young_escaped, old_escaped |
| Save emission: True |
| Per particle emission: True |
|=============================================================================================================|
Imaging¶
The last step before we can make any images is to define a PhotometricImager with the resolution and filters attached. We will also need the FOV, or width, of the images. Like many other inputs to Synthesizer, both the resolution and FOV must have units associated to them so the code can internally transform all quantities to a consistent unit system.
[3]:
from synthesizer.instruments import PhotometricImager
# Define the width of the image
width = 30 * kpc
# Define image resolution (here we arbitrarily set it to 200
# pixels along an axis)
resolution = width / 200
# Create an imaging instrument
instrument = PhotometricImager(
label="DemoInstrument", resolution=resolution, filters=filters
)
print(f"Image width is {width:.2f} with {resolution:.2f} resolution")
Image width is 30.00 kpc with 0.15 kpc resolution
Now we have everything we need to make images. Although it is possible to work with the low-level ImageCollection and Image methods, the main public interface is the high-level galaxy API. There are two helper methods, get_images_luminosity for luminosity images and get_images_flux for flux images. We will focus on the former here.
The image helper methods both take the image properties we previously defined, and a spectra type. These types can be any spectra for which you have calculated photometry, e.g. "incident", "intrinsic", or "attenuated". You can pass any number of these labels to get multiple images at once.
Images can either be simple 2D histograms, or the particles can be smoothed over their kernels. What type of image is made is controlled by the img_type argument. Below we demonstrate both approaches. However, for the latter we also need to define a kernel, which we have already imported from the kernel_functions module.
[4]:
# Get the SPH kernel
sph_kernel = Kernel()
kernel_data = sph_kernel.get_kernel()
# Get the image
hist_imgs = gal.get_images_luminosity(
"attenuated",
instrument=instrument,
fov=width,
img_type="hist",
kernel=kernel_data,
kernel_threshold=1,
cosmo=cosmo,
)
# Get the image
smooth_imgs = gal.get_images_luminosity(
"attenuated",
instrument=instrument,
fov=width,
img_type="smoothed",
kernel=kernel_data,
kernel_threshold=1,
cosmo=cosmo,
)
Images generated using these getter methods are both returned and attached to the Galaxy and components depending on where the photometry was calculated. In the example above, the “attenuated” spectra were calculated on the Stars component, so the images are attached to gal.stars.
Additionally, the getter methods also return either a single ImageCollection or a dictionary of ImageCollections depending on how many emission types were requested.
When making images in multiple bands, the image arrays themselves are stored on the returned ImageCollection in a dictionary called imgs, of the form {filter_code: Image}. Here, an Image is a container including the image array itself (arr), unit information (units) and the resolution and fov of the image. An Image object also includes a number of different methods for manipulating and visualising individual images.
Below we will extract this dictionary and plot each of the images we have made.
[5]:
# Lets set up a simple normalisation across all images
vmax = 0
for img in hist_imgs.values():
up = np.percentile(img.arr, 99.9)
if up > vmax:
vmax = up
hist_norm = cm.Normalize(vmin=0, vmax=vmax)
vmax = 0
for img in smooth_imgs.values():
up = np.percentile(img.arr, 99.9)
if up > vmax:
vmax = up
smooth_norm = cm.Normalize(vmin=0, vmax=vmax)
# Set up plot
fig = plt.figure(figsize=(4 * len(filters), 4 * 2))
gs = gridspec.GridSpec(2, len(filters), hspace=0.0, wspace=0.0)
# Create top row
axes = []
for i in range(len(filters)):
axes.append(fig.add_subplot(gs[0, i]))
# Loop over images plotting them
for ax, fcode in zip(axes, filter_codes):
ax.imshow(hist_imgs[fcode].arr, norm=hist_norm, cmap="Greys_r")
ax.set_title(fcode)
ax.tick_params(
axis="both",
which="both",
left=False,
right=False,
labelleft=False,
labelright=False,
bottom=False,
top=False,
labelbottom=False,
labeltop=False,
)
# Set y axis label on left most plot
axes[0].set_ylabel("Histogram")
# Create bottom row
axes = []
for i in range(len(filters)):
axes.append(fig.add_subplot(gs[1, i]))
# Loop over images plotting them
for ax, fcode in zip(axes, filter_codes):
ax.imshow(smooth_imgs[fcode].arr, norm=smooth_norm, cmap="Greys_r")
ax.tick_params(
axis="both",
which="both",
left=False,
right=False,
labelleft=False,
labelright=False,
bottom=False,
top=False,
labelbottom=False,
labeltop=False,
)
# Set y axis label on left most plot
axes[0].set_ylabel("Smoothed")
# Plot the image
plt.show()
plt.close(fig)
Applying a Point Spread Function (PSF)¶
To properly model observations from a particular instrument we must take into account the point spread function (PSF).
To apply a PSF we configure the PSFs on the PhotometricImager, then call the instrument-owned apply_psfs method on the returned ImageCollection. Here we will just create a fake gaussian PSF for all filters, but PSFs can be sourced however the user wishes, as long as a numpy array is passed for each filter when constructing the instrument.
To get the most accurate result from the PSF convolution it is recommended to do the convolution on a super-sampled image (i.e. much higher resolution than the PSF). Although-here we just supersample the images we have already made, it is recommended to first make the images at the super-sampled resolution, and then subsequently downsample after the fact.
[6]:
# Create a fake PSF for each filter
psf = np.outer(
signal.windows.gaussian(100, 3), signal.windows.gaussian(100, 3)
)
psfs = {f: psf for f in filters.filter_codes}
# Supersample the image
smooth_imgs.supersample(2)
# Attach the PSFs to an imaging instrument
psf_instrument = PhotometricImager(
label="DemoInstrumentPSF",
resolution=resolution,
filters=filters,
psfs=psfs,
)
# Apply the PSFs
psf_imgs = psf_instrument.apply_psfs(smooth_imgs)
# And downsample back to the native resolution
smooth_imgs.downsample(0.5)
psf_imgs.downsample(0.5)
apply_psfs returns a new ImageCollection containing the newly convolved Image objects. We can now use the plotting helper function to plot these images, with some normalisation.
[7]:
# Lets set up a simple normalisation across all images
vmax = 0
for img in psf_imgs.values():
up = np.percentile(img.arr, 99.9)
if up > vmax:
vmax = up
# Get the plot
fig, ax = psf_imgs.plot_images(show=True, vmin=0, vmax=vmax)
plt.close(fig)
Applying noise¶
The final ingredient for a fully forward modelled synthetic image is a noise field. We include 5 different approaches to include noise:
apply_noise_array: Add an existing noise field / array.apply_noise_from_std: Derive a noise distribution, centered on 0, given a user specified standard deviation, and then generate and add a noise array.apply_noise_from_snr(aperture): Derive a noise distribution from a Signal-to-Noise Ratio (SNR), defined in an aperture with sizeaperture_radiusand a specifieddepth. This will derive the standard deviation of the noise distribution assuming \(SNR= S / \sigma\) for an aperture, before deriving the per pixel noise, computing the noise array and adding it.apply_noise_from_snr(point source): Derive a noise distribution from a SNR and depth. This will derive the standard deviation of the noise distribution assuming \(SNR= S / \sigma\) for a pixel before computing the noise array and adding it. This behaviour can be achieved by omittingaperture_radiusin the call toapply_noise_from_snrapply_correlated_noise: Model the spatial correlation structure of an observed noise field (stored on aPhotometricImagerundernoise_source_maps) and generate a zero-mean noise realisation with the same power spectrum, which is then added to the image. The correlation function is cached on the instrument so the expensive FFT step runs only once per filter.
As with applying a PSF, these methods have singular versions (as listed above) which can be used on an individual Image, and pluralised versions which can be used on an ImageCollection. For most noise methods the pluralised versions take dictionaries keyed by filter code; apply_correlated_noise instead takes the PhotometricImager directly since the correlated-noise source maps are already stored there.
If an image has units then the passed noise_arr or noise_std must also have units.
Applying noise with any of the methods described above will return a new ImageCollection / Image containing the noisy image. In addition to the noisy image (stored under Image.arr) the new Image (or the new Image objects within an ImageCollection) will contain the noise array stored in the noise_arr attribute, and the weight map stored in the weight_map attribute on the Image.
Below we demonstrate each method using the ImageCollection interface. The noise and weight maps are stored in Image.noise_arr and Image.weight_map.
Noise arrays¶
[8]:
img_start = time.time()
# Get a noise array for each filter
noises = {
f: np.random.normal(loc=0, scale=10**26.0, size=tuple(psf_imgs.npix))
* erg
/ s
/ Hz
for f in psf_imgs.keys()
}
# Apply the noise array
noise_array_imgs = psf_imgs.apply_noise_arrays(noises)
print("Noisy images made, took:", time.time() - img_start)
# Get the plot
fig, ax = noise_array_imgs.plot_images(show=True, vmin=0, vmax=vmax)
plt.close(fig)
Noisy images made, took: 0.011090517044067383
Noise from standard deviations¶
[9]:
img_start = time.time()
# Get a noise standard deviation for each filter
noise_stds = {f: 10**26.3 * erg / s / Hz for f in psf_imgs.keys()}
# Apply the noise array
noise_std_imgs = psf_imgs.apply_noise_from_stds(noise_stds)
print("Noisy images made, took:", time.time() - img_start)
# Get the plot
fig, ax = noise_std_imgs.plot_images(show=True, vmin=0, vmax=vmax)
plt.close(fig)
Noisy images made, took: 0.009202241897583008
Noise from an aperture depth¶
[10]:
# Get dictionaries with the noise properties for each filter
snrs = {f: 5 for f in psf_imgs.keys()}
depths = {f: 10**28.0 * erg / s / Hz for f in psf_imgs.keys()}
# Apply the noise array
noise_app_imgs = psf_imgs.apply_noise_from_snrs(
snrs=snrs,
depths=depths,
aperture_radius=0.5 * kpc,
)
print("Noisy images made, took:", time.time() - img_start)
# Get the plot
fig, ax = noise_app_imgs.plot_images(show=True, vmin=0, vmax=vmax)
plt.close(fig)
Noisy images made, took: 0.13922572135925293
Noise from point source depth¶
[11]:
# Get dictionaries with the noise properties for each filter
snrs = {f: 5 for f in psf_imgs.keys()}
depths = {f: 10**27.0 * erg / s / Hz for f in psf_imgs.keys()}
# Apply the noise array
noise_ps_imgs = psf_imgs.apply_noise_from_snrs(snrs=snrs, depths=depths)
print("Noisy images made, took:", time.time() - img_start)
# Get the plot
fig, ax = noise_ps_imgs.plot_images(show=True, vmin=0, vmax=vmax)
plt.close(fig)
Noisy images made, took: 0.27416110038757324
Correlated Noise¶
Real noise fields often have correlations between pixels caused by drizzling, PSF convolution, or detector read patterns. The apply_correlated_noise method models the spatial correlation structure from an observed noise template by estimating its power spectrum, generates a new zero-mean noise realisation with the same statistical properties, and adds it to each image.
The noise template for each filter is stored on the PhotometricImager under noise_source_maps. This has two benefits:
Instrument coherence — noise characteristics belong to the imaging-instrument definition alongside PSFs, filters, and resolution.
Caching — the correlation function (CF) is computed once per
(filter_code, options)combination and cached on the instrument. When the same instrument is reused across many images (e.g. in anImageCollection), the expensive FFT step runs only once per filter regardless of how many images are processed.
This is distinct from noise_maps, which are fixed noise arrays applied directly. noise_source_maps are instead used to measure the correlated noise statistics and generate a fresh realisation when invoked.
Below we simulate a correlated noise template by smoothing Gaussian white noise in Fourier space, construct a PhotometricImager with those source maps attached, and then apply it to each filter.
[12]:
npix = psf_imgs.npix[0]
rng = np.random.default_rng(0)
white_noise = rng.normal(size=(npix, npix))
# Gaussian smoothing kernel in Fourier space (correlation length ~ 3 px)
corr_length_pix = 3.0
freq = np.fft.fftfreq(npix)
fx, fy = np.meshgrid(freq, freq)
kernel_fft = np.exp(-2 * np.pi**2 * corr_length_pix**2 * (fx**2 + fy**2))
noise_source = np.real(np.fft.ifft2(np.fft.fft2(white_noise) * kernel_fft))
# Scale to ~20 % of the image peak and attach physical units
noise_source *= 0.2 * vmax / np.max(noise_source)
# Store one correlated-noise source map per filter. In a
# real workflow these would be blank-sky cutouts from actual observations.
# Here we reuse the same template for every filter, but each could be
# different.
noise_source_maps = {f: noise_source * erg / (Hz * s) for f in psf_imgs.keys()}
correlated_noise_instrument = PhotometricImager(
label="DemoInstrumentCorrelatedNoise",
resolution=resolution,
filters=filters,
noise_source_maps=noise_source_maps,
)
# Apply correlated noise. The correlation function is computed from
# correlated_noise_instrument.noise_source_maps[filter_code] on the first
# call for each
# filter and cached on the instrument — subsequent calls (e.g. for
# additional images) skip the FFT step entirely.
correlated_noise_imgs = psf_imgs.apply_correlated_noise(
correlated_noise_instrument
)
# Get the plot
fig, ax = correlated_noise_imgs.plot_images(show=True, vmin=0, vmax=vmax)
plt.close(fig)
Making an RGB image¶
Finally we can use the RGB image method on the ImageCollection to make quick RGB images by simply providing a dictionary detailing which filters we want in which bands, and optional weights for each filter (which we will ignore here).
[13]:
# Create the rgb image from the psf image above example using our 3 filters
rgb_img = noise_ps_imgs.make_rgb_image(
rgb_filters={
"R": [
"JWST/NIRCam.F200W",
],
"G": [
"JWST/NIRCam.F150W",
],
"B": [
"JWST/NIRCam.F090W",
],
},
)
# Set up minima and maxima
vmin = -np.percentile(rgb_img, 32)
vmax = np.percentile(rgb_img, 99.9)
norm = cm.Normalize(vmin=vmin, vmax=vmax, clip=True)
print("Scaling to:", vmin, "->", vmax)
# Normalise the image.
rgb_img = norm(rgb_img)
print(rgb_img.shape)
# Plot the image
fig = plt.figure()
ax = fig.add_subplot(111)
ax.imshow(rgb_img, origin="lower", interpolation="nearest")
ax.axis("off")
plt.show()
plt.close(fig)
Scaling to: 7.092094494583863e+25 -> 3.618121590664665e+27
(200, 200, 3)
Cartesian vs Angular imaging¶
Synthesizer enables imaging in both Cartesian and angular coordinates. Which of these is used is controlled by the units of the fov and resolution parameters. If these are in angular units, then the images will be made in angular coordinates and the component coordinates and smoothing lengths (if applicable) will be converted to angular units.
Below we define the same image properties as above but this time in angular units.
[14]:
# Define the width of the image
width = 30 * cosmo.arcsec_per_kpc_proper(gal.redshift).value * arcsecond
# Define image resolution (here we arbitrarily set it to 200
# pixels along an axis)
resolution = width / 500
print(f"Image width is {width:.2f} with {resolution:.2f} resolution")
Image width is 15.73 arcsec with 0.03 arcsec resolution
The only extra thing we need to do angular imaging is pass an astropy Cosmology object to the image generator method (below get_images_luminosity). This is used to convert the coordinates and smoothing lengths of the particles to angular units. If you forget to do this, Synthesizer will raise an error.
[15]:
# Get the image
smooth_imgs = gal.get_images_luminosity(
"attenuated",
instrument=instrument,
fov=width,
img_type="smoothed",
kernel=kernel_data,
kernel_threshold=1,
cosmo=cosmo,
)
smooth_imgs.plot_images(show=True)
[15]:
(<Figure size 1400x350 with 3 Axes>, <Axes: >)
Unknown Pleasures¶
This functionality needs very little explanation. If you like Joy Division you already know.
[16]:
psf_imgs["JWST/NIRCam.F200W"].plot_unknown_pleasures()
[16]:
(<Figure size 800x800 with 1 Axes>, <Axes: >)
