hfpytrace.utils

Package

Common utility functions for config, routes, interpolation, and data conversion.

Key Methods

Function read_params_2D()
Function resolve_config_path()
Function get_installed_config_path()
Function get_default_config_name()
Function load_config()
Function load_config_1D()
Function load_config_2D()
Function load_config_3D()
Function to_namespace()
Function build_route_from_cfg()
Function build_heights_from_cfg()
Function build_elevations_from_cfg()
Function build_freqs_from_cfg()

API

hfpytrace.utils

Utility functions and physical constants for hfpytrace.

Constants

dict

Physical constants used throughout the package:

============= ========================================= ================ Key Description Value / Units ============= ========================================= ================ kconst Conversion constant (plasma frequency) 80.6 (unitless) boltz Boltzmann constant 1.38066×10⁻²³ J/K h Planck constant 6.626×10⁻³⁴ J·s c Speed of light 2.9979×10⁸ m/s avo Avogadro's number 6.023×10²³ mol⁻¹ Re Earth radius 6.371×10⁶ m amu Atomic mass unit 1.6605×10⁻²⁷ kg q_e Electron charge 1.602×10⁻¹⁹ C m_e Electron mass 9.109×10⁻³¹ kg g Surface gravitational acceleration 9.81 m/s² eps0 Permittivity of free space ~8.854×10⁻¹² F/m R Ideal gas constant 8.31 J/(mol·K) ============= ========================================= ================

Key Functions

interpolate_by_altitude Altitude-interpolate (or extrapolate) a 1D vertical profile. extrap1d Linear-extrapolation wrapper around scipy.interpolate.interp1d. load_config / load_config_1D / load_config_2D / load_config_3D Load a JSON config file and return it as a SimpleNamespace. to_namespace Recursively convert dict/list to SimpleNamespace/list. resolve_config_path Resolve a config path from a user override or the installed package default.

get_gridded_parameters(q, xparam='beam', yparam='slist', zparam='v', r=0, rounding=True)

Method converts scans to "beam" and "slist" or gate

Source code in hfpytrace/utils.py

def get_gridded_parameters(
    q, xparam="beam", yparam="slist", zparam="v", r=0, rounding=True
):
    """
    Method converts scans to "beam" and "slist" or gate
    """
    plotParamDF = q[[xparam, yparam, zparam]]
    if rounding:
        plotParamDF.loc[:, xparam] = np.round(plotParamDF[xparam].tolist(), r)
        plotParamDF.loc[:, yparam] = np.round(plotParamDF[yparam].tolist(), r)
    plotParamDF = plotParamDF.groupby([xparam, yparam]).mean().reset_index()
    plotParamDF = plotParamDF[[xparam, yparam, zparam]].pivot(
        index=xparam, columns=yparam
    )
    x = plotParamDF.index.values
    y = plotParamDF.columns.levels[1].values
    X, Y = np.meshgrid(x, y)
    # Mask the nan values! pcolormesh can't handle them well!
    Z = np.ma.masked_where(
        np.isnan(plotParamDF[zparam].values), plotParamDF[zparam].values
    )
    return X, Y, Z

get_folder(rad, beam, date, model, base)

Get folder by date

Source code in hfpytrace/utils.py

def get_folder(rad, beam, date, model, base):
    """
    Get folder by date
    """
    import os

    fold = os.path.join(base, date.strftime("%Y-%m-%d"), rad, "%02d" % beam, model)
    os.makedirs(fold, exist_ok=True)
    return fold

get_hamsci_folder(source, date, model, base, call_sign=None)

Get folder by date

Source code in hfpytrace/utils.py

def get_hamsci_folder(source, date, model, base, call_sign=None):
    """
    Get folder by date
    """
    import os

    fold = os.path.join(base, date.strftime("%Y-%m-%d"), source, model)
    fold = os.path.join(fold, call_sign) if call_sign else fold
    os.makedirs(fold, exist_ok=True)
    return fold

get_installed_config_path(config_name)

Return the installed package path for a config file under hfpytrace/cfg.

Source code in hfpytrace/utils.py

def get_installed_config_path(config_name: str) -> Path:
    """
    Return the installed package path for a config file under hfpytrace/cfg.
    """
    try:
        cfg = resources.files("hfpytrace").joinpath("cfg", config_name)
        return Path(str(cfg))
    except Exception:
        return Path(__file__).resolve().parent / "cfg" / config_name

resolve_config_path(config_path, default_config_name)

Resolve config path from user-provided value or package default.

Rules: - If config_path is provided and exists, use it. - If config_path is provided but missing, try hfpytrace/cfg/. - If config_path is None, use hfpytrace/cfg/.

Source code in hfpytrace/utils.py

def resolve_config_path(
    config_path: str | Path | None,
    default_config_name: str,
) -> Path:
    """
    Resolve config path from user-provided value or package default.

    Rules:
    - If `config_path` is provided and exists, use it.
    - If `config_path` is provided but missing, try hfpytrace/cfg/<basename>.
    - If `config_path` is None, use hfpytrace/cfg/<default_config_name>.
    """
    if config_path is not None:
        p = Path(config_path).expanduser()
        if p.exists():
            return p.resolve()
        pkg_try = get_installed_config_path(p.name)
        if pkg_try.exists():
            return pkg_try.resolve()
        raise FileNotFoundError(
            f"Config not found: {p}. Also tried package default: {pkg_try}"
        )

    default_path = get_installed_config_path(default_config_name)
    if not default_path.exists():
        raise FileNotFoundError(f"Package default config not found: {default_path}")
    return default_path.resolve()

get_default_config_name(config_dim)

Return bundled config filename for dimension selector.

Accepted values: - 1, "1d" - 2, "2d" - 3, "3d"

Source code in hfpytrace/utils.py

def get_default_config_name(config_dim: str | int) -> str:
    """
    Return bundled config filename for dimension selector.

    Accepted values:
    - 1, "1d"
    - 2, "2d"
    - 3, "3d"
    """
    key = str(config_dim).strip().lower()
    names = {
        "1": "config1D.json",
        "1d": "config1D.json",
        "2": "config2D.json",
        "2d": "config2D.json",
        "3": "config3D.json",
        "3d": "config3D.json",
    }
    if key not in names:
        raise ValueError("config_dim must be one of: 1/1d, 2/2d, 3/3d")
    return names[key]

load_config(config_dim, config_path=None)

Load hfpytrace config namespace for requested dimension.

Source code in hfpytrace/utils.py

def load_config(
    config_dim: str | int,
    config_path: str | Path | None = None,
):
    """
    Load hfpytrace config namespace for requested dimension.
    """
    default_name = get_default_config_name(config_dim)
    resolved = resolve_config_path(config_path, default_name)
    return read_params(resolved)

to_namespace(obj)

Recursively convert dict/list containers to SimpleNamespace/list.

Source code in hfpytrace/utils.py

def to_namespace(obj):
    """
    Recursively convert dict/list containers to SimpleNamespace/list.
    """
    if isinstance(obj, dict):
        return SimpleNamespace(**{k: to_namespace(v) for k, v in obj.items()})
    if isinstance(obj, list):
        return [to_namespace(v) for v in obj]
    return obj

extrap1d(x, y, kind='linear')

This method is used to extrapolate 1D paramteres

Source code in hfpytrace/utils.py

def extrap1d(x, y, kind="linear"):
    """This method is used to extrapolate 1D paramteres"""
    interpolator = interp1d(x, y, kind=kind)
    xs = interpolator.x
    ys = interpolator.y

    def pointwise(x):
        if x < xs[0]:
            return ys[0] + (x - xs[0]) * (ys[1] - ys[0]) / (xs[1] - xs[0])
        elif x > xs[-1]:
            return ys[-1] + (x - xs[-1]) * (ys[-1] - ys[-2]) / (xs[-1] - xs[-2])
        else:
            return interpolator(x)

    def ufunclike(xs):
        return np.array(list(map(pointwise, np.array(xs))))

    return ufunclike

interpolate_by_altitude(h, hx, param, scale='log', kind='cubic', method='intp')

Altitude-interpolate (or extrapolate) a 1D vertical profile.

Parameters

array-like

Source altitude axis [km] — must be strictly increasing.

array-like

Target altitude levels [km].

array-like

Profile values at h (must be > 0 when scale='log').

{'log', 'linear'}

Interpolation scale. 'log' operates on log10(param) so exponentially varying profiles (e.g. Ne) are well-behaved.

str

Interpolation kind passed to scipy.interpolate.interp1d (e.g. 'linear', 'cubic'). Default 'cubic'.

{'intp', 'extp'}

'intp' uses interp1d (raises if hx is outside h); 'extp' uses :func:extrap1d which linearly extends past the source bounds.

Returns

np.ndarray Interpolated (or extrapolated) profile at hx.

Source code in hfpytrace/utils.py

def interpolate_by_altitude(h, hx, param, scale="log", kind="cubic", method="intp"):
    """Altitude-interpolate (or extrapolate) a 1D vertical profile.

    Parameters
    ----------
    h : array-like
        Source altitude axis [km] — must be strictly increasing.
    hx : array-like
        Target altitude levels [km].
    param : array-like
        Profile values at ``h`` (must be > 0 when ``scale='log'``).
    scale : {'log', 'linear'}
        Interpolation scale.  ``'log'`` operates on ``log10(param)`` so
        exponentially varying profiles (e.g. Ne) are well-behaved.
    kind : str
        Interpolation kind passed to ``scipy.interpolate.interp1d``
        (e.g. ``'linear'``, ``'cubic'``).  Default ``'cubic'``.
    method : {'intp', 'extp'}
        ``'intp'`` uses ``interp1d`` (raises if ``hx`` is outside ``h``);
        ``'extp'`` uses :func:`extrap1d` which linearly extends past the
        source bounds.

    Returns
    -------
    np.ndarray
        Interpolated (or extrapolated) profile at ``hx``.
    """
    if scale == "linear":
        pnew = (
            interp1d(h, param, kind=kind)(hx)
            if method == "intp"
            else extrap1d(h, param, kind=kind)(hx)
        )
    if scale == "log":
        pnew = (
            10 ** interp1d(h, np.log10(param), kind=kind)(hx)
            if method == "intp"
            else 10 ** extrap1d(h, np.log10(param), kind=kind)(hx)
        )
    return pnew

create_movie(folder, outfile, pat, fps=3, out_fold=None)

Create movies from pngs

Source code in hfpytrace/utils.py

def create_movie(folder, outfile, pat, fps=3, out_fold=None):
    """
    Create movies from pngs
    """
    files = glob.glob(f"{folder}/{pat}")
    files.sort()

    out_fold = out_fold if out_fold else folder
    fourcc = cv2.VideoWriter_fourcc(*"XVID")
    img = cv2.imread(files[0])
    height, width, layers = img.shape
    size = (int(width / 2) * 2, int(height / 2) * 2)
    video = cv2.VideoWriter(out_fold + "/" + outfile, fourcc, fps, size)
    for idx in range(len(files)):
        im = cv2.imread(files[idx])
        im = cv2.resize(im, size)
        video.write(im)
        del im
    video.release()
    return

loadmatlabfiles(filename)

Load a .mat file and convert mat-objects to nested dictionaries.

Source code in hfpytrace/utils.py

def loadmatlabfiles(filename):
    """
    Load a .mat file and convert mat-objects to nested dictionaries.
    """
    data = spio.loadmat(filename, struct_as_record=False, squeeze_me=True)
    return _check_key_(data)

great_circle_distance(lat1, lon1, lat2, lon2, R=6371)

Calculate the great-circle distance between two points on Earth.

Parameters:

Name	Type	Description	Default
lat1	float	Latitude of the first point in degrees.	required
lon1	float	Longitude of the first point in degrees.	required
lat2	float	Latitude of the second point in degrees.	required
lon2	float	Longitude of the second point in degrees.	required

Returns:

Name	Type	Description
float		The great-circle distance in kilometers.

Source code in hfpytrace/utils.py

def great_circle_distance(
    lat1: float,
    lon1: float,
    lat2: float,
    lon2: float,
    R=6371,  # Radius of the Earth in kilometers
):
    """
    Calculate the great-circle distance between two points on Earth.

    Args:
        lat1 (float): Latitude of the first point in degrees.
        lon1 (float): Longitude of the first point in degrees.
        lat2 (float): Latitude of the second point in degrees.
        lon2 (float): Longitude of the second point in degrees.

    Returns:
        float: The great-circle distance in kilometers.
    """

    # Convert degrees to radians
    lat1 = math.radians(lat1)
    lon1 = math.radians(lon1)
    lat2 = math.radians(lat2)
    lon2 = math.radians(lon2)

    # Haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = (
        math.sin(dlat / 2) ** 2
        + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    )
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    distance = R * c

    return distance

build_route_from_cfg(cfg, n_range)

Build route samples along a great-circle path using config route fields.

Priority: 1) If cfg.route.end exists, use start -> end and compute initial bearing. 2) Otherwise use cfg.route.bearing from cfg.route.start.

Returns:

Type	Description
	lats (np.ndarray), lons (np.ndarray), bearing_deg (float), total_km (float)

Source code in hfpytrace/utils.py

def build_route_from_cfg(cfg, n_range: int):
    """
    Build route samples along a great-circle path using config route fields.

    Priority:
    1) If cfg.route.end exists, use start -> end and compute initial bearing.
    2) Otherwise use cfg.route.bearing from cfg.route.start.

    Returns:
        lats (np.ndarray), lons (np.ndarray), bearing_deg (float), total_km (float)
    """
    from geopy.distance import great_circle

    if n_range < 2:
        raise ValueError("n_range must be >= 2")

    start = cfg.route.start
    lat1 = float(start.lat)
    lon1 = float(start.lon)

    has_end = (
        hasattr(cfg.route, "end")
        and cfg.route.end is not None
        and hasattr(cfg.route.end, "lat")
        and hasattr(cfg.route.end, "lon")
    )

    if has_end:
        end = cfg.route.end
        lat2 = float(end.lat)
        lon2 = float(end.lon)
        bearing_deg, _ = calculate_bearing(lat1, lon1, lat2, lon2)
        total_km = great_circle_distance(lat1, lon1, lat2, lon2)
    else:
        if not hasattr(cfg.route, "bearing"):
            raise ValueError("cfg.route must include `bearing` when `end` is absent")
        if not hasattr(cfg, "max_ground_range_km"):
            raise ValueError(
                "cfg must include `max_ground_range_km` when using route bearing"
            )
        bearing_deg = float(cfg.route.bearing)
        total_km = float(cfg.max_ground_range_km)

    dists_km = np.linspace(0.0, total_km, n_range)
    lats = np.empty(n_range, dtype=float)
    lons = np.empty(n_range, dtype=float)
    origin = (lat1, lon1)
    for i, d_km in enumerate(dists_km):
        p = great_circle(kilometers=float(d_km)).destination(
            origin, bearing=bearing_deg
        )
        lats[i] = p.latitude
        lons[i] = p.longitude

    return lats, lons, float(bearing_deg), float(total_km)

build_heights_from_cfg(cfg)

Build height grid [km] from cfg start/end/increment fields.

Source code in hfpytrace/utils.py

def build_heights_from_cfg(cfg) -> np.ndarray:
    """Build height grid [km] from cfg start/end/increment fields."""
    h0 = float(cfg.start_height_km)
    h1 = float(cfg.end_height_km)
    dh = float(getattr(cfg, "height_incriment_km", 1.0))
    return np.arange(h0, h1, dh, dtype=float)

build_elevations_from_cfg(cfg)

Build elevation fan [deg] from cfg start/end/increment fields.

Source code in hfpytrace/utils.py

def build_elevations_from_cfg(cfg) -> np.ndarray:
    """Build elevation fan [deg] from cfg start/end/increment fields."""
    e0 = float(cfg.start_elevation)
    e1 = float(cfg.end_elevation)
    de = float(getattr(cfg, "elevation_inctiment", 0.1))
    return np.arange(e0, e1 + (0.5 * de), de, dtype=float)

build_freqs_from_cfg(cfg, elevs)

Build frequency vector [MHz] with one entry per elevation.

Source code in hfpytrace/utils.py

def build_freqs_from_cfg(cfg, elevs: np.ndarray) -> np.ndarray:
    """Build frequency vector [MHz] with one entry per elevation."""
    return np.full_like(elevs, float(cfg.frequency), dtype=float)

Source Code

#!/usr/bin/env python

"""Utility functions and physical constants for hfpytrace.

Constants
---------
pconst : dict
    Physical constants used throughout the package:

    ============= ========================================= ================
    Key           Description                               Value / Units
    ============= ========================================= ================
    ``kconst``    Conversion constant (plasma frequency)    80.6 (unitless)
    ``boltz``     Boltzmann constant                        1.38066×10⁻²³ J/K
    ``h``         Planck constant                           6.626×10⁻³⁴ J·s
    ``c``         Speed of light                            2.9979×10⁸ m/s
    ``avo``       Avogadro's number                         6.023×10²³ mol⁻¹
    ``Re``        Earth radius                              6.371×10⁶ m
    ``amu``       Atomic mass unit                          1.6605×10⁻²⁷ kg
    ``q_e``       Electron charge                           1.602×10⁻¹⁹ C
    ``m_e``       Electron mass                             9.109×10⁻³¹ kg
    ``g``         Surface gravitational acceleration        9.81 m/s²
    ``eps0``      Permittivity of free space                ~8.854×10⁻¹² F/m
    ``R``         Ideal gas constant                        8.31 J/(mol·K)
    ============= ========================================= ================

Key Functions
-------------
interpolate_by_altitude
    Altitude-interpolate (or extrapolate) a 1D vertical profile.
extrap1d
    Linear-extrapolation wrapper around ``scipy.interpolate.interp1d``.
load_config / load_config_1D / load_config_2D / load_config_3D
    Load a JSON config file and return it as a ``SimpleNamespace``.
to_namespace
    Recursively convert ``dict``/``list`` to ``SimpleNamespace``/``list``.
resolve_config_path
    Resolve a config path from a user override or the installed package default.
"""

__author__ = "Chakraborty, S."
__copyright__ = ""
__credits__ = []
__license__ = "MIT"
__version__ = "1.0."
__maintainer__ = "Chakraborty, S."
__email__ = "chakras4@erau.edu"
__status__ = "Research"

import glob
import json
import math
import os
from importlib import resources
from pathlib import Path
from types import SimpleNamespace

import numpy as np
import scipy.io as spio
from loguru import logger
from scipy.interpolate import interp1d

pconst = {
    "kconst": 80.6,  # Converstion constant
    "boltz": 1.38066e-23,  # Boltzmann constant  in Jule K^-1
    "h": 6.626e-34,  # Planks constant  in ergs s
    "c": 2.9979e08,  # in m s^-1
    "avo": 6.023e23,  # avogadro's number
    "Re": 6371.0e3,
    "amu": 1.6605e-27,
    "q_e": 1.602e-19,  # Electron charge in C
    "m_e": 9.109e-31,  # Electron mass in kg
    "g": 9.81,  # Gravitational acceleration on the surface of the Earth
    "eps0": 1e-9 / (36 * np.pi),
    "R": 8.31,  # J mol^-1 K^-1
}


def get_gridded_parameters(
    q, xparam="beam", yparam="slist", zparam="v", r=0, rounding=True
):
    """
    Method converts scans to "beam" and "slist" or gate
    """
    plotParamDF = q[[xparam, yparam, zparam]]
    if rounding:
        plotParamDF.loc[:, xparam] = np.round(plotParamDF[xparam].tolist(), r)
        plotParamDF.loc[:, yparam] = np.round(plotParamDF[yparam].tolist(), r)
    plotParamDF = plotParamDF.groupby([xparam, yparam]).mean().reset_index()
    plotParamDF = plotParamDF[[xparam, yparam, zparam]].pivot(
        index=xparam, columns=yparam
    )
    x = plotParamDF.index.values
    y = plotParamDF.columns.levels[1].values
    X, Y = np.meshgrid(x, y)
    # Mask the nan values! pcolormesh can't handle them well!
    Z = np.ma.masked_where(
        np.isnan(plotParamDF[zparam].values), plotParamDF[zparam].values
    )
    return X, Y, Z


def get_folder(rad, beam, date, model, base):
    """
    Get folder by date
    """
    import os

    fold = os.path.join(base, date.strftime("%Y-%m-%d"), rad, "%02d" % beam, model)
    os.makedirs(fold, exist_ok=True)
    return fold


def get_hamsci_folder(source, date, model, base, call_sign=None):
    """
    Get folder by date
    """
    import os

    fold = os.path.join(base, date.strftime("%Y-%m-%d"), source, model)
    fold = os.path.join(fold, call_sign) if call_sign else fold
    os.makedirs(fold, exist_ok=True)
    return fold


def read_params(
    fname: str | Path | None = None,
    default_config_name: str | None = "config2D.json",
):
    if fname is None:
        if default_config_name is None:
            raise ValueError("Provide `fname` or `default_config_name`.")
        fname = resolve_config_path(None, default_config_name)
    logger.info(f"Load config files: {fname}")
    with open(fname, "r") as f:
        param = json.load(f, object_hook=lambda x: SimpleNamespace(**x))
    return param


def get_installed_config_path(config_name: str) -> Path:
    """
    Return the installed package path for a config file under hfpytrace/cfg.
    """
    try:
        cfg = resources.files("hfpytrace").joinpath("cfg", config_name)
        return Path(str(cfg))
    except Exception:
        return Path(__file__).resolve().parent / "cfg" / config_name


def resolve_config_path(
    config_path: str | Path | None,
    default_config_name: str,
) -> Path:
    """
    Resolve config path from user-provided value or package default.

    Rules:
    - If `config_path` is provided and exists, use it.
    - If `config_path` is provided but missing, try hfpytrace/cfg/<basename>.
    - If `config_path` is None, use hfpytrace/cfg/<default_config_name>.
    """
    if config_path is not None:
        p = Path(config_path).expanduser()
        if p.exists():
            return p.resolve()
        pkg_try = get_installed_config_path(p.name)
        if pkg_try.exists():
            return pkg_try.resolve()
        raise FileNotFoundError(
            f"Config not found: {p}. Also tried package default: {pkg_try}"
        )

    default_path = get_installed_config_path(default_config_name)
    if not default_path.exists():
        raise FileNotFoundError(f"Package default config not found: {default_path}")
    return default_path.resolve()


def get_default_config_name(config_dim: str | int) -> str:
    """
    Return bundled config filename for dimension selector.

    Accepted values:
    - 1, "1d"
    - 2, "2d"
    - 3, "3d"
    """
    key = str(config_dim).strip().lower()
    names = {
        "1": "config1D.json",
        "1d": "config1D.json",
        "2": "config2D.json",
        "2d": "config2D.json",
        "3": "config3D.json",
        "3d": "config3D.json",
    }
    if key not in names:
        raise ValueError("config_dim must be one of: 1/1d, 2/2d, 3/3d")
    return names[key]


def load_config(
    config_dim: str | int,
    config_path: str | Path | None = None,
):
    """
    Load hfpytrace config namespace for requested dimension.
    """
    default_name = get_default_config_name(config_dim)
    resolved = resolve_config_path(config_path, default_name)
    return read_params(resolved)


def load_config_1D(config_path: str | Path | None = None):
    return load_config("1d", config_path=config_path)


def load_config_2D(config_path: str | Path | None = None):
    return load_config("2d", config_path=config_path)


def load_config_3D(config_path: str | Path | None = None):
    return load_config("3d", config_path=config_path)


def to_namespace(obj):
    """
    Recursively convert dict/list containers to SimpleNamespace/list.
    """
    if isinstance(obj, dict):
        return SimpleNamespace(**{k: to_namespace(v) for k, v in obj.items()})
    if isinstance(obj, list):
        return [to_namespace(v) for v in obj]
    return obj


def clean():
    files = glob.glob(str(Path.home() / "matlab_crash_dump*"))
    for f in files:
        if os.path.exists(f):
            os.remove(f)
    return


def extrap1d(x, y, kind="linear"):
    """This method is used to extrapolate 1D paramteres"""
    interpolator = interp1d(x, y, kind=kind)
    xs = interpolator.x
    ys = interpolator.y

    def pointwise(x):
        if x < xs[0]:
            return ys[0] + (x - xs[0]) * (ys[1] - ys[0]) / (xs[1] - xs[0])
        elif x > xs[-1]:
            return ys[-1] + (x - xs[-1]) * (ys[-1] - ys[-2]) / (xs[-1] - xs[-2])
        else:
            return interpolator(x)

    def ufunclike(xs):
        return np.array(list(map(pointwise, np.array(xs))))

    return ufunclike


def interpolate_by_altitude(h, hx, param, scale="log", kind="cubic", method="intp"):
    """Altitude-interpolate (or extrapolate) a 1D vertical profile.

    Parameters
    ----------
    h : array-like
        Source altitude axis [km] — must be strictly increasing.
    hx : array-like
        Target altitude levels [km].
    param : array-like
        Profile values at ``h`` (must be > 0 when ``scale='log'``).
    scale : {'log', 'linear'}
        Interpolation scale.  ``'log'`` operates on ``log10(param)`` so
        exponentially varying profiles (e.g. Ne) are well-behaved.
    kind : str
        Interpolation kind passed to ``scipy.interpolate.interp1d``
        (e.g. ``'linear'``, ``'cubic'``).  Default ``'cubic'``.
    method : {'intp', 'extp'}
        ``'intp'`` uses ``interp1d`` (raises if ``hx`` is outside ``h``);
        ``'extp'`` uses :func:`extrap1d` which linearly extends past the
        source bounds.

    Returns
    -------
    np.ndarray
        Interpolated (or extrapolated) profile at ``hx``.
    """
    if scale == "linear":
        pnew = (
            interp1d(h, param, kind=kind)(hx)
            if method == "intp"
            else extrap1d(h, param, kind=kind)(hx)
        )
    if scale == "log":
        pnew = (
            10 ** interp1d(h, np.log10(param), kind=kind)(hx)
            if method == "intp"
            else 10 ** extrap1d(h, np.log10(param), kind=kind)(hx)
        )
    return pnew


def smooth(x, window_len=11, window="hanning"):
    if x.ndim != 1:
        raise ValueError("smooth only accepts 1 dimension arrays.")
    if x.size < window_len:
        raise ValueError("Input vector needs to be bigger than window size.")
    if window_len < 3:
        return x
    if not window in ["flat", "hanning", "hamming", "bartlett", "blackman"]:
        raise ValueError(
            "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
        )
    s = np.r_[x[window_len - 1 : 0 : -1], x, x[-2 : -window_len - 1 : -1]]
    if window == "flat":
        w = np.ones(window_len, "d")
    else:
        w = eval("np." + window + "(window_len)")
    y = np.convolve(w / w.sum(), s, mode="valid")
    d = window_len - 1
    y = y[int(d / 2) : -int(d / 2)]
    return y


def create_movie(folder, outfile, pat, fps=3, out_fold=None):
    """
    Create movies from pngs
    """
    files = glob.glob(f"{folder}/{pat}")
    files.sort()

    out_fold = out_fold if out_fold else folder
    fourcc = cv2.VideoWriter_fourcc(*"XVID")
    img = cv2.imread(files[0])
    height, width, layers = img.shape
    size = (int(width / 2) * 2, int(height / 2) * 2)
    video = cv2.VideoWriter(out_fold + "/" + outfile, fourcc, fps, size)
    for idx in range(len(files)):
        im = cv2.imread(files[idx])
        im = cv2.resize(im, size)
        video.write(im)
        del im
    video.release()
    return


def create_mp4(folder, pat, outfile, fps=3, out_fold=None):
    import os

    out_fold = out_fold if out_fold else folder
    cmd = f"""
    ffmpeg -framerate {fps} -pattern_type glob -i '{folder}/{pat}' -c:v libx264 {out_fold}/{outfile}
    """
    os.system(cmd)
    return


def loadmatlabfiles(filename):
    """
    Load a .mat file and convert mat-objects to nested dictionaries.
    """
    data = spio.loadmat(filename, struct_as_record=False, squeeze_me=True)
    return _check_key_(data)


def _check_key_(data_dict):
    """
    Recursively check and convert mat-objects in a dictionary to nested dictionaries.
    """
    for key, value in data_dict.items():
        if isinstance(value, spio.matlab.mio5_params.mat_struct):
            data_dict[key] = _todict_(value)
    return data_dict


def _todict_(matobj):
    """
    Recursively construct nested dictionaries from mat-objects.
    """
    return {
        strg: (
            _todict_(getattr(matobj, strg))
            if isinstance(getattr(matobj, strg), spio.matlab.mio5_params.mat_struct)
            else getattr(matobj, strg)
        )
        for strg in matobj._fieldnames
    }


def setsize(size: float = 8):
    import matplotlib as mpl

    mpl.rcParams.update(
        {"xtick.labelsize": size, "ytick.labelsize": size, "font.size": size}
    )
    return


def calculate_bearing(lat1: float, lon1: float, lat2: float, lon2: float):
    lat1, lon1, lat2, lon2 = map(math.radians, [lat1, lon1, lat2, lon2])

    dlon = lon2 - lon1
    y = math.sin(dlon) * math.cos(lat2)
    x = math.cos(lat1) * math.sin(lat2) - math.sin(lat1) * math.cos(lat2) * math.cos(
        dlon
    )
    bearing = math.atan2(y, x)

    return (math.degrees(bearing) + 360) % 360, math.degrees(bearing)


def great_circle_distance(
    lat1: float,
    lon1: float,
    lat2: float,
    lon2: float,
    R=6371,  # Radius of the Earth in kilometers
):
    """
    Calculate the great-circle distance between two points on Earth.

    Args:
        lat1 (float): Latitude of the first point in degrees.
        lon1 (float): Longitude of the first point in degrees.
        lat2 (float): Latitude of the second point in degrees.
        lon2 (float): Longitude of the second point in degrees.

    Returns:
        float: The great-circle distance in kilometers.
    """

    # Convert degrees to radians
    lat1 = math.radians(lat1)
    lon1 = math.radians(lon1)
    lat2 = math.radians(lat2)
    lon2 = math.radians(lon2)

    # Haversine formula
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = (
        math.sin(dlat / 2) ** 2
        + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2
    )
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    distance = R * c

    return distance


def build_route_from_cfg(cfg, n_range: int):
    """
    Build route samples along a great-circle path using config route fields.

    Priority:
    1) If cfg.route.end exists, use start -> end and compute initial bearing.
    2) Otherwise use cfg.route.bearing from cfg.route.start.

    Returns:
        lats (np.ndarray), lons (np.ndarray), bearing_deg (float), total_km (float)
    """
    from geopy.distance import great_circle

    if n_range < 2:
        raise ValueError("n_range must be >= 2")

    start = cfg.route.start
    lat1 = float(start.lat)
    lon1 = float(start.lon)

    has_end = (
        hasattr(cfg.route, "end")
        and cfg.route.end is not None
        and hasattr(cfg.route.end, "lat")
        and hasattr(cfg.route.end, "lon")
    )

    if has_end:
        end = cfg.route.end
        lat2 = float(end.lat)
        lon2 = float(end.lon)
        bearing_deg, _ = calculate_bearing(lat1, lon1, lat2, lon2)
        total_km = great_circle_distance(lat1, lon1, lat2, lon2)
    else:
        if not hasattr(cfg.route, "bearing"):
            raise ValueError("cfg.route must include `bearing` when `end` is absent")
        if not hasattr(cfg, "max_ground_range_km"):
            raise ValueError(
                "cfg must include `max_ground_range_km` when using route bearing"
            )
        bearing_deg = float(cfg.route.bearing)
        total_km = float(cfg.max_ground_range_km)

    dists_km = np.linspace(0.0, total_km, n_range)
    lats = np.empty(n_range, dtype=float)
    lons = np.empty(n_range, dtype=float)
    origin = (lat1, lon1)
    for i, d_km in enumerate(dists_km):
        p = great_circle(kilometers=float(d_km)).destination(
            origin, bearing=bearing_deg
        )
        lats[i] = p.latitude
        lons[i] = p.longitude

    return lats, lons, float(bearing_deg), float(total_km)


def build_heights_from_cfg(cfg) -> np.ndarray:
    """Build height grid [km] from cfg start/end/increment fields."""
    h0 = float(cfg.start_height_km)
    h1 = float(cfg.end_height_km)
    dh = float(getattr(cfg, "height_incriment_km", 1.0))
    return np.arange(h0, h1, dh, dtype=float)


def build_elevations_from_cfg(cfg) -> np.ndarray:
    """Build elevation fan [deg] from cfg start/end/increment fields."""
    e0 = float(cfg.start_elevation)
    e1 = float(cfg.end_elevation)
    de = float(getattr(cfg, "elevation_inctiment", 0.1))
    return np.arange(e0, e1 + (0.5 * de), de, dtype=float)


def build_freqs_from_cfg(cfg, elevs: np.ndarray) -> np.ndarray:
    """Build frequency vector [MHz] with one entry per elevation."""
    return np.full_like(elevs, float(cfg.frequency), dtype=float)