galfits.images

Image Processing and Photometry Module

This module provides tools for astronomical image processing, source detection, background estimation, PSF extraction, and photometric measurements. It integrates with photutils and astropy for standard astronomical data reduction tasks.

Key Features

  • Source detection: Detect and segment sources in images

  • Background estimation: 2D background modeling and subtraction

  • PSF extraction: Build empirical PSFs from star stacks

  • Aperture photometry: Measure fluxes within various apertures

  • Image manipulation: Resampling, matching, convolution

  • Cutout extraction: Create image cutouts around targets

  • Visualization: Image display with astronomical conventions

Functions

  • normimg : Normalize image for display

  • grow_segmentation : Expand segmentation maps

  • match_psf : Match PSF between images

  • extract_psf : Build PSF from stellar sources

  • measure_background : Estimate and subtract background

  • detect_sources : Source detection and cataloging

The module is designed to work seamlessly with the main fitting pipeline, preparing data for galaxy modeling and SED fitting.

class galfits.images.image(filename=None, IFS=False, hdu=0, wcs=None, unit=None, zero_point=None, pixel_scales=None, wcs_rotation=None, mask=None, psfFWHM_map=None, fitSED=True)[source]

Bases: object

A single image object.

Functions

  • Read from fits file use CCDData.

  • get_size : Get the image size.

  • plot : Plot the image.

  • sigma_clipped_stats : Calculate the basic statistics of the image.

  • set_data : Load from numpy array.

  • set_mask : Set image mask.

  • set_pixel_scales : Set the pixel scales along two axes.

  • set_zero_point : Set magnitude zero point.

cal_sigma_image(expt, data=None, gain=1.0, fac1=1.0, fac2=1.0, filepath=None)[source]

Construct uncertainty map from photon statistics and sky background.

Calculates per-pixel uncertainties following GALFIT’s error model: quadrature sum of Poisson noise from source photons and sky background variance. Assumes input data is in count rate (e-/s) units.

Parameters:

  • expt (float) – Exposure time (seconds). Converts count rate to total counts.

  • data (ndarray, optional) – Image data in count rate units (e-/s). If None, uses self.data.data.

  • gain (float, default: 1.0) – Detector gain (e-/ADU). For pre-calibrated data in electrons, set gain=1.

  • fac1 (float, default: 1.0) – Multiplicative factor for Poisson noise term. Can account for sub-pixel resampling effects (set > 1) or correlation (set < 1).

  • fac2 (float, default: 1.0) – Multiplicative factor for sky background variance. Adjust for non-ideal background estimation.

  • filepath (str, optional) – If provided, save sigma map as FITS file to this path.

Returns:

sigmap – Uncertainty map in count rate units (e-/s). Same shape as input data.

Return type:

ndarray

Notes

Uncertainty Model:

\[\sigma_{{\rm total}} = \sqrt{{f_1 N_{{\rm source}} + f_2 \sigma_{{\rm sky}}^2}}\]

where:

  • \(N_{{\rm source}} = |I| \times G \times t\) (total electrons from source)

  • \(\sigma_{{\rm sky}}\) = sky RMS in electrons (from sigma-clipped stats)

  • \(f_1, f_2\) = empirical correction factors (fac1, fac2)

  • \(G\) = gain, \(t\) = exposure time

Final uncertainty is converted back to count rate by dividing by \(G \times t\).

The absolute value of intensity handles both positive flux and negative residuals during iterative fitting.

Typical Values:

  • HST/ACS: gain ~1-2 e-/ADU, already in e-/s → set gain=1

  • Ground-based CCD: gain ~2-5 e-/ADU, units in ADU/s → provide gain

  • Drizzled images: set fac1 ~1.5 to account for pixel correlation

This method sets self.sigma_image for later use by fitting routines.

Examples

Calculate uncertainties for HST/WFC3 (gain=1, units already e-/s):

>>> img = image('hst_wfc3.fits')
>>> sigma = img.cal_sigma_image(expt=1000.0, gain=1.0)
>>> print(f"Median S/N: {np.median(img.data.data/sigma):.1f}")

Ground-based image with gain=2.5 e-/ADU:

>>> sigma = img.cal_sigma_image(expt=600.0, gain=2.5, fac1=1.2)

Save uncertainty map for inspection:

>>> sigma = img.cal_sigma_image(expt=1000.0, filepath='sigma.fits')

See also

cal_sigma_image_of_cutout

Calculate sigma for cutout region only

read_sigmap

Load externally computed sigma map

cal_sigma_image_of_cutout(expt, gain=1.0, fac1=1.0, fac2=1.0, filepath=None)[source]

Construct sigma map following the same procedure as Galfit (quadruture sum of sigma at each pixel from source and sky background). .. note:

'GAIN' keyword must be available in the image header and ADU x GAIN = electron
assume a counts rate unit of image
Parameters:

filepath – Whether and where to save sigma map

gaussian_psf(size=15)[source]

Generate analytic Gaussian PSF kernel.

Creates a 2D Gaussian PSF normalized to unit sum, using the image’s FWHM value (self.psfFWHM_arcsec). Useful as fallback when empirical PSF cannot be built from stars.

Parameters:

size (int, default: 15) – PSF kernel size (pixels). Should be odd for symmetric centering. Typically 5-7 times FWHM in pixels is sufficient for accurate convolution.

Returns:

psf – 2D Gaussian kernel with shape (size, size), normalized to sum=1.

Return type:

ndarray

Notes

The Gaussian width is set by:

\[\sigma = {{\rm FWHM}} / (2\sqrt{{2\ln 2}}) \approx {{\rm FWHM}} / 2.355\]

where FWHM is converted from arcsec to pixels using self.pixel_scales[0].

The kernel is generated via astropy’s Gaussian2DKernel with default centering and normalization.

This method sets self.PSF for use by convolution routines in galaxy modeling.

Examples

Create Gaussian PSF for HST/WFC3 (FWHM ~ 0.15 arcsec):

>>> img.psfFWHM_arcsec = 0.15
>>> psf = img.gaussian_psf(size=21)
>>> print(psf.shape, psf.sum())
(21, 21) 1.0

Use larger kernel for ground-based seeing (FWHM ~ 0.8 arcsec):

>>> img.psfFWHM_arcsec = 0.8
>>> psf = img.gaussian_psf(size=51)  # ~6 * FWHM/pixscale

See also

generate_PSFs

Build empirical PSF from stellar sources

read_PSF

Load pre-computed PSF from file

set_psf_fwhm

Set FWHM value for PSF generation

generate_PSFs(psfFWHM, data=None, mheader=None, sigmaimg=None, nsigma=2.0, npixels=50, a_TV=1000000.0, a_sky=1.0, learning_rate=0.0001, fluxthreshould=2, isolate=0.1, qt=0.8, medianfwhm=None, sigfwhm=2, starcord_list=None, central_tol=0.0, figurepath=None, size=31, oversampling=1, plot=True, crossmatch_tolerance=None, filepaths=None, droppsf=[], use_sep=False, notfit=False, crossmatch=False, star_savepath=None, makestackpsf=False, photutils_method=False, num_steps=14000, verbose=False, referimg=None, referthre=10.0)[source]

Build empirical point spread function from stellar sources in image.

Detects point sources, extracts cutouts, and constructs high-fidelity PSF model using either photutils EPSFBuilder or custom JAX-based optimizer with TV regularization. The PSF can be oversampled for sub-pixel accuracy.

Parameters:

  • psfFWHM (float) – Expected PSF FWHM (arcsec). Used for source detection kernel and star selection criteria.

  • data (ndarray, optional) – Image data for source detection. If None, uses self.data.data.

  • sigmaimg (ndarray, optional) – Uncertainty map. If provided, used for weighted PSF fitting.

  • nsigma (float, default: 2.0) – Detection threshold (sigma above background) for candidate stars.

  • npixels (int, default: 50) – Minimum connected pixels above threshold to constitute a source.

  • a_TV (float, default: 1e6) – Total variation regularization strength (JAX method only). Higher values enforce smoother PSF. Typical range: 1e4 - 1e7.

  • a_sky (float, default: 1.0) – Sky background penalty weight (JAX method only). Prevents PSF wings from becoming negative.

  • learning_rate (float, default: 1e-4) – Adam optimizer learning rate (JAX method only).

  • fluxthreshould (float, default: 2) – Flux selection: keep stars with flux > fluxthreshould * mean(all_fluxes).

  • isolate (float, default: 0.1) – Isolation criterion: reject stars if > isolate fraction of cutout is masked by neighbors or boundaries.

  • qt (float, default: 0.8) – Axis ratio threshold: keep stars with q = b/a > qt (rejects elongated sources likely to be galaxies or blends).

  • medianfwhm (float, optional) – Override median FWHM for star selection. If None, computed from detections.

  • sigfwhm (float, default: 2) – FWHM selection: keep stars within sigfwhm * std of median FWHM, or set None to disable.

  • starcord_list (list of [RA, Dec], optional) – Manually specified star coordinates (degrees). If provided, uses these instead of automatic detection.

  • central_tol (float, default: 0.0) – Fraction of central region (0.25*size radius) allowed to be masked. Rejects stars with bad pixels near peak.

  • figurepath (str, optional) – Directory path for saving diagnostic plots. If None and plot=True, shows interactively.

  • size (int, default: 31) – Cutout size (pixels) for star extraction. Should encompass PSF wings.

  • oversampling (int, default: 1) – PSF oversampling factor. E.g., oversampling=2 creates PSF with 2x finer pixel grid. Useful for dithered observations.

  • plot (bool, default: True) – If True, display cutouts of all stars and final PSF.

  • filepaths (str, optional) – FITS file path to save final PSF.

  • droppsf (list of int, optional) – Indices of stars to exclude from PSF building (after initial filtering).

  • notfit (bool, default: False) – If True, return None after showing star cutouts (diagnostic mode).

  • crossmatch (bool, default: False) – If True, cross-match detections with Gaia DR3 to confirm stars. Requires astroquery.

  • star_savepath (str, optional) – Directory path to save individual star cutouts as FITS files (for debugging).

  • makestackpsf (bool, default: False) – If True, return simple median-stacked PSF instead of fitting.

  • photutils_method (bool, default: False) – If True, use photutils EPSFBuilder. If False, use JAX-based fitter.

  • num_steps (int, default: 14000) – Number of optimization iterations (JAX method only).

  • verbose (bool, default: False) – If True, print loss every 100 iterations (JAX method).

  • referimg (ndarray, optional) – Reference PSF for quality check. Stars with chi-square > referthre when fit with referimg are rejected.

  • referthre (float, default: 10.0) – Reduced chi-square threshold for reference PSF filtering.

Returns:

epsf – Empirical PSF with shape (size*oversampling, size*oversampling). Normalized to unit sum.

Return type:

ndarray

Notes

Star Selection Pipeline:

  1. Detection: Use photutils detect_sources with Gaussian kernel

  2. Gaia cross-match (optional): Confirm point sources

  3. Flux cut: Keep bright sources (> fluxthreshould * mean_flux)

  4. Morphology cuts: - Axis ratio q > qt (reject elongated sources) - FWHM within sigfwhm sigma of median (reject outliers)

  5. Isolation check: Fraction masked < isolate

  6. Central quality: Central region masking < central_tol

  7. Reference fit (optional): Chi-square with referimg < referthre

PSF Building Methods:

  • Photutils (photutils_method=True): Anderson & King (2000) ePSF algorithm with 150 iterations

  • JAX (photutils_method=False): Simultaneous fit of PSF + positions + normalizations using Adam optimizer with total variation regularization:

    Loss = 0.5 * sum((data - model)^2 / sigma^2)
       + a_TV * TV(PSF)
       + a_sky * sum(PSF[outside_circle]^2)

The JAX method typically produces sharper PSFs by avoiding smoothing inherent in interpolation, but requires careful tuning of a_TV and a_sky.

This method sets self.PSF for use in galaxy model convolution.

Examples

Build PSF from automatic star detection:

>>> img = image('science.fits')
>>> img.cal_sigma_image(expt=1000.0, gain=1.0)
>>> psf = img.generate_PSFs(psfFWHM=0.15, sigmaimg=img.sigma_image,
...                         size=31, oversampling=2)
>>> print(psf.shape)
(62, 62)

Use manual star list (e.g., from catalog):

>>> star_coords = [[150.123, 2.456], [150.145, 2.478], ...]
>>> psf = img.generate_PSFs(psfFWHM=0.15, starcord_list=star_coords,
...                         crossmatch=False)

Build PSF with stricter selection and save diagnostics:

>>> psf = img.generate_PSFs(psfFWHM=0.12, fluxthreshould=3, qt=0.9,
...                         sigfwhm=1.5, filepaths='psf.fits',
...                         figurepath='./psf_diagnostics/')

Quick stacked PSF (no fitting):

>>> psf = img.generate_PSFs(psfFWHM=0.15, makestackpsf=True)

See also

gaussian_psf

Create analytic Gaussian PSF

fittingPSF

Low-level JAX PSF fitter

photutils.psf.EPSFBuilder

Photutils ePSF construction

generate_cutmask(psfFWHM, data=None, sdmask=None, nsigma=3.0, magnification=3.0, npixels=5, nlevels=32, contrast=0.001, deblend=False, sky_level=None, source_dia=0.5, max_grow_size=100, growmethod=True, nearbygrow=0, addgrow=0, left_source=True)[source]

Generate pixel mask for contaminating sources in cutout region.

Detects and masks background/foreground sources to isolate the target object for fitting. Uses photutils segmentation to identify sources, then optionally grows segmentation regions adaptively based on local flux levels.

Parameters:

  • psfFWHM (float) – Point spread function FWHM (arcsec). Used to set detection kernel size and minimum source separation.

  • data (ndarray, optional) – Image data to analyze. If None, uses self.cut_image.

  • sdmask (ndarray of bool, optional) – Pre-existing mask (e.g., bad pixels). Pixels marked True are excluded from source detection.

  • nsigma (float, default: 3.0) – Detection threshold in units of background RMS. Sources must have peak flux > nsigma * sky_std.

  • magnification (float, default: 3.0) – Factor by which to expand detected source regions when growmethod=False. Mask semi-major axis becomes magnification * semimajor_sigma.

  • npixels (int, default: 5) – Minimum connected pixels above threshold to constitute a source.

  • nlevels (int, default: 32) – Number of multi-thresholding levels for deblending (if deblend=True).

  • contrast (float, default: 0.001) – Deblending contrast parameter. Smaller values → more aggressive separation.

  • deblend (bool, default: False) – If True, attempt to separate blended sources using multi-level segmentation.

  • sky_level (float, optional) – Flux threshold for adaptive growth termination. If None, set to 3*sky_std. Growth stops when median flux of new pixels < sky_level.

  • source_dia (float, default: 0.5) – Protected region diameter (arcsec) around target source position. Sources within this distance are treated specially (not masked if left_source=True).

  • max_grow_size (int, default: 100) – Maximum number of dilation iterations for adaptive growth.

  • growmethod (bool, default: True) – If True, use adaptive growth algorithm (grow_segmentation). If False, simply expand all sources by magnification factor using their ellipse parameters.

  • nearbygrow (int, default: 0) – Number of dilation iterations applied to sources near target (within source_dia). Set > 0 to also mask nearby companions.

  • addgrow (int, default: 0) – Additional dilation iterations applied to all contaminating sources after adaptive growth converges.

  • left_source (bool, default: True) – If True, protect the central source at self.sources_skycord_cut from being masked. If False, mask all detected sources.

Returns:

mask – Boolean mask with shape matching data. True indicates pixels to exclude from fitting (contaminating sources, bad pixels, NaN/Inf).

Return type:

ndarray of bool

Notes

Adaptive Growth Algorithm (when growmethod=True):

  1. Detect sources using photutils.detect_sources with Gaussian kernel

  2. Identify target source label at self.sources_skycord_cut

  3. For each contaminating source: - Iteratively dilate segmentation by 1 pixel - Stop when median of new pixels < sky_level OR iteration > max_grow_size - Apply additional addgrow iterations

  4. Combine all masks, excluding target source region

Elliptical Masking (when growmethod=False):

  • Measure ellipse parameters (a, b, PA) for each source

  • Mask ellipse with semi-major axis = magnification * a

  • Preserve target source by checking overlap with protected circle

The adaptive method is recommended for crowded fields as it follows source isophotes more accurately than fixed ellipses.

Examples

Basic contamination masking for isolated galaxy:

>>> img.img_cut(ra=150.123, dec=2.456, cutsize=20.0)
>>> mask = img.generate_cutmask(psfFWHM=0.15, nsigma=2.0)
>>> print(f"Masked {np.sum(mask)/mask.size*100:.1f}% of pixels")

Crowded field with aggressive masking:

>>> mask = img.generate_cutmask(psfFWHM=0.12, nsigma=2.5, magnification=5.0,
...                             max_grow_size=50, addgrow=5)

Mask all sources including target (for background estimation):

>>> mask_all = img.generate_cutmask(psfFWHM=0.15, left_source=False)

See also

grow_segmentation

Adaptive segmentation growth algorithm

sky_subtraction

Background subtraction using source masking

get_data_info()[source]

Data information to generate model image.

Returns:

d –

shape(ny, nx)

Image array shape.

pixel_scale(pixelscale_x, pixelscale_y), default units: arcsec

Pixel scales.

wcs_rotationangle, default units: radian

WCS rotation, east of north.

Return type:

dict

get_size(units='pixel')[source]

Get the size of the image.

Parameters:

units (string) – Units of the size (pixel or angular units).

Returns:

x, y – Size along X and Y axes.

Return type:

float

img_cut(ra, dec, cutsize, cutposition=None, IFS=False, usemaxpix=False, sigma_clipped=True, referpix=None, zoomsize=100, spcorrection=[0.0, 0.0])[source]

Extract image cutout around a celestial position.

This method creates a cutout from the full image centered on specified RA/Dec coordinates, handles coordinate transformation between world and pixel coordinates, and calculates transformation parameters needed for model generation. It can optionally recenter on the brightest pixel within a search region.

Parameters:

  • ra (float) – Right ascension of the source position (degrees).

  • dec (float) – Declination of the source position (degrees).

  • cutsize (float or array-like) – Cut size in arcsec. If float, symmetric box. If array [xmin, xmax, ymin, ymax], asymmetric cutout around source position.

  • cutposition (list of int, optional) – Direct pixel coordinates [minx, maxx, miny, maxy]. If provided, overrides cutsize calculation.

  • IFS (bool, default: False) – If True, treat data as IFS (Integral Field Spectrograph) cube with shape (nwave, ny, nx).

  • usemaxpix (bool, default: False) – If True, recenter cutout on brightest pixel within zoomsize of initial position. Useful when catalog position is inaccurate.

  • sigma_clipped (bool, default: True) – If True, estimate sky median using sigma-clipped statistics. Otherwise use simple median.

  • referpix (list of float, optional) – Reference pixel [x, y] for centering search when usemaxpix=True. If None, uses WCS-derived position from RA/Dec.

  • zoomsize (float, default: 100) – Search radius in arcsec for brightest pixel when usemaxpix=True.

  • spcorrection (list of float, default: [0., 0.]) – Manual position correction [dx, dy] in pixels to apply to WCS-derived position before cutout.

Returns:

  • imcut (jax.numpy.ndarray) – Cutout image with shape (ny, nx) or (nwave, ny, nx) for IFS data.

  • cutposition (list of int) – Pixel bounds [minx, maxx, miny, maxy] of the cutout in the full image.

Notes

This method sets several instance attributes used for modeling:

  • cut_image : The extracted cutout

  • sources_skycord : [x, y] position in full image pixel coordinates

  • sources_skycord_cut : [x, y] position in cutout pixel coordinates

  • coordinates_transfer_para : Dictionary with coordinate transformation info:

    {
    'x0': float,        # X position in cutout (pixels)
    'x0shift': float,   # X alignment offset (pixels)
    'y0': float,        # Y position in cutout (pixels)
    'y0shift': float,   # Y alignment offset (pixels)
    'dxra': float,      # factor of ra contribute to x axis
    'dxdec': float,     # factor of dec contribute to x axis
    'dyra': float,      # factor of ra contribute to y axis
    'dydec': float,     # factor of dec contribute to y axis
    'pixsc': float,     # Pixel scale (arcsec/pixel)
    'delta_ang': float  # Position angle east of north (degrees)
}

  • sky_median : Background level estimated from cutout

  • cut_mask_image : Mask applied to cutout (if data has mask)

  • cut_sigma_image : Uncertainty map for cutout (if sigma_image exists)

  • dcord : [dx, dy] offset of CRPIX from cutout origin

The transformation parameters are calculated from WCS information and allow conversion between sky offsets (arcsec) and pixel coordinates in arbitrary orientation, essential for modeling rotated galaxies.

Examples

Extract 10 arcsec cutout around a galaxy:

>>> img = image('data.fits')
>>> cutout, bounds = img.img_cut(ra=150.123, dec=2.456, cutsize=10.0)
>>> print(cutout.shape)
(100, 100)  # if pixel scale ~0.1 arcsec

Asymmetric cutout extending further in +x direction:

>>> cutout, bounds = img.img_cut(ra=150.123, dec=2.456,
...                              cutsize=[5, 15, 10, 10])

Recenter on brightest pixel within 2 arcsec:

>>> cutout, bounds = img.img_cut(ra=150.123, dec=2.456, cutsize=20.0,
...                              usemaxpix=True, zoomsize=2.0)

See also

write_imcut

Save cutout and ancillary data to FITS file

read_imcut

Load previously saved cutout

img_cut_align(ra, dec, cutsize, cutposition=None, IFS=False, usemaxpix=False, sigma_clipped=True, referpix=None, zoomsize=100, spcorrection=[0.0, 0.0])[source]

Extract a cutout around a celestial position, accurately resampled to a strict North-Up, East-Left standard WCS grid.

plot(stretch='asinh', units='arcsec', vmin=None, vmax=None, a=None, ax=None, plain=False, **kwargs)[source]

Plot an image.

Parameters:

  • stretch (string (default: 'asinh')) – Choice of stretch: asinh, linear, sqrt, log.

  • units (string (default: 'arcsec')) – Units of pixel scale.

  • (optional) (ax) – Minimal value of imshow.

  • (optional) – Maximal value of imshow.

  • (optional) – Scale factor of some stretch function.

  • (optional) – Axis to plot the image.

  • plain (bool (default: False)) – If False, tune the image.

  • **kwargs (Additional parameters goes into plt.imshow())

Returns:

ax – Axis to plot the image.

Return type:

matplotlib.Axis

plot_direction(ax, xy=(0, 0), len_E=None, len_N=None, color='k', fontsize=20, linewidth=2, frac_len=0.1, units='arcsec', backextend=0.05)[source]

Plot the direction arrow. Only applied to plots using WCS. :param ax: Axis to plot the direction. :type ax: Axis :param xy: Coordinate of the origin of the arrows. :type xy: (x, y) :param length: Length of the arrows, units: pixel. :type length: float :param units: Units of xy. :type units: string (default: arcsec)

read_PSF(filepath)[source]

read PSF image from “filepath”

read_imcut(read_hdu, transpar, cutimage=0, sigmaimage=None, sigma_clipped=True)[source]
read_sigmap(filepath)[source]

read sigma image from “filepath”

read_ss_image(filepath)[source]

read sky subtracted image from “filepath”

reproject(ra, dec, r, oversampling=1)[source]

Extract and resample image cutout to arbitrary pixel scale and orientation.

Uses reproject.reproject_interp to create cutout centered on specified coordinates with custom pixel scale. Useful for matching images to common grid, zooming into regions, or aligning data from different instruments.

Parameters:

  • ra (float) – Right ascension of cutout center (degrees).

  • dec (float) – Declination of cutout center (degrees).

  • r (int) – Half-width of cutout (pixels in output frame). Output shape will be (2*r+1, 2*r+1).

  • oversampling (float, default: 1) – Pixel scale reduction factor. Output pixel scale = input pixel scale / oversampling. E.g., oversampling=2 makes pixels twice as fine (0.1” → 0.05”).

Returns:

imcut – Reprojected image cutout with shape (2*r+1, 2*r+1). Pixels outside original image footprint are NaN.

Return type:

ndarray

Notes

WCS Construction:

Creates target WCS with: - CRVAL = (RA, Dec) - cutout center - CRPIX = (r, r) - places center at image center - CDELT = ±pixel_scale / oversampling (degrees/pixel) - CTYPE = “RA—TAN”, “DEC–TAN” (tangent plane projection)

Interpolation Method:

Uses reproject_interp which: - Transforms pixel coordinates between input and output WCS - Performs bi-linear interpolation of flux values - Conserves surface brightness (not flux) - Handles rotation and scaling simultaneously

This is not a simple cutout - it resamples the data, so noise properties change. For analysis requiring independent pixels, consider drizzling instead.

Currently only supports 2D images (not IFS cubes).

Examples

Extract 5 arcsec cutout at native resolution:

>>> img = image('hst_wfc3.fits')  # pixel scale 0.06 arcsec
>>> cutout = img.reproject(ra=150.123, dec=2.456, r=42, oversampling=1)
>>> print(cutout.shape)
(85, 85)  # 85 pixels * 0.06 arcsec/pix = 5.1 arcsec

Zoom in 4x to resolve fine structure:

>>> cutout_zoom = img.reproject(ra=150.123, dec=2.456, r=84,
...                             oversampling=4)
>>> print(cutout_zoom.shape)
(169, 169)  # pixel scale now 0.015 arcsec

Match JWST NIRCam (0.03”/pix) to HST/ACS (0.05”/pix):

>>> jwst = image('jwst_nircam.fits')
>>> cutout = jwst.reproject(ra=150.123, dec=2.456, r=100,
...                         oversampling=0.6)  # 0.03/0.6 = 0.05

See also

img_cut

Fast native-resolution cutout extraction

reproject.reproject_interp

Underlying reprojection function

set_data(data, unit)[source]
Parameters:

  • data (2D array) – Image data.

  • unit (string) – Unit for CCDData.

set_mask(mask)[source]

Set mask for the extension.

Parameters:

mask (2D array) – The mask.

set_pixel_scales(pixel_scales)[source]
Parameters:

(optional) (pixel_scales) – Pixel scales along the first and second directions, units: arcsec.

set_psf_fwhm(fwhm=None)[source]
set_zero_point(zp)[source]

Set magnitude zero point.

sigma_clipped_stats(**kwargs)[source]

Run astropy.stats.sigma_clipped_stats to get the basic statistics of the image.

Parameters:

astropy.stats.sigma_clipped_stats(). (All of the parameters go to)

Returns:

mean, median, stddev – The mean, median, and standard deviation of the sigma-clipped data.

Return type:

float

sky_subtraction(sdmask=None, polyfit=True, order=3, filepath=None, psfFWHM=1.5, nsigma=3.0, npixels=5, max_grow_size=100, addgrow=0, usemedian=False)[source]

Estimate and subtract sky background from image.

Masks sources using segmentation, then fits 2D polynomial or computes median/2D background model from unmasked pixels. Useful for removing spatially varying sky before galaxy modeling.

Parameters:

  • sdmask (ndarray of bool, optional) – Pre-existing mask (e.g., bad pixels). Combined with detected source mask.

  • polyfit (bool, default: True) – If True, fit 2D polynomial to unmasked sky pixels. If False, use photutils Background2D with median filter boxes (unless usemedian=True).

  • order (int, default: 3) – Polynomial order for polyfit=True. Order 1 = plane, 2 = paraboloid, etc. Higher orders can capture complex gradients but may overfit.

  • filepath (str, optional) – If provided, save sky-subtracted image to this FITS file.

  • psfFWHM (float, default: 1.5) – PSF FWHM (arcsec) for source detection kernel.

  • nsigma (float, default: 3.0) – Detection threshold (sigma) for masking sources.

  • npixels (int, default: 5) – Minimum connected pixels to constitute a source.

  • max_grow_size (int, default: 100) – Maximum dilation iterations for growing source masks.

  • addgrow (int, default: 0) – Additional dilation iterations after adaptive growth.

  • usemedian (bool, default: False) – If True, subtract simple median of unmasked pixels (fastest, least accurate).

Returns:

maskplus – Combined mask of sources, bad pixels, and boundary artifacts.

Return type:

ndarray of bool

Notes

Background Subtraction Methods:

  1. Median (usemedian=True): - Background = median(unmasked pixels) - Fastest, suitable for flat fields

  2. 2D Polynomial (polyfit=True): - Fits polynomial surface to unmasked pixels via least squares - Good for gradients from scattered light, moon, etc. - Equation for order=2: \(B(x,y) = c_0 + c_1 x + c_2 y + c_3 x^2 + c_4 xy + c_5 y^2\)

  3. 2D Median Filter (polyfit=False, usemedian=False): - Uses photutils Background2D with 30x30 pixel boxes - Adapts to local sky variations - Slower but handles complex backgrounds

Masking Strategy:

  • Detects sources with detect_sources and detect_threshold

  • Grows segmentation masks to exclude faint wings and artifacts

  • Additionally masks pixels with values exactly equal to the 4 image corners (often indicates bad data or padding)

The sky-subtracted data overwrites self.data.data. Original data is lost unless saved beforehand.

Examples

Polynomial sky subtraction (most common):

>>> img = image('science.fits')
>>> mask = img.sky_subtraction(psfFWHM=0.15, polyfit=True, order=2)
>>> print(f"Sky median after subtraction: {np.median(img.data.data):.6f}")

2D median filter for variable background:

>>> mask = img.sky_subtraction(psfFWHM=0.15, polyfit=False, usemedian=False)

Simple median subtraction for quick processing:

>>> mask = img.sky_subtraction(psfFWHM=0.15, usemedian=True)

Save sky-subtracted image:

>>> mask = img.sky_subtraction(psfFWHM=0.15, order=3,
...                            filepath='science_skysub.fits')

See also

generate_cutmask

Source masking algorithm

polynomialfit

2D polynomial fitting function

photutils.background.Background2D

Photutils 2D background estimator

source_detection_individual(psfFWHM, nsigma=3.0, sc_key='')[source]
Parameters:

  • psfFWHM (float) – FWHM of the imaging point spread function

  • nsigma (float) – source detection threshold

source_detection_residual(psfFWHM, nsigma=3.0, sc_key='')[source]

Copied from source_detection_individual but detects both positive and negative sources

write_imcut(ra, dec, outputname, header=None, addpixsc=False, addmodel=False, Alamda=None)[source]

Write cutout image and ancillary data to multi-extension FITS file.

Saves the cutout data, mask, uncertainty, and optionally model to a FITS file with proper WCS keywords. Useful for archiving intermediate results or manual inspection of fitting regions.

Parameters:

  • ra (float, optional) – Right ascension of source (degrees). If None, uses stored self.skycord[0].

  • dec (float, optional) – Declination of source (degrees). If None, uses stored self.skycord[1].

  • outputname (str) – Output FITS file path.

  • header (astropy.io.fits.Header, optional) – Base header for primary HDU. If None, uses original image header.

  • addpixsc (bool, default: False) – If True, add CDELT1/CDELT2 keywords with pixel scale (degrees/pixel).

  • addmodel (bool, default: False) – If True, write 5-extension FITS with residual (data-model) as primary HDU, plus mask, sigma, model, and original data extensions.

  • Structure (FITS)

  • --------------

  • extensions) (With model (5) –

    • HDU 0: Cut image data (PRIMARY, float32)

    • HDU 1: Cut mask (IMAGE, float32, bool converted to float)

    • HDU 2: Cut sigma/uncertainty (IMAGE, float32)

  • extensions) –

    • HDU 0: Residual = data - model (PRIMARY, float32)

    • HDU 1: Cut mask (IMAGE, float32)

    • HDU 2: Cut sigma/uncertainty (IMAGE, float32)

    • HDU 3: Model image (IMAGE, float32)

    • HDU 4: Original cut image (IMAGE, float32)

Notes

The primary HDU header is updated with correct WCS reference pixel (CRPIX1/2) to reflect the cutout’s position within the original image. This allows proper coordinate conversion even for cutouts.

Examples

Save cutout with default 3-extension format:

>>> img.img_cut(ra=150.123, dec=2.456, cutsize=10.0)
>>> img.write_imcut(ra=150.123, dec=2.456, outputname='cutout.fits')

Save cutout with model and updated pixel scale:

>>> img.model_image = create_model(...)  # after fitting
>>> img.write_imcut(ra=150.123, dec=2.456, outputname='cutout_model.fits',
...                 addpixsc=True, addmodel=True)

See also

img_cut

Extract cutout from full image

read_imcut

Load previously saved cutout

write_imcut_align(ra=None, dec=None, outputname='cutout.fits', header=None, addpixsc=False, addmodel=False, Alamda=None)[source]

Write cutout image and ancillary data to multi-extension FITS file. (Updated for strictly aligned standard WCS).

write_model(filepath)[source]

Write model image to a fits file.

class galfits.images.image_atlas(image_list=None, band_list=None, grism_image_list=None, grism_band_list=None, same_grid=True, refere_agn=0, Nlogv=2000, name='')[source]

Bases: object

Container for multi-band imaging and spectroscopy observations.

Groups multiple image objects (photometric bands) and optional grism spectra for simultaneous multi-band SED fitting. Handles cross-matching catalogs, coordinated masking, PSF building, and residual tracking across all bands.

Parameters:

  • image_list (list of image, optional) – List of image objects, one per photometric band. Should all observe the same target region (though pixel grids may differ).

  • band_list (list of str, optional) – Band identifiers (e.g., [‘F435W’, ‘F606W’, ‘F814W’]). Must match length of image_list. Check instrument_info for standard naming.

  • grism_image_list (list of image, optional) – List of image objects for grism/IFU spectroscopy. Spectral dimension handled separately from photometry.

  • grism_band_list (list of str, optional) – Identifiers for grism observations.

  • same_grid (bool, default: True) – If True, all images share the same pixel grid (pre-aligned). Enables faster processing. If False, coordinate transformations handled per-image.

  • refere_agn (int, default: 0) – Index of reference band for AGN modeling (if fitting AGN+host decomposition).

  • Nlogv (int, default: 2000) – Number of logarithmic velocity bins for emission line fitting.

  • name (str, default: '') – Dataset identifier for tracking multiple image atlases.

image_list

Photometric band images.

Type:

list of image

band_list

Band names.

Type:

list of str

common_catalog

Cross-matched source catalog from make_common_catalog(). Contains positions, shapes, and master (maximum across bands) semi-major/minor axes.

Type:

astropy.table.Table

spectra

Attached spectra from add_spectrum(). Can be fiber, longslit, or IFU.

Type:

list

specapermask

Aperture masks for spectral extraction, one per spectrum.

Type:

list of ndarray

img_residual

Chi-square residuals per band after fitting, shape (N_bands, 1).

Type:

ndarray

img_dof

Degrees of freedom per band, shape (N_bands, 1).

Type:

ndarray

spec_residual

Chi-square residuals per spectrum.

Type:

list

spec_dof

Degrees of freedom per spectrum.

Type:

list

add_spectrum(spectrum, referimag, sposition, sapersize, ...)[source]

Attach spectroscopic data with aperture definition.

source_detection(nsigma)[source]

Detect sources in all bands.

make_common_catalog(CM_separation, magnification)[source]

Cross-match multi-band catalogs.

master_mask(magnification, applylist)[source]

Apply unified masking based on common catalog.

sky_subtraction(order, filepaths)[source]

Background subtraction across all bands.

generate_PSFs(equivalent_radius, size, oversampling, ...)[source]

Build PSF from stars common to all bands.

Notes

Typical Workflow:

  1. Create image objects for each band:

    img_f435w = image('f435w.fits', zero_point=25.68, ...)
img_f606w = image('f606w.fits', zero_point=26.51, ...)
atlas = image_atlas([img_f435w, img_f606w], ['F435W', 'F606W'])

  2. Sky subtraction and source detection:

    atlas.sky_subtraction(order=3)
atlas.source_detection(nsigma=3.0)

  3. Cross-match and mask contaminating sources:

    atlas.make_common_catalog(CM_separation=2.5)
atlas.master_mask(magnification=3.0)

  4. Build PSF from common stars:

    atlas.generate_PSFs(equivalent_radius=1.5, size=31, oversampling=2)

  5. Cut target regions and attach to fitter:

    for img in atlas.image_list:
    img.img_cut(ra=150.123, dec=2.456, cutsize=10.0)

Residual Tracking:

After fitting, img_residual and spec_residual store chi-square values for each band/spectrum, enabling identification of poor fits or missing model components.

Examples

Create atlas from HST/ACS multi-band images:

>>> from galfits.images import image, image_atlas
>>> f435w = image('acs_f435w.fits', zero_point=25.68)
>>> f606w = image('acs_f606w.fits', zero_point=26.51)
>>> f814w = image('acs_f814w.fits', zero_point=25.94)
>>> atlas = image_atlas([f435w, f606w, f814w], ['F435W', 'F606W', 'F814W'])
>>> print(len(atlas))
3

Access individual bands:

>>> img_v = atlas['F606W']
>>> img_v.sigma_clipped_stats()

See also

image

Single-band image class

galfitS_data

Multi-atlas container for handling multiple targets/instruments

add_spectrum(spectrum, referimag=-1, sposition=[0, 0], sapersize=2.0, prefix='', stype='fiber', PA=0, a=1, b=1, hires=False, norm=1.0, lsf=None)[source]

Attach spectroscopic observation with aperture definition.

Links spectrum to photometric reference image and defines extraction aperture geometry. Enables joint image+spectrum SED fitting with spatially-matched constraints.

Parameters:

  • spectrum (array-like or spectrum object) – Spectroscopic data (wavelength, flux, error arrays or spectrum object).

  • referimag (int, default: -1) – Index of reference image in self.image_list for aperture definition. If -1, no aperture mask created (whole image).

  • sposition (list of float, default: [0, 0]) – Aperture center in sky offset coordinates (arcsec) from reference image cutout center. [0,0] = use galaxy center from photometry.

  • sapersize (float, default: 2.0) – Aperture size parameter (arcsec). Interpretation depends on stype: - ‘fiber’: radius - ‘longslit’: half-width perpendicular to slit

  • prefix (str, default: '') – Identifier prefix for this spectrum (e.g., MOS1_, IFU_).

  • stype ({'fiber', 'longslit'}, default: 'fiber') – Aperture geometry type: - ‘fiber’: circular aperture - ‘longslit’: rectangular slit

  • PA (float, default: 0) – Position angle of longslit (degrees, east of north). Only used if stype=’longslit’.

  • a (float, default: 1) – Longslit half-length along slit direction (arcsec). Only used if stype=’longslit’.

  • b (float, default: 1) – Longslit half-width perpendicular to slit (arcsec). Only used if stype=’longslit’.

  • hires (bool, default: False) – If True, use high-resolution spectral templates (e.g., for echelle data).

  • norm (float, default: 1.0) – Flux normalization factor applied to spectrum before fitting.

  • lsf (array-like or None, optional) – 2D array defining line spread function (LSF) for spectral [wave, R]

Notes

Aperture Masking:

Creates boolean mask with shape matching reference image cutout. Mask is flattened and stored in self.specapermask for efficient aperture integration during model generation.

For fiber aperture:

\[{{\rm mask}} = \sqrt{{(x - x_c)^2 + (y - y_c)^2}} < r_{{\rm aper}}\]

For longslit:

\[{{\rm mask}} = (|{{\bf r}} \cdot {{\bf \hat{{s}}}}| < a) \land (|{{\bf r}} \cdot {{\bf \hat{{n}}}}| < b)\]

where \({{\bf \hat{{s}}}}\) is unit vector along slit (PA) and \({{\bf \hat{{n}}}}\) is perpendicular.

Residual Tracking:

Initializes spec_residual and spec_dof entries for this spectrum, updated during fitting.

Examples

Add fiber spectrum centered on galaxy:

>>> from galfits.spectrum import spectrum
>>> spec = spectrum('muse_spectrum.fits')
>>> atlas.add_spectrum(spec, referimag=0, sposition=[0,0],
...                    sapersize=1.5, stype='fiber')

Add off-center longslit along PA=45°:

>>> atlas.add_spectrum(spec, referimag=0, sposition=[0.5, -0.3],
...                    stype='longslit', PA=45, a=3.0, b=0.5)

Add spectrum without spatial matching (no aperture):

>>> atlas.add_spectrum(spec, referimag=-1)

See also

image.img_cut

Define cutout region for aperture mapping

generate_PSFs(equivalent_radius=1.5, threshold=50.0, axrit=0.75, size=20.0, oversampling=1, plot=None, filepaths=None, photutils=False)[source]

Generate effective point spread fuctions (ePSFs) for each image :param equivalent_radius: radius criteria to indentify star: fwhm of sourve < equivalent_radius*psfFWHM :type equivalent_radius: float, :param size: use what size box to extract stars :type size: float, unit pixel :param oversampling: oversample the ePSF :type oversampling: int :param plot: list like [1,2,3] to plot rgb image :type plot: None for not plot stars & ePSF :param filepaths: :type filepaths: filepath to store the ePSFs

init_sed()[source]

Initialize the SED model

make_common_catalog(CM_separation=2.5, magnification=3.0, applylist=None)[source]

Do multi-band source detection :param CM_separation: :type CM_separation: float :param angular separation used to do sky coordinates crossmatching: :param unit in deg: :param magnification: magnification for generating mask foe each image :type magnification: float, or a array with same size as image_atlas :param applylist: :type applylist: [list of index] :param None for all images:

master_mask(magnification=3.0, applylist=None)[source]

Do multi-band source masking :param magnification: magnification for generating mask foe each image :type magnification: float, or a array with same size as image_atlas :param applylist: :type applylist: [list of index] :param None for all images:

sky_subtraction(order=3, filepaths=None)[source]

Do multi-band sky subtration :param order (optional): order of the polynomial :type order (optional): int :param filepaths: :type filepaths: filepath to store the sky subtracted images

source_detection(nsigma=3.0)[source]

Do multi-band source detection :param nsigma: source detection threshold :type nsigma: float, or a array with same size as image_atlas