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: .. code-block:: python 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):** .. code-block:: python 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):** .. code-block:: python 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):** .. code-block:: python 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 | Large datasets | Radio astronomy | | | images | and compression | visibilities | +-----------------+------------------+------------------+------------------+ | **Header** | FITS headers | Zarr attributes | MS metadata | | **Support** | (automatic) | (automatic) | (automatic) | +-----------------+------------------+------------------+------------------+ | **Chunking** | Dask support | Native chunks | Advanced chunks | | | | and compression | and parallel I/O | +-----------------+------------------+------------------+------------------+ | **Performance** | Good | Excellent | Excellent | | | (memory-mapped) | (parallel I/O) | (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:** .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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: .. code-block:: python # 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** .. code-block:: python # 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: .. code-block:: python # 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:** .. code-block:: python # 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 ) ---- :doc:`parameter_management` | :doc:`io_operations`