BackStats

BackStats - Background Statistics Calculator

Space Telescope Science Institute

Command Line Usage

BackStats.py [-np 8] [-all] [-rm] [-sigma_clipped] field T01 T02 ...

Required Arguments:

field

Field name (e.g., ‘LMC_c42’)

T01 T02 …

One or more tile identifiers (unless -all is specified)

Optional Arguments:

-all

Process all tiles (T01 through T16)

-rm

Remove files created by BackPrep after successful completion. Only occurs at the end of the program so that if something fails, no background files will have been deleted.

-np N

Perform calculations in parallel using N threads (default: 1)

-sigma_clipped

Diagnostic mode that calculates sigma-clipped versions of mean, median, and std. Greatly increases run time but provides robust statistics.

Output

Results written to Summary/{field}_{tile}_xxx.txt containing:

  • Overlapping file pairs

  • Number of overlapping pixels

  • Mean, median, std of differences

  • Sigma-clipped statistics (if requested)

  • Skewness and kurtosis

  • Mode (half-sample algorithm)

  • Filter and exposure time metadata

Performance Notes This routine can take considerable time because:

  • Each tile involves many files (often 100-1000)

  • Must check all pairwise combinations (N × (N-1) matches)

  • For 100 files: ~10,000 comparisons

  • For 1000 files: ~1,000,000 comparisons

Optimization:

Running with maximum available threads is strongly recommended. Use -np 0 for auto-detection or explicitly set thread count based on your system.

Background Level Determination The mode (calculated using half-sample algorithm) is used for setting background levels between images. This is more robust than mean or median for skewed distributions common in astronomical images.

Notes

This routine can take considerable time because:

  • Each tile involves many files (often 100-1000)

  • Must check all pairwise combinations (N × (N-1) matches)

  • For 100 files: ~10,000 comparisons

  • For 1000 files: ~1,000,000 comparisons

Optimization:

Running with maximum available threads is strongly recommended. Use -np 0 for auto-detection or explicitly set thread count based on your system.

Background Level Determination The mode (calculated using half-sample algorithm) is used for setting background levels between images. This is more robust than mean or median for skewed distributions common in astronomical images.

Version History

230505 ksl

Coding begun and routine parallelized

230617 ksl

Adapted to new version of routines

230704 ksl

Mode and additional statistics added. Mode now used for background levels.

Author Space Telescope Science Institute

Attributes

BACKDIR

ierror

Functions

calculate_one(file, xmatch_files, indir[, ...])

Calculate background statistics for one file vs multiple matches.

do_one_tile([field, tile, nproc, calc_sigma_clipped])

Process all overlapping file pairs for one tile.

halfsamplemode(inputData[, axis])

Compute mode of an array using the half-sample algorithm.

quartile_stats(xdata)

Compute skewness and kurtosis with outlier rejection.

steer(argv)

Parse command-line arguments and execute background statistics.

Module Contents

BackStats.BACKDIR = 'DECam_BACK'
BackStats.calculate_one(file, xmatch_files, indir, calc_sigma_clipped=False, npix_min=100)

Calculate background statistics for one file vs multiple matches.

Computes statistics on the difference in overlap regions between one base file and a list of comparison files.

Parameters

filestr

Base FITS filename.

xmatch_fileslist of str

List of FITS filenames to compare against base file.

indirstr

Directory containing the FITS files.

calc_sigma_clippedbool, optional

If True, calculate sigma-clipped statistics. Increases runtime significantly. Default: False.

npix_minint, optional

Minimum number of overlapping pixels required. Pairs with fewer pixels are skipped. Default: 100.

Returns

list of list

List of records, where each record is a list containing:

  • file1, file2 : str - Filenames

  • npix : int - Number of overlapping pixels

  • mean, med, std : float - Unbiased statistics

  • mean_clipped, med_clipped, std_clipped : float - Sigma-clipped stats

  • skew, kurt : float - Higher moments

  • Delta : float - Mode of difference

Notes

Processing Steps:

  1. Open base file

  2. For each comparison file:

    • Find overlapping pixels (both non-zero)

    • Calculate difference in overlap region

    • Compute statistics (mean, median, std)

    • Optionally compute sigma-clipped stats

    • Calculate skewness and kurtosis

    • Calculate mode using half-sample algorithm

  3. Close files and clean up memory

Memory Management:

Uses explicit garbage collection and file closure to handle large datasets without memory leaks.

Error Handling:

  • Sets global ierror flag on failures

  • Logs errors for problematic file pairs

  • Skips pairs with insufficient overlap or NaN results

Examples

>>> records = calculate_one('base.fits', ['comp1.fits', 'comp2.fits'],
...                         'DECam_BACK/LMC_c42/T07')
>>> print(f"Processed {len(records)} comparisons")
BackStats.do_one_tile(field='LMC_c42', tile='T07', nproc=1, calc_sigma_clipped=False)

Process all overlapping file pairs for one tile.

Orchestrates background statistics calculation for all file pairs in a tile, with optional parallel processing.

Parameters

fieldstr, optional

Field name. Default: ‘LMC_c42’.

tilestr, optional

Tile identifier. Default: ‘T07’.

nprocint, optional

Number of parallel processes. If 1, uses serial processing. Default: 1.

calc_sigma_clippedbool, optional

If True, calculate sigma-clipped statistics. Default: False.

Returns

Table

Astropy Table containing background statistics for all file pairs.

Raises

IOError

If required input files or directories cannot be found.

Notes

Input Files:

  • Summary/{field}_{tile}_overlap.txt - Overlap file pairs

  • Summary/{field}_{tile}.txt - Image metadata

  • DECam_BACK/{field}/{tile}/ - Background FITS files

Processing:

  1. Read overlap file listing all file pairs

  2. Create work queue for all base files

  3. Process in parallel (if nproc > 1) or serially

  4. Combine results into single table

  5. Rescale by pixel area (0.2631² → 2² arcsec²)

  6. Join with metadata (filter, exposure time)

  7. Write to Summary directory

Memory Management:

Uses explicit garbage collection and careful table handling to avoid memory leaks in astropy join operations.

Output Scaling:

Results are scaled by (0.2631/2)² to account for pixel size difference between native and binned images.

Examples

>>> # Serial processing
>>> tab = do_one_tile('LMC_c42', 'T07')

>>> # Parallel processing with 8 threads
>>> tab = do_one_tile('LMC_c42', 'T07', nproc=8)

>>> # With sigma-clipped statistics
>>> tab = do_one_tile('LMC_c42', 'T07', nproc=8, calc_sigma_clipped=True)
BackStats.halfsamplemode(inputData, axis=None)

Compute mode of an array using the half-sample algorithm.

This is a robust estimator of the mode that works well for approximately unimodal distributions.

Parameters

inputDatandarray

Input array of values.

axisint, optional

Axis along which to compute the mode. If None, computes mode of flattened array. Default: None.

Returns

float

Mode estimate.

Notes

Algorithm:

  1. Sort the data

  2. Find the narrowest half-sample (shortest interval containing half the data)

  3. Recursively apply to the narrowest half

  4. Continue until 1-2 values remain

This algorithm is more robust than simple binning approaches and works well for skewed distributions common in background differences.

Computational Complexity:

O(n log n) due to initial sort, then O(log n) iterations.

Examples

>>> data = np.random.normal(5, 1, 1000)
>>> mode = halfsamplemode(data)
>>> print(f"Mode: {mode:.2f}")
BackStats.ierror = False
BackStats.quartile_stats(xdata)

Compute skewness and kurtosis with outlier rejection.

Uses quartile-based outlier rejection before computing higher-order moments, providing robust statistics for distributions with outliers.

Parameters

xdatandarray

Input data array.

Returns

skewfloat

Skewness (third standardized moment).

kurtfloat

Kurtosis (fourth standardized moment).

Notes

Outlier Rejection:

Uses the IQR (interquartile range) method:

  • Q1 = 25th percentile

  • Q3 = 75th percentile

  • IQR = Q3 - Q1

  • Outliers: values < Q1 - 1.5×IQR or > Q3 + 1.5×IQR

After outlier removal, computes standard moments.

Reference:

See: https://towardsdatascience.com/skewness-and-kurtosis-with-outliers-f43167532c69

Interpretation:

  • Skewness = 0: symmetric distribution

  • Skewness > 0: right-skewed (tail extends right)

  • Skewness < 0: left-skewed (tail extends left)

  • Kurtosis = 3: normal distribution

  • Kurtosis > 3: heavy tails

  • Kurtosis < 3: light tails

Examples

>>> data = np.random.normal(0, 1, 10000)
>>> skew, kurt = quartile_stats(data)
>>> print(f"Skewness: {skew:.3f}, Kurtosis: {kurt:.3f}")
BackStats.steer(argv)

Parse command-line arguments and execute background statistics.

Main entry point for command-line execution. Handles argument parsing and orchestrates tile processing.

Parameters

argvlist

Command-line arguments (typically sys.argv).

Returns

None

Results written to files and logged.

Notes

Processing:

  1. Parse arguments

  2. Expand -all to tiles T01-T16 if specified

  3. Open log file

  4. Process each tile with timing

  5. Optionally remove BackPrep files on success

  6. Close log

Error Handling:

If global ierror flag is set, does not remove BackPrep files even with -rm option.

Examples

From command line:

python BackStats.py LMC_c42 T07 T08
python BackStats.py -all -np 8 LMC_c42
python BackStats.py -rm -sigma_clipped LMC_c42 T07