Atmospheric Modeling

Overview

The atmospheric modeling system is the foundation of accurate balloon trajectory prediction. It provides comprehensive atmospheric state data including temperature, pressure, density, winds, and humidity from the surface to the edge of space. The system employs a sophisticated hierarchy of data sources, from real-time weather APIs to validated atmospheric models, ensuring global coverage with optimal accuracy.

Data Source Hierarchy

The simulator implements a sophisticated fallback system to ensure atmospheric data is always available:

Data Source Selection Algorithm

def select_atmospheric_data_source(lat, lon, alt, time):
    """
    Intelligent data source selection based on availability,
    location, and quality metrics
    """
    # Try primary real-time sources
    if openweather_api.is_available():
        data = openweather_api.fetch(lat, lon, time)
        if data.quality_score > 0.8:
            return data, "OpenWeather API"

    # US-specific high-quality data
    if is_within_us(lat, lon) and noaa_api.is_available():
        data = noaa_api.fetch(lat, lon, time)
        if data.quality_score > 0.9:
            return data, "NOAA Weather API"

    # Global forecast model
    if gfs_nomads.is_available():
        data = gfs_nomads.fetch(lat, lon, time)
        if data.quality_score > 0.7:
            return data, "GFS NOMADS"

    # Fallback to standard atmosphere
    return standard_atmosphere.calculate(lat, lon, alt), "US Standard Atmosphere"

US Standard Atmosphere 1976

The US Standard Atmosphere 1976 provides a mathematical model of Earth's atmosphere based on average conditions at mid-latitudes. It defines temperature, pressure, and density as functions of altitude through a series of atmospheric layers.

Temperature Profile by Layer

The atmosphere is divided into seven distinct layers, each with its own temperature gradient:

Pressure Calculation

Pressure is calculated using the barometric formula, which differs for isothermal and gradient layers:

For Isothermal Layers (constant temperature):

For Gradient Layers (linear temperature variation):

Density Calculation

Air density is calculated from the ideal gas law:

Implementation Example

class USStandardAtmosphere1976:
    def __init__(self):
        # Layer definitions: (h_base[km], T_base[K], lapse_rate[K/km])
        self.layers = [
            (0,  288.15, -6.5),   # Troposphere
            (11, 216.65,  0.0),   # Tropopause
            (20, 216.65,  1.0),   # Lower Stratosphere
            (32, 228.65,  2.8),   # Upper Stratosphere
            (47, 270.65,  0.0),   # Stratopause
            (51, 270.65, -2.8),   # Lower Mesosphere
            (71, 214.65, -2.0),   # Upper Mesosphere
        ]

        # Physical constants
        self.g0 = 9.80665      # m/s²
        self.M = 0.0289644     # kg/mol
        self.R = 8.3144598     # J/(mol·K)
        self.P0 = 101325       # Pa (sea level)

    def get_temperature(self, altitude_km):
        """Calculate temperature at given altitude"""
        for i, (h_base, T_base, lapse_rate) in enumerate(self.layers):
            if i == len(self.layers) - 1 or altitude_km < self.layers[i+1][0]:
                return T_base + lapse_rate * (altitude_km - h_base)

        # Above 86 km
        return self.layers[-1][1]

    def get_pressure(self, altitude_m):
        """Calculate pressure at given altitude"""
        h_km = altitude_m / 1000.0

        # Find current layer
        layer_idx = 0
        for i in range(len(self.layers) - 1):
            if h_km >= self.layers[i+1][0]:
                layer_idx = i + 1
            else:
                break

        # Calculate pressure from sea level up through each layer
        P = self.P0
        h_current = 0

        for i in range(layer_idx + 1):
            h_base, T_base, lapse_rate = self.layers[i]
            h_top = self.layers[i+1][0] if i < len(self.layers)-1 else h_km
            h_top = min(h_top, h_km)

            if abs(lapse_rate) < 1e-10:  # Isothermal
                P *= np.exp(-self.g0 * self.M * (h_top - h_base) * 1000 / (self.R * T_base))
            else:  # Gradient layer
                T_top = T_base + lapse_rate * (h_top - h_base)
                P *= (T_top / T_base) ** (-self.g0 * self.M / (self.R * lapse_rate * 0.001))

            if h_top >= h_km:
                break

        return P

Real-Time Data Integration

OpenWeather API Integration

OpenWeather provides global atmospheric data with 3-hour forecasts and current conditions:

class OpenWeatherAPI:
    def __init__(self, api_key):
        self.api_key = api_key
        self.base_url = "https://api.openweathermap.org/data/2.5"

    def fetch_atmospheric_profile(self, lat, lon, time=None):
        """Fetch complete atmospheric profile from OpenWeather"""

        # Get current conditions
        current = self._fetch_current(lat, lon)

        # Get forecast data for vertical interpolation
        forecast = self._fetch_forecast(lat, lon)

        # Build atmospheric profile
        profile = AtmosphericProfile()

        # Surface conditions
        profile.add_level(
            altitude=current['main']['grnd_level'] if 'grnd_level' in current['main'] else 0,
            temperature=current['main']['temp'],
            pressure=current['main']['pressure'] * 100,  # hPa to Pa
            humidity=current['main']['humidity'],
            wind_speed=current['wind']['speed'],
            wind_direction=current['wind']['deg']
        )

        # Interpolate to standard pressure levels
        return self._interpolate_profile(profile, forecast)

NOAA Weather API Integration

NOAA provides high-resolution data for US locations through the National Weather Service API:

Data Quality Metrics

Each data source is evaluated for quality using multiple metrics:

GFS NOMADS Processing

The Global Forecast System (GFS) provides comprehensive atmospheric data in GRIB2 format through NOAA's NOMADS servers.

Pressure Levels

GFS data is provided at 26 standard pressure levels:

Variables Extracted

Pressure to Altitude Conversion

Converting pressure levels to geometric altitude uses the barometric formula:

GRIB2 Processing Implementation

class GFSNOMADSFetcher:
    def __init__(self):
        self.base_url = "https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl"
        self.pressure_levels = [1000, 975, 950, 925, 900, 850, 800, 750, 700,
                               650, 600, 550, 500, 450, 400, 350, 300, 250,
                               200, 150, 100, 70, 50, 30, 20, 10]

    def fetch_atmospheric_data(self, lat, lon, forecast_hour=0):
        """Fetch GFS data for specific location"""

        # Download GRIB2 file for the region
        grib_file = self._download_grib_subset(lat, lon, forecast_hour)

        # Parse GRIB2 data
        atmospheric_data = {}

        with pygrib.open(grib_file) as grbs:
            for grb in grbs:
                if grb.level in self.pressure_levels:
                    level_hPa = grb.level

                    if grb.name == 'Temperature':
                        atmospheric_data.setdefault(level_hPa, {})['temperature'] = grb.values
                    elif grb.name == 'U component of wind':
                        atmospheric_data.setdefault(level_hPa, {})['u_wind'] = grb.values
                    elif grb.name == 'V component of wind':
                        atmospheric_data.setdefault(level_hPa, {})['v_wind'] = grb.values
                    # ... process other variables

        # Convert to altitude-based profile
        return self._convert_to_altitude_profile(atmospheric_data, lat, lon)

Derived Parameters

From the basic GFS variables, we calculate additional parameters:

Data Fusion Algorithms

When multiple data sources are available, sophisticated fusion algorithms combine them optimally:

Weighted Average Fusion

For overlapping data sources, we use quality-weighted averaging:

where weights are determined by:

Vertical Interpolation

Atmospheric data often comes at discrete levels. We use cubic spline interpolation in log-pressure coordinates:

Gap Filling

When data is missing at certain levels, physics-based interpolation is used:

def fill_atmospheric_gaps(profile):
    """Fill gaps using physical constraints"""

    # Temperature: use moist adiabatic lapse rate if unstable
    if is_unstable(profile):
        lapse_rate = calculate_moist_adiabatic_lapse_rate(
            profile.temperature,
            profile.pressure,
            profile.humidity
        )
    else:
        lapse_rate = -6.5  # K/km standard

    # Fill temperature gaps
    profile.temperature = interpolate_with_lapse_rate(
        profile.altitude,
        profile.temperature,
        lapse_rate
    )

    # Pressure: use hydrostatic equation
    for i in range(len(profile.altitude) - 1):
        if np.isnan(profile.pressure[i+1]):
            dz = profile.altitude[i+1] - profile.altitude[i]
            T_avg = (profile.temperature[i] + profile.temperature[i+1]) / 2
            profile.pressure[i+1] = profile.pressure[i] * np.exp(
                -g * M_air * dz / (R * T_avg)
            )

Quality Control

Physical Consistency Checks

All atmospheric data undergoes rigorous quality control:

Statistical Quality Metrics

Statistical measures identify outliers and anomalies:

def calculate_quality_metrics(profile):
    """Calculate comprehensive quality metrics"""

    metrics = {}

    # Vertical resolution quality
    dz = np.diff(profile.altitude)
    metrics['vertical_resolution'] = 1.0 / (1.0 + np.std(dz) / np.mean(dz))

    # Temperature smoothness
    d2T_dz2 = np.diff(np.diff(profile.temperature))
    metrics['temperature_smoothness'] = 1.0 / (1.0 + np.std(d2T_dz2))

    # Data completeness
    metrics['completeness'] = np.sum(~np.isnan(profile.temperature)) / len(profile.temperature)

    # Physical consistency
    violations = 0
    violations += check_hydrostatic_balance(profile)
    violations += check_temperature_limits(profile)
    violations += check_wind_shear(profile)
    metrics['physical_consistency'] = 1.0 / (1.0 + violations)

    # Overall quality score
    metrics['overall'] = np.mean([
        metrics['vertical_resolution'],
        metrics['temperature_smoothness'],
        metrics['completeness'],
        metrics['physical_consistency']
    ])

    return metrics

Implementation Details

Caching Strategy

To minimize API calls and improve performance, atmospheric data is cached intelligently:

class AtmosphericCache:
    def __init__(self, cache_dir='data/atmospheric_cache'):
        self.cache_dir = cache_dir
        self.cache_duration = {
            'current': 600,      # 10 minutes
            'forecast': 3600,    # 1 hour
            'gfs': 21600,       # 6 hours
        }

    def get_cached_data(self, source, lat, lon, time):
        """Retrieve cached data if valid"""
        cache_key = self._generate_key(source, lat, lon, time)
        cache_file = os.path.join(self.cache_dir, f"{cache_key}.pkl")

        if os.path.exists(cache_file):
            age = time.time() - os.path.getmtime(cache_file)
            if age < self.cache_duration.get(source, 3600):
                with open(cache_file, 'rb') as f:
                    return pickle.load(f)

        return None

Error Handling and Fallbacks

Robust error handling ensures the simulation continues even when data sources fail:

def get_atmospheric_profile_with_fallback(lat, lon, alt, time):
    """Get atmospheric profile with automatic fallback"""

    errors = []

    # Try each source in priority order
    for source in ['openweather', 'noaa', 'gfs', 'standard']:
        try:
            if source == 'standard':
                # Final fallback - always works
                return standard_atmosphere.get_profile(lat, lon, alt)

            profile = fetch_from_source(source, lat, lon, time)

            # Validate quality
            if profile.quality_score > 0.5:
                return profile
            else:
                errors.append(f"{source}: Low quality score {profile.quality_score}")

        except Exception as e:
            errors.append(f"{source}: {str(e)}")
            continue

    # Log all errors for debugging
    logger.warning(f"All atmospheric sources failed: {errors}")

    # Return standard atmosphere as last resort
    return standard_atmosphere.get_profile(lat, lon, alt)

Performance Optimization

Several techniques optimize atmospheric data processing:

async def fetch_all_sources_parallel(lat, lon, time):
    tasks = [
        fetch_openweather_async(lat, lon, time),
        fetch_noaa_async(lat, lon, time),
        fetch_gfs_async(lat, lon, time)
    ]
    results = await asyncio.gather(*tasks, return_exceptions=True)
    return process_results(results)