Astronomical Data Processing on Alpine Linux 🌌

Processing astronomical data is like being a digital astronomer exploring the universe from your computer! 🔭 Alpine Linux provides a powerful platform for analyzing telescope images, tracking celestial objects, and discovering the wonders of space. Let’s build a complete astronomical data processing system! 🚀

Understanding Astronomical Data 🤔

Astronomical data includes:

• FITS files - Flexible Image Transport System format
• Telescope images - Raw data from observatories
• Time series data - Variable star observations
• Spectral data - Light wavelength analysis
• Catalog data - Star and galaxy databases

Think of it as decoding messages from the stars! ✨

Setting Up the Environment 📦

Install the necessary tools:

# Update package list
sudo apk update

# Install Python and development tools
sudo apk add python3 py3-pip python3-dev
sudo apk add gcc g++ musl-dev linux-headers
sudo apk add gfortran openblas-dev lapack-dev

# Install image processing tools
sudo apk add imagemagick ffmpeg
sudo apk add libjpeg-turbo-dev libpng-dev

# Install additional utilities
sudo apk add git curl wget

Installing Astronomy Libraries 🌟

Set up Python packages for astronomical analysis:

# Create virtual environment
python3 -m venv astroenv
source astroenv/bin/activate

# Install core astronomy packages
pip install numpy scipy matplotlib
pip install astropy astroquery
pip install photutils ccdproc

# Install visualization tools
pip install aplpy astroplan
pip install jupyter notebook

# Install additional analysis tools
pip install pandas scikit-image
pip install reproject regions

Working with FITS Files 📸

Process telescope images:

# FITS file processor
cat > process_fits.py << 'EOF'
#!/usr/bin/env python3
import numpy as np
from astropy.io import fits
from astropy.visualization import PercentileInterval, AsinhStretch
from astropy.visualization.mpl_normalize import ImageNormalize
import matplotlib.pyplot as plt

class FITSProcessor:
    def __init__(self, filename):
        self.filename = filename
        self.data = None
        self.header = None
        
    def load_fits(self):
        """Load FITS file data"""
        with fits.open(self.filename) as hdul:
            self.data = hdul[0].data
            self.header = hdul[0].header
            
        print(f"Loaded FITS file: {self.filename}")
        print(f"Image dimensions: {self.data.shape}")
        print(f"Object: {self.header.get('OBJECT', 'Unknown')}")
        print(f"Telescope: {self.header.get('TELESCOP', 'Unknown')}")
        
    def display_image(self, save_path=None):
        """Display FITS image with proper scaling"""
        # Create figure
        fig, ax = plt.subplots(figsize=(10, 10))
        
        # Apply stretch and normalization
        transform = AsinhStretch() + PercentileInterval(99.5)
        norm = ImageNormalize(self.data, transform)
        
        # Display image
        im = ax.imshow(self.data, cmap='gray', norm=norm, origin='lower')
        
        # Add colorbar
        plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
        
        # Add title
        obj_name = self.header.get('OBJECT', 'Unknown Object')
        ax.set_title(f'{obj_name} - {self.filename}', fontsize=16)
        
        # Remove axes for cleaner look
        ax.set_xticks([])
        ax.set_yticks([])
        
        if save_path:
            plt.savefig(save_path, dpi=150, bbox_inches='tight')
            print(f"Saved image to {save_path}")
        
        plt.show()
        
    def calculate_statistics(self):
        """Calculate image statistics"""
        stats = {
            'mean': np.mean(self.data),
            'median': np.median(self.data),
            'std': np.std(self.data),
            'min': np.min(self.data),
            'max': np.max(self.data),
            'total_flux': np.sum(self.data)
        }
        
        print("\nImage Statistics:")
        for key, value in stats.items():
            print(f"  {key}: {value:.2f}")
        
        return stats

# Example usage
if __name__ == "__main__":
    # Process a FITS file
    processor = FITSProcessor("galaxy.fits")
    processor.load_fits()
    processor.calculate_statistics()
    processor.display_image("galaxy_processed.png")
EOF

chmod +x process_fits.py

Star Detection and Photometry ⭐

Find and measure stars in images:

# Star detection script
cat > star_detection.py << 'EOF'
#!/usr/bin/env python3
import numpy as np
from astropy.io import fits
from astropy.stats import sigma_clipped_stats
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture
import matplotlib.pyplot as plt

class StarDetector:
    def __init__(self, fits_file):
        self.fits_file = fits_file
        self.data = None
        self.stars = None
        
    def load_data(self):
        """Load FITS data"""
        with fits.open(self.fits_file) as hdul:
            self.data = hdul[0].data.astype(float)
            
    def find_stars(self, fwhm=3.0, threshold_sigma=5.0):
        """Detect stars in the image"""
        # Calculate background statistics
        mean, median, std = sigma_clipped_stats(self.data, sigma=3.0)
        
        print(f"Background statistics:")
        print(f"  Mean: {mean:.2f}")
        print(f"  Median: {median:.2f}")
        print(f"  Std Dev: {std:.2f}")
        
        # Find stars
        daofind = DAOStarFinder(fwhm=fwhm, threshold=threshold_sigma*std)
        self.stars = daofind(self.data - median)
        
        print(f"\nFound {len(self.stars)} stars")
        
        # Show brightest stars
        self.stars.sort('peak')
        self.stars.reverse()
        
        print("\nBrightest 10 stars:")
        print("  ID   X       Y      Peak")
        for star in self.stars[:10]:
            print(f"  {star['id']:3d} {star['xcentroid']:7.2f} {star['ycentroid']:7.2f} {star['peak']:8.1f}")
            
    def plot_stars(self, save_path=None):
        """Plot detected stars on image"""
        fig, ax = plt.subplots(figsize=(12, 10))
        
        # Display image
        mean, median, std = sigma_clipped_stats(self.data, sigma=3.0)
        im = ax.imshow(self.data, cmap='gray_r', 
                      vmin=median - 3*std, 
                      vmax=median + 10*std,
                      origin='lower')
        
        # Mark detected stars
        positions = np.transpose((self.stars['xcentroid'], 
                                 self.stars['ycentroid']))
        apertures = CircularAperture(positions, r=10)
        apertures.plot(color='red', lw=2, alpha=0.7)
        
        # Add star IDs for brightest stars
        for i, star in enumerate(self.stars[:20]):
            ax.text(star['xcentroid'] + 15, star['ycentroid'], 
                   f"{star['id']}", 
                   color='yellow', fontsize=10)
        
        ax.set_title(f'Detected Stars: {len(self.stars)}', fontsize=16)
        plt.colorbar(im, ax=ax, fraction=0.046)
        
        if save_path:
            plt.savefig(save_path, dpi=150, bbox_inches='tight')
            
        plt.show()
        
    def perform_photometry(self):
        """Measure star brightness"""
        from photutils.aperture import aperture_photometry
        
        # Define apertures
        positions = np.transpose((self.stars['xcentroid'], 
                                 self.stars['ycentroid']))
        apertures = CircularAperture(positions, r=5)
        
        # Perform photometry
        phot_table = aperture_photometry(self.data, apertures)
        
        # Add to stars table
        self.stars['flux'] = phot_table['aperture_sum']
        
        # Calculate magnitudes (instrumental)
        self.stars['mag_inst'] = -2.5 * np.log10(self.stars['flux'])
        
        return self.stars

# Example usage
if __name__ == "__main__":
    detector = StarDetector("star_field.fits")
    detector.load_data()
    detector.find_stars()
    detector.plot_stars("detected_stars.png")
    stars = detector.perform_photometry()
EOF

chmod +x star_detection.py

Celestial Coordinate Systems 🧭

Convert between coordinate systems:

# Coordinate converter
cat > coordinates.py << 'EOF'
#!/usr/bin/env python3
from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time
import astropy.units as u
import numpy as np

class CelestialCoordinates:
    def __init__(self, observer_location=None):
        if observer_location:
            self.location = EarthLocation(
                lat=observer_location['lat']*u.deg,
                lon=observer_location['lon']*u.deg,
                height=observer_location['height']*u.m
            )
        else:
            # Default: Greenwich Observatory
            self.location = EarthLocation.of_site('Greenwich')
            
    def ra_dec_to_altaz(self, ra, dec, obstime=None):
        """Convert RA/Dec to Alt/Az"""
        if obstime is None:
            obstime = Time.now()
        
        # Create coordinate object
        coord = SkyCoord(ra=ra*u.deg, dec=dec*u.deg, frame='icrs')
        
        # Convert to horizontal coordinates
        altaz = coord.transform_to(AltAz(obstime=obstime, location=self.location))
        
        return altaz.alt.deg, altaz.az.deg
    
    def find_best_observation_time(self, target, date=None):
        """Find when object is highest in sky"""
        if date is None:
            date = Time.now()
        
        # Create time array for 24 hours
        times = date + np.linspace(0, 24, 144)*u.hour
        
        # Calculate altitude for each time
        altitudes = []
        for time in times:
            alt, _ = self.ra_dec_to_altaz(target['ra'], target['dec'], time)
            altitudes.append(alt)
        
        # Find maximum altitude
        max_idx = np.argmax(altitudes)
        best_time = times[max_idx]
        max_alt = altitudes[max_idx]
        
        return best_time, max_alt
    
    def calculate_airmass(self, altitude):
        """Calculate airmass for given altitude"""
        if altitude <= 0:
            return np.inf
        
        zenith_angle = 90 - altitude
        airmass = 1 / np.cos(np.radians(zenith_angle))
        
        return airmass

# Example usage
if __name__ == "__main__":
    # Create coordinate converter
    converter = CelestialCoordinates({
        'lat': 37.7749,  # San Francisco
        'lon': -122.4194,
        'height': 100
    })
    
    # Convert coordinates for M31 (Andromeda Galaxy)
    m31 = {'ra': 10.6847, 'dec': 41.2687}
    alt, az = converter.ra_dec_to_altaz(m31['ra'], m31['dec'])
    
    print(f"M31 Andromeda Galaxy:")
    print(f"  Current Alt/Az: {alt:.1f}°, {az:.1f}°")
    print(f"  Airmass: {converter.calculate_airmass(alt):.2f}")
    
    # Find best observation time
    best_time, max_alt = converter.find_best_observation_time(m31)
    print(f"  Best observation: {best_time.iso}")
    print(f"  Maximum altitude: {max_alt:.1f}°")
EOF

chmod +x coordinates.py

Time Series Analysis 📈

Analyze variable stars:

# Variable star analysis
cat > variable_stars.py << 'EOF'
#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt
from astropy.timeseries import LombScargle
from astropy.time import Time
import pandas as pd

class VariableStarAnalyzer:
    def __init__(self):
        self.data = None
        self.periods = None
        
    def load_observations(self, filename):
        """Load time series data"""
        # Expected format: JD, magnitude, error
        self.data = pd.read_csv(filename, names=['jd', 'mag', 'err'])
        
        print(f"Loaded {len(self.data)} observations")
        print(f"Time span: {self.data['jd'].max() - self.data['jd'].min():.1f} days")
        
    def plot_light_curve(self):
        """Plot the light curve"""
        plt.figure(figsize=(12, 6))
        
        plt.errorbar(self.data['jd'], self.data['mag'], 
                    yerr=self.data['err'], 
                    fmt='o', markersize=4, alpha=0.7)
        
        plt.xlabel('Julian Date')
        plt.ylabel('Magnitude')
        plt.title('Variable Star Light Curve')
        plt.gca().invert_yaxis()  # Brighter = lower magnitude
        plt.grid(True, alpha=0.3)
        plt.show()
        
    def find_period(self, min_period=0.1, max_period=100):
        """Find periodicity using Lomb-Scargle"""
        # Prepare data
        t = self.data['jd'].values
        y = self.data['mag'].values
        dy = self.data['err'].values
        
        # Create frequency grid
        frequency = np.linspace(1/max_period, 1/min_period, 10000)
        
        # Compute periodogram
        ls = LombScargle(t, y, dy)
        power = ls.power(frequency)
        
        # Find best period
        best_frequency = frequency[np.argmax(power)]
        best_period = 1/best_frequency
        
        # Calculate false alarm probability
        fap = ls.false_alarm_probability(power.max())
        
        print(f"\nPeriod Analysis:")
        print(f"  Best period: {best_period:.4f} days")
        print(f"  False alarm probability: {fap:.2e}")
        
        # Plot periodogram
        plt.figure(figsize=(10, 6))
        periods = 1/frequency
        plt.semilogx(periods, power)
        plt.axvline(best_period, color='red', linestyle='--', 
                   label=f'Best period: {best_period:.4f} days')
        plt.xlabel('Period (days)')
        plt.ylabel('Lomb-Scargle Power')
        plt.title('Period Analysis')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.show()
        
        self.best_period = best_period
        return best_period
        
    def phase_fold(self, period=None):
        """Create phase-folded light curve"""
        if period is None:
            period = self.best_period
            
        # Calculate phases
        phases = ((self.data['jd'] - self.data['jd'].min()) % period) / period
        
        # Plot phase-folded curve
        plt.figure(figsize=(10, 6))
        
        # Plot data twice for clarity
        for offset in [0, 1]:
            plt.errorbar(phases + offset, self.data['mag'], 
                        yerr=self.data['err'],
                        fmt='o', markersize=4, alpha=0.6)
        
        plt.xlim(0, 2)
        plt.xlabel('Phase')
        plt.ylabel('Magnitude')
        plt.title(f'Phase-folded Light Curve (P = {period:.4f} days)')
        plt.gca().invert_yaxis()
        plt.grid(True, alpha=0.3)
        plt.show()

# Example usage
if __name__ == "__main__":
    analyzer = VariableStarAnalyzer()
    
    # Create sample data
    # In practice, load real observations
    t = np.linspace(0, 100, 500) + np.random.normal(0, 0.1, 500)
    period = 2.3456
    mag = 12.0 + 0.5 * np.sin(2*np.pi*t/period) + np.random.normal(0, 0.02, 500)
    err = np.full_like(mag, 0.02)
    
    # Save sample data
    sample_data = pd.DataFrame({'jd': t, 'mag': mag, 'err': err})
    sample_data.to_csv('variable_star.csv', index=False, header=False)
    
    # Analyze
    analyzer.load_observations('variable_star.csv')
    analyzer.plot_light_curve()
    analyzer.find_period()
    analyzer.phase_fold()
EOF

chmod +x variable_stars.py

Sky Visualization 🌠

Create sky maps and charts:

# Sky map generator
cat > sky_map.py << 'EOF'
#!/usr/bin/env python3
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
from astropy.time import Time
from astroplan import Observer
from astroplan.plots import plot_sky
import astropy.units as u

class SkyMapper:
    def __init__(self, observer_location):
        self.observer = Observer(
            location=EarthLocation(
                lat=observer_location['lat']*u.deg,
                lon=observer_location['lon']*u.deg,
                height=observer_location['height']*u.m
            ),
            name=observer_location.get('name', 'Observer')
        )
        
    def create_sky_chart(self, time=None, targets=None):
        """Create current sky chart"""
        if time is None:
            time = Time.now()
            
        # Create figure
        fig, ax = plt.subplots(figsize=(10, 10), 
                              subplot_kw=dict(projection='polar'))
        
        # Plot sky
        plot_sky(targets, self.observer, time, ax=ax)
        
        # Customize
        ax.set_title(f"Sky Chart - {self.observer.name}\n{time.iso}", 
                    fontsize=16, pad=20)
        
        plt.tight_layout()
        return fig
        
    def plot_object_visibility(self, target, date=None):
        """Plot object visibility over night"""
        from astroplan import plots
        
        if date is None:
            date = Time.now()
            
        # Create figure
        fig, ax = plt.subplots(figsize=(12, 6))
        
        # Plot airmass
        plots.plot_airmass(target, self.observer, date, ax=ax)
        
        ax.set_title(f"Visibility of {target.name}")
        plt.tight_layout()
        
        return fig

# Example usage
if __name__ == "__main__":
    # Create sky mapper
    mapper = SkyMapper({
        'name': 'Alpine Observatory',
        'lat': 46.5197,
        'lon': 6.6323,
        'height': 372
    })
    
    # Define targets
    m42 = SkyCoord.from_name('M42')  # Orion Nebula
    m31 = SkyCoord.from_name('M31')  # Andromeda
    
    # Create charts
    mapper.create_sky_chart(targets=[m42, m31])
    plt.savefig('current_sky.png')
EOF

Data Pipeline Automation 🔄

Automate processing workflow:

# Automated pipeline
cat > astro_pipeline.sh << 'EOF'
#!/bin/sh
# Astronomical Data Processing Pipeline

DATA_DIR="$HOME/astrodata"
PROCESSED_DIR="$DATA_DIR/processed"
LOG_FILE="$DATA_DIR/pipeline.log"

# Create directories
mkdir -p "$DATA_DIR/raw" "$PROCESSED_DIR" "$DATA_DIR/calibration"

log() {
    echo "[$(date '+%Y-%m-%d %H:%M:%S')] $1" | tee -a "$LOG_FILE"
}

# Download new data
download_data() {
    log "Downloading new observations..."
    # Add your data source here
    # wget/curl commands
}

# Process FITS files
process_fits() {
    log "Processing FITS files..."
    
    for fits_file in "$DATA_DIR/raw"/*.fits; do
        if [ -f "$fits_file" ]; then
            basename=$(basename "$fits_file" .fits)
            log "Processing $basename..."
            
            python process_fits.py "$fits_file" \
                --output "$PROCESSED_DIR/${basename}_processed.fits"
        fi
    done
}

# Run star detection
detect_stars() {
    log "Running star detection..."
    
    for fits_file in "$PROCESSED_DIR"/*_processed.fits; do
        if [ -f "$fits_file" ]; then
            python star_detection.py "$fits_file"
        fi
    done
}

# Generate reports
generate_reports() {
    log "Generating reports..."
    
    # Create summary report
    python generate_astro_report.py \
        --input "$PROCESSED_DIR" \
        --output "$DATA_DIR/reports/$(date +%Y%m%d)_report.pdf"
}

# Main pipeline
main() {
    log "Starting astronomical data pipeline"
    
    download_data
    process_fits
    detect_stars
    generate_reports
    
    log "Pipeline completed successfully"
}

# Run pipeline
main
EOF

chmod +x astro_pipeline.sh

Best Practices 📌

1. Calibrate images - Use dark frames and flats
2. Handle large files - Use memory mapping for big FITS
3. Document observations - Keep detailed logs
4. Version control - Track analysis code changes
5. Validate results - Cross-check with known sources

Troubleshooting 🔧

Memory Issues with Large FITS

# Use memory mapping
from astropy.io import fits
hdul = fits.open('large_file.fits', memmap=True)

Coordinate Conversion Errors

# Always specify frames explicitly
coord = SkyCoord(ra=10*u.deg, dec=20*u.deg, frame='icrs')

Quick Reference 📋

# Activate environment
source astroenv/bin/activate

# Process FITS file
python process_fits.py image.fits

# Detect stars
python star_detection.py image.fits

# Analyze variable star
python variable_stars.py lightcurve.csv

# Check coordinates
python coordinates.py

Conclusion 🎯

You’ve built a comprehensive astronomical data processing system on Alpine Linux! From processing telescope images to analyzing variable stars and creating sky maps, you now have the tools to explore the universe digitally. Remember to always validate your results and enjoy your journey through the cosmos! 🌌✨