Image Processing

Image Processing#

Overview#

The Image class extends Parameter with astronomical-specific functionality, providing comprehensive support for astronomical image handling, FITS header management, coordinate systems, and unit conversions.

Note

The Image class automatically parses FITS headers, manages coordinate systems, and provides astronomical-specific properties and methods.

Key Features#

  • Coordinate Systems: Automatic parsing of FITS headers for coordinate information

  • Phase Center Management: Automatic detection and management of phase center coordinates

  • FITS Header Support: Comprehensive FITS header creation and parsing

  • Unit Conversions: Built-in intensity unit transformations

  • Noise Estimation: Robust noise calculation using sigma-clipped statistics

Creating Images#

Initialize images from various data sources:

from pyralysis.reconstruction import Image
import numpy as np
import astropy.units as u

# From numpy array with cellsize
data = np.random.random((512, 512))
image = Image(
    data=data,
    cellsize=0.1*u.arcsec
)

# From existing parameter
from pyralysis.reconstruction import Parameter
param = Parameter(data=data, cellsize=0.1*u.arcsec)
image = Image(data=param.data, cellsize=param.cellsize)

# From FITS file with automatic header parsing
from pyralysis.io import FITS
fits_io = FITS(input_name="image.fits", use_dask=True)
image = fits_io.read()  # Automatically creates Image with parsed header

# From Zarr array with Dask support
from pyralysis.io import ZarrArray
zarr_io = ZarrArray(input_name="image.zarr", use_dask=True)
image = zarr_io.read()  # Automatically creates Image with parsed attributes

# From Measurement Set (for visibility data)
from pyralysis.io import DaskMS
ms_io = DaskMS(input_name="data.ms")
dataset = ms_io.read()  # Returns Dataset object, not Image

Image I/O Operations#

Pyralysis provides multiple I/O options for reading and writing images and datasets:

FITS Files (Recommended for Astronomical Data):

from pyralysis.io import FITS

# Read FITS image with automatic header parsing
fits_io = FITS(input_name="image.fits", use_dask=True)
image = fits_io.read()

# Write image to FITS with header preservation
fits_io = FITS(output_name="output.fits")
fits_io.write(image, overwrite=True)

# Read with custom chunking for large files
fits_io = FITS(
    input_name="large_image.fits",
    chunks=(512, 512),
    use_dask=True
)
image = fits_io.read()

Zarr Arrays (Efficient for Large Datasets):

from pyralysis.io import ZarrArray

# Read Zarr image with Dask support
zarr_io = ZarrArray(input_name="image.zarr", use_dask=True)
image = zarr_io.read()

# Write image to Zarr with compression
zarr_io = ZarrArray(output_name="output.zarr")
zarr_io.write(image, overwrite=True)

# Read with custom chunking
zarr_io = ZarrArray(
    input_name="large_image.zarr",
    chunks=(256, 256),
    use_dask=True
)
image = zarr_io.read()

Measurement Sets (Radio Interferometry Data):

from pyralysis.io import DaskMS

# Read Measurement Set with Dask (default row chunk size: 100_000)
ms_io = DaskMS(input_name="data.ms")
dataset = ms_io.read()

# Read with custom chunking strategy
ms_io = DaskMS(
    input_name="data.ms",
    chunks={
        'row': (40000, 60000, 40000, 60000),
        'chan': (16, 16, 16, 16),
        'corr': (1, 2, 1)
    }
)
dataset = ms_io.read()

I/O Format Comparison:

Feature

FITS

Zarr

Measurement Set

Best For

Astronomical images

Large datasets and compression

Radio astronomy visibilities

Header Support

FITS headers (automatic)

Zarr attributes (automatic)

MS metadata (automatic)

Chunking

Dask support

Native chunks and compression

Advanced chunks and parallel I/O

Performance

Good (memory-mapped)

Excellent (parallel I/O)

Excellent (Dask-based)

File Size

Standard

Compressed

Large

Recommendations: - FITS: Best for standard astronomical images with headers - Zarr: Best for large datasets requiring compression and parallel I/O - Measurement Set: Best for radio interferometry visibility data

Best Practices for I/O Operations:

# 1. Choose appropriate format for your data
if image_size < 1024*1024:
    # Small images: Use FITS for compatibility
    io_handler = FITS(output_name="small_image.fits")
else:
    # Large images: Use Zarr for performance
    io_handler = ZarrArray(output_name="large_image.zarr", use_dask=True)

# 2. Configure chunking for large datasets
if image_size > 2048*2048:
    io_handler.chunks = (256, 256)  # Optimal chunking

# 3. Enable Dask for out-of-core processing
io_handler = FITS(
    input_name="huge_image.fits",
    use_dask=True,
    chunks=(512, 512)
)

# 4. Preserve metadata and attributes
io_handler.preserve_attributes = True
io_handler.auto_parse_header = True

Automatic Header Parsing#

The Image class automatically extracts information from FITS headers:

# Create image from FITS data
image = Image(data=fits_data_with_header)

# Access parsed information
print(f"Phase center: {image.phase_center}")
print(f"Center pixel: {image.center_pixel}")
print(f"Cellsize: {image.cellsize}")

# Check if header was parsed
if image.attrs:
    print("FITS header was successfully parsed")

Coordinate System Properties#

Access astronomical coordinate information:

# Phase center coordinates
if image.phase_center:
    ra = image.phase_center.ra
    dec = image.phase_center.dec
    frame = image.phase_center.frame
    print(f"Phase center: {ra}, {dec} ({frame})")

# Center pixel coordinates
if image.center_pixel is not None:
    center_y, center_x = image.center_pixel
    print(f"Center pixel: ({center_y}, {center_x})")

# Image dimensions
if image.imsize:
    height, width = image.imsize
    print(f"Image size: {height} x {width} pixels")

Spatial Properties#

Calculate and access spatial information:

# Field of view
if image.field_of_view:
    fov_y, fov_x = image.field_of_view
    print(f"Field of view: {fov_y} x {fov_x}")

# Pixel coordinates
if image.pixel_coordinates:
    y_coords, x_coords = image.pixel_coordinates
    print(f"Pixel coordinate grids: {y_coords.shape} x {x_coords.shape}")

# Sky coordinates
if image.sky_coordinates:
    sky_coords = image.sky_coordinates
    print(f"Sky coordinate grid: {sky_coords.shape}")

FITS Header Management#

Create and manage FITS headers:

# Create basic header
image.create_header(
    projection="SIN",
    bunit="JY/BEAM",
    dataset=dataset,
    add_beam=True
)

# Access header
header = image.header
if header:
    print(f"Header has {len(header)} keywords")

# Update header attributes
image.attrs["BUNIT"] = "JY/PIXEL"
image.attrs["OBJECT"] = "Test Source"

Unit Conversions#

Convert between different intensity unit systems:

# Convert intensity units
flux_jy = image.transform_intensity_units(
    intensity=1.0*u.Jy/u.beam,
    beam_area=beam_area,
    unit=u.pixel
)

print(f"Flux in Jy/pixel: {flux_jy}")

# Convert to different units
flux_mjy = flux_jy.to(u.mJy/u.pixel)
print(f"Flux in mJy/pixel: {flux_mjy}")

Noise Estimation#

Calculate noise estimates using robust statistics:

# Calculate noise using sigma-clipped MAD
noise = image.calculate_noise(sigma=3.0, maxiters=5)
print(f"Noise estimate: {noise}")

# Custom noise calculation parameters
noise_robust = image.calculate_noise(sigma=5.0, maxiters=10)
print(f"Robust noise estimate: {noise_robust}")

Advanced Image Features#

  • Multi-dimensional Support: Handle spectral and polarization data

  • Automatic Validation: Validate cellsize and coordinate information

  • Header Caching: Efficient FITS header management and updates

  • Coordinate Transformations: Automatic coordinate system handling

Multi-dimensional Support#

Handle images with spectral and polarization dimensions:

# 3D spectral image
spectral_data = np.random.random((512, 512, 50))
spectral_image = Image(
    data=spectral_data,
    cellsize=0.1*u.arcsec
)

# Access spectral properties
if spectral_image.ndim >= 3:
    n_channels = spectral_image.shape[-3]
    print(f"Number of channels: {n_channels}")

# 4D polarization image
pol_data = np.random.random((512, 512, 50, 4))
pol_image = Image(
    data=pol_data,
    cellsize=0.1*u.arcsec
)

Coordinate Validation#

Automatic validation of coordinate information:

# Validate cellsize
if image.cellsize is None:
    print("Warning: Cellsize not set")
else:
    # Check units
    if not check_units(image.cellsize, u.rad):
        print("Warning: Cellsize should have angular units")

# Validate phase center
if image.phase_center is None:
    print("Warning: Phase center not available")

# Validate center pixel
if image.center_pixel is None:
    print("Warning: Center pixel not set")

Performance Optimizations#

Optimize image processing performance:

# Use appropriate chunking
image.chunks = (256, 256)  # Good for most astronomical images

# Enable persistence for repeated access
image.persist()

# Compute when needed
if image.is_dask():
    image.compute()  # Convert to numpy for faster access

Best Practices#

  1. Set Cellsize: Always provide cellsize for proper astronomical calculations

  2. Validate Headers: Check that FITS headers are properly parsed

  3. Use Appropriate Chunking: Configure chunking based on image size and memory constraints

  4. Monitor Memory: Track memory usage for large images

  5. Validate Coordinates: Ensure coordinate systems are properly configured

Example: Complete Image Workflow

# Create image from FITS data
image = Image(data=fits_data_with_header)

# Validate and set properties
if image.cellsize is None:
    image.cellsize = [0.1, 0.1] * u.arcsec

# Create comprehensive header
image.create_header(
    projection="SIN",
    bunit="JY/BEAM",
    dataset=dataset,
    add_beam=True
)

# Calculate noise
noise = image.calculate_noise()
print(f"Image noise: {noise}")

# Access astronomical properties
print(f"Field of view: {image.field_of_view}")
print(f"Phase center: {image.phase_center}")

# Convert units
flux_pixel = image.transform_intensity_units(
    intensity=1.0*u.Jy/u.beam,
    beam_area=beam_area,
    unit=u.pixel
)

Troubleshooting#

Common issues and solutions:

# Issue: Cellsize not set
if image.cellsize is None:
    print("Set cellsize: image.cellsize = [0.1, 0.1] * u.arcsec")

# Issue: Phase center not available
if image.phase_center is None:
    print("Check FITS header for CRVAL1, CRVAL2 keywords")

# Issue: Center pixel not set
if image.center_pixel is None:
    print("Check FITS header for CRPIX1, CRPIX2 keywords")

# Issue: Header not parsed
if not image.attrs:
    print("Ensure data contains FITS header information")

I/O Troubleshooting:

# Issue: FITS file too large for memory
# Solution: Use Dask chunking
fits_io = FITS(
    input_name="large_image.fits",
    use_dask=True,
    chunks=(256, 256)
)

# Issue: Zarr file not found
# Solution: Check if directory exists
if not zarr_io.exists():
    print("Zarr directory does not exist")

# Issue: Measurement Set chunking too large
# Solution: Reduce chunk size
ms_io = DaskMS(
    input_name="data.ms",
    chunks={'row': 1000}  # Smaller chunks
)

Parameter Management | I/O Operations & Data Formats