Skip to content

hfpytrace.density.gemini

Package

GEMINI model adapter for 2D and 3D electron density grid workflows.

Key Classes

Class GEMINI2d

Key Methods

Method GEMINI2d.fetch_dataset()
Method GEMINI2d.fetch_dataset_3d()

API

hfpytrace.density.gemini

GEMINI electron density reader.

Reads GEMINI (Geospace Environment Model of Ion-Neutral Interactions) HDF5 output files and interpolates the electron density onto arbitrary lat/lon/altitude grids for use in HF ray-tracing.

Requires

h5py, pandas : for HDF5 I/O and grid operations.

Classes

GEMINI2d Loads GEMINI HDF5 output files and provides fetch_dataset / fetch_dataset_3d methods that resample Ne onto user-specified grids.

GEMINI2d

Bases: object

Electron density from GEMINI simulation output.

Reads GEMINI HDF5 output and interpolates electron density to requested (lat, lon, alt) grids. The grid coordinate file (cfg.grid_coordinate_file) is loaded once; individual time-step files are matched by wall-clock time.

Parameters
SimpleNamespace

Config with density_file_location, density_file_name (glob pattern for data files), grid_coordinate_file (HDF5 coordinate file), scale, and kind (interpolation scale/kind).

datetime.datetime

Reference event time; the nearest file timestamp is used.

Source code in hfpytrace/density/gemini.py 
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
class GEMINI2d(object):
    """Electron density from GEMINI simulation output.

    Reads GEMINI HDF5 output and interpolates electron density to requested
    (lat, lon, alt) grids.  The grid coordinate file (``cfg.grid_coordinate_file``)
    is loaded once; individual time-step files are matched by wall-clock time.

    Parameters
    ----------
    cfg : SimpleNamespace
        Config with ``density_file_location``, ``density_file_name`` (glob
        pattern for data files), ``grid_coordinate_file`` (HDF5 coordinate
        file), ``scale``, and ``kind`` (interpolation scale/kind).
    event : datetime.datetime
        Reference event time; the nearest file timestamp is used.
    """

    def __init__(
        self,
        cfg,
        event,
    ):
        self.cfg = cfg
        self.event = event
        self.load_grid()
        self.search_mat_files()
        return

    def search_mat_files(self):
        """
        Search all files
        """
        self.files = glob.glob(os.path.join(self.cfg.density_file_location, "*.mat"))
        self.files.remove(self.ccord_file)
        self.files.sort()
        self.dates = []
        for fname in self.files:
            day = dt.datetime.strptime(fname.split("_")[0].split("/")[-1], "%Y%m%d")
            seconds = int(fname.split("_")[1].split(".")[0])
            day = day + dt.timedelta(seconds=seconds)
            self.dates.append(day)
        return

    def load_grid(self):
        """
        Load grid file
        """
        logger.info(f"Loading GEMINI Grid files")
        self.ccord_file = self.cfg.density_file_location + self.cfg.grid_coordinate_file
        self.grid = pd.DataFrame()
        with h5py.File(self.ccord_file, "r") as fkey:
            self.grid["glat"] = np.array(fkey.get("glat")[0]).tolist()
            self.grid["glon"] = np.mod((fkey.get("glon")[0] + 180), 360) - 180
            self.grid["alt"] = fkey.get("galt")[0, :]
            self.grid["glat"] = fkey.get("glat")[0, :]
        return

    def load_data(self, fname):
        """Load electron density array from a single GEMINI HDF5 time-step file.

        Parameters
        ----------
        fname : str
            Path to the HDF5 data file.

        Returns
        -------
        pandas.DataFrame
            Grid DataFrame extended with a ``nsall`` column (electron density in m⁻³).
        """
        with h5py.File(fname, "r") as fkey:
            o = self.grid.copy()
            o["nsall"] = fkey.get("nsall")[0, :]
            logger.info(f"Load dataset: {fname}")
        return o

    def fetch_dataset(
        self,
        time: dt.datetime,
        lats,
        lons,
        alts,
        to_file: str = None,
        dlat: float = 0.2,
        dlon: float = 0.2,
        intp_edens_xlim_index: float = 0,
    ):
        """Fetch 2D electron density along a route at a given time.

        The nearest GEMINI time-step file is selected.  For each route point
        ``(lat, lon)``, grid cells within ``±dlat`` / ``±dlon`` are averaged
        and then altitude-interpolated.

        Parameters
        ----------
        time : datetime.datetime
            Requested time snapshot.
        lats, lons : array-like, shape (npts,)
            Geographic coordinates of route points [°].
        alts : array-like, shape (nalt,)
            Target altitude levels [km].
        to_file : str, optional
            If given, save the result to a ``.mat`` file at this path.
        dlat, dlon : float
            Half-width of the spatial averaging box [°].  Default ``0.2``.
        intp_edens_xlim_index : int
            If > 0, all columns from this index onward are filled by repeating
            the column at this index (useful for extending the ionosphere past
            the GEMINI domain boundary).

        Returns
        -------
        param : np.ndarray, shape (nalt, npts)
            Electron density in cm⁻³.
        alts : array-like
            The ``alts`` argument (returned for caller convenience).
        """
        self.param = np.zeros((len(alts), len(lats)))
        i = np.argmin([np.abs((t - time).total_seconds()) for t in self.dates])
        file = self.files[i]
        df = self.load_data(file)
        df = df[df.alt >= 0]
        self.param = np.zeros((len(alts), len(lats)))
        for ix, lat, lon in zip(range(len(lats)), lats, lons):
            uf = df[
                (df.glat >= lat - dlat)
                & (df.glat <= lat + dlat)
                & (df.glon >= lon - dlon)
                & (df.glon <= lon + dlon)
            ]
            uf = uf.groupby("alt").mean().reset_index()
            if (len(uf) > 0) and (uf.alt.max() / 1e3 >= max(alts) * 0.7):
                method = "intp" if uf.alt.max() / 1e3 > max(alts) else "extp"
                self.param[:, ix] = (
                    utils.interpolate_by_altitude(
                        np.array(uf.alt) / 1e3,
                        alts,
                        np.array(uf["nsall"]),
                        self.cfg.scale,
                        self.cfg.kind,
                        method=method,
                    )
                    * 1e-6
                )
        logger.info(f"Index: {intp_edens_xlim_index}")
        if intp_edens_xlim_index:
            for j in range(intp_edens_xlim_index, len(lats)):
                self.param[:, j] = self.param[:, intp_edens_xlim_index]
        if to_file:
            savemat(to_file, dict(ne=self.param))
        return self.param, alts

    def _fetch_profile_from_df(self, df, lat, lon, alts, dlat=0.2, dlon=0.2):
        """Return a single altitude-interpolated Ne profile at (lat, lon).

        Parameters
        ----------
        df : pandas.DataFrame
            Loaded grid/data frame with columns ``glat``, ``glon``, ``alt``,
            and ``nsall``.
        lat, lon : float
            Target geographic coordinates [°].
        alts : array-like
            Target altitude levels [km].
        dlat, dlon : float
            Spatial averaging half-width [°].

        Returns
        -------
        np.ndarray, shape (nalt,)
            Electron density in cm⁻³.
        """
        uf = df[
            (df.glat >= lat - dlat)
            & (df.glat <= lat + dlat)
            & (df.glon >= lon - dlon)
            & (df.glon <= lon + dlon)
        ]
        uf = uf.groupby("alt").mean().reset_index()
        p = np.zeros(len(alts), dtype=float)
        if (len(uf) > 0) and (uf.alt.max() / 1e3 >= max(alts) * 0.7):
            method = "intp" if uf.alt.max() / 1e3 > max(alts) else "extp"
            p = (
                utils.interpolate_by_altitude(
                    np.array(uf.alt) / 1e3,
                    alts,
                    np.array(uf["nsall"]),
                    self.cfg.scale,
                    self.cfg.kind,
                    method=method,
                )
                * 1e-6
            )
        return np.asarray(p, dtype=float)

    def fetch_dataset_3d(
        self,
        time: dt.datetime,
        lats,
        lons,
        alts,
        to_file: str = None,
        workers: int = 1,
        dlat: float = 0.2,
        dlon: float = 0.2,
    ):
        """
        Fetch 3D electron density cube on (lat, lon, alt) with optional
        point-wise parallelism.
        """
        lats = np.asarray(lats, dtype=float)
        lons = np.asarray(lons, dtype=float)
        alts = np.asarray(alts, dtype=float)
        i = np.argmin([np.abs((t - time).total_seconds()) for t in self.dates])
        file = self.files[i]
        df = self.load_data(file)
        df = df[df.alt >= 0]

        out = np.zeros((lats.size, lons.size, alts.size), dtype=float)
        ij = [(ii, jj) for ii in range(lats.size) for jj in range(lons.size)]

        def _job(ii, jj):
            p = self._fetch_profile_from_df(
                df, lats[ii], lons[jj], alts, dlat=dlat, dlon=dlon
            )
            return ii, jj, p

        n_workers = max(1, int(workers))
        if n_workers == 1:
            for ii, jj in ij:
                _, _, p = _job(ii, jj)
                out[ii, jj, :] = p
        else:
            with ThreadPoolExecutor(max_workers=n_workers) as ex:
                for ii, jj, p in ex.map(lambda t: _job(*t), ij):
                    out[ii, jj, :] = p

        self.param3d = out
        if to_file:
            savemat(to_file, dict(ne=self.param3d))
        return self.param3d, alts

    def load_from_file(self, to_file: str):
        """Load a previously saved electron density array from a ``.mat`` file.

        Parameters
        ----------
        to_file : str
            Path to the ``.mat`` file containing key ``ne``.

        Returns
        -------
        np.ndarray
            Electron density array as stored.
        """
        logger.info(f"Load from file {to_file.split('/')[-1]}")
        self.param = loadmat(to_file)["ne"]
        return self.param

search_mat_files()

Search all files

Source code in hfpytrace/density/gemini.py 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def search_mat_files(self):
    """
    Search all files
    """
    self.files = glob.glob(os.path.join(self.cfg.density_file_location, "*.mat"))
    self.files.remove(self.ccord_file)
    self.files.sort()
    self.dates = []
    for fname in self.files:
        day = dt.datetime.strptime(fname.split("_")[0].split("/")[-1], "%Y%m%d")
        seconds = int(fname.split("_")[1].split(".")[0])
        day = day + dt.timedelta(seconds=seconds)
        self.dates.append(day)
    return

load_grid()

Load grid file

Source code in hfpytrace/density/gemini.py 
75
76
77
78
79
80
81
82
83
84
85
86
87
def load_grid(self):
    """
    Load grid file
    """
    logger.info(f"Loading GEMINI Grid files")
    self.ccord_file = self.cfg.density_file_location + self.cfg.grid_coordinate_file
    self.grid = pd.DataFrame()
    with h5py.File(self.ccord_file, "r") as fkey:
        self.grid["glat"] = np.array(fkey.get("glat")[0]).tolist()
        self.grid["glon"] = np.mod((fkey.get("glon")[0] + 180), 360) - 180
        self.grid["alt"] = fkey.get("galt")[0, :]
        self.grid["glat"] = fkey.get("glat")[0, :]
    return

load_data(fname)

Load electron density array from a single GEMINI HDF5 time-step file.

Parameters
str

Path to the HDF5 data file.

Returns

pandas.DataFrame Grid DataFrame extended with a nsall column (electron density in m⁻³).

Source code in hfpytrace/density/gemini.py 
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
def load_data(self, fname):
    """Load electron density array from a single GEMINI HDF5 time-step file.

    Parameters
    ----------
    fname : str
        Path to the HDF5 data file.

    Returns
    -------
    pandas.DataFrame
        Grid DataFrame extended with a ``nsall`` column (electron density in m⁻³).
    """
    with h5py.File(fname, "r") as fkey:
        o = self.grid.copy()
        o["nsall"] = fkey.get("nsall")[0, :]
        logger.info(f"Load dataset: {fname}")
    return o

fetch_dataset(time, lats, lons, alts, to_file=None, dlat=0.2, dlon=0.2, intp_edens_xlim_index=0)

Fetch 2D electron density along a route at a given time.

The nearest GEMINI time-step file is selected. For each route point (lat, lon), grid cells within ±dlat / ±dlon are averaged and then altitude-interpolated.

Parameters
datetime.datetime

Requested time snapshot.

lats, lons : array-like, shape (npts,) Geographic coordinates of route points [°].

array-like, shape (nalt,)

Target altitude levels [km].

str, optional

If given, save the result to a .mat file at this path.

dlat, dlon : float Half-width of the spatial averaging box [°]. Default 0.2.

int

If > 0, all columns from this index onward are filled by repeating the column at this index (useful for extending the ionosphere past the GEMINI domain boundary).

Returns
np.ndarray, shape (nalt, npts)

Electron density in cm⁻³.

array-like

The alts argument (returned for caller convenience).

Source code in hfpytrace/density/gemini.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
def fetch_dataset(
    self,
    time: dt.datetime,
    lats,
    lons,
    alts,
    to_file: str = None,
    dlat: float = 0.2,
    dlon: float = 0.2,
    intp_edens_xlim_index: float = 0,
):
    """Fetch 2D electron density along a route at a given time.

    The nearest GEMINI time-step file is selected.  For each route point
    ``(lat, lon)``, grid cells within ``±dlat`` / ``±dlon`` are averaged
    and then altitude-interpolated.

    Parameters
    ----------
    time : datetime.datetime
        Requested time snapshot.
    lats, lons : array-like, shape (npts,)
        Geographic coordinates of route points [°].
    alts : array-like, shape (nalt,)
        Target altitude levels [km].
    to_file : str, optional
        If given, save the result to a ``.mat`` file at this path.
    dlat, dlon : float
        Half-width of the spatial averaging box [°].  Default ``0.2``.
    intp_edens_xlim_index : int
        If > 0, all columns from this index onward are filled by repeating
        the column at this index (useful for extending the ionosphere past
        the GEMINI domain boundary).

    Returns
    -------
    param : np.ndarray, shape (nalt, npts)
        Electron density in cm⁻³.
    alts : array-like
        The ``alts`` argument (returned for caller convenience).
    """
    self.param = np.zeros((len(alts), len(lats)))
    i = np.argmin([np.abs((t - time).total_seconds()) for t in self.dates])
    file = self.files[i]
    df = self.load_data(file)
    df = df[df.alt >= 0]
    self.param = np.zeros((len(alts), len(lats)))
    for ix, lat, lon in zip(range(len(lats)), lats, lons):
        uf = df[
            (df.glat >= lat - dlat)
            & (df.glat <= lat + dlat)
            & (df.glon >= lon - dlon)
            & (df.glon <= lon + dlon)
        ]
        uf = uf.groupby("alt").mean().reset_index()
        if (len(uf) > 0) and (uf.alt.max() / 1e3 >= max(alts) * 0.7):
            method = "intp" if uf.alt.max() / 1e3 > max(alts) else "extp"
            self.param[:, ix] = (
                utils.interpolate_by_altitude(
                    np.array(uf.alt) / 1e3,
                    alts,
                    np.array(uf["nsall"]),
                    self.cfg.scale,
                    self.cfg.kind,
                    method=method,
                )
                * 1e-6
            )
    logger.info(f"Index: {intp_edens_xlim_index}")
    if intp_edens_xlim_index:
        for j in range(intp_edens_xlim_index, len(lats)):
            self.param[:, j] = self.param[:, intp_edens_xlim_index]
    if to_file:
        savemat(to_file, dict(ne=self.param))
    return self.param, alts

fetch_dataset_3d(time, lats, lons, alts, to_file=None, workers=1, dlat=0.2, dlon=0.2)

Fetch 3D electron density cube on (lat, lon, alt) with optional point-wise parallelism.

Source code in hfpytrace/density/gemini.py 
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def fetch_dataset_3d(
    self,
    time: dt.datetime,
    lats,
    lons,
    alts,
    to_file: str = None,
    workers: int = 1,
    dlat: float = 0.2,
    dlon: float = 0.2,
):
    """
    Fetch 3D electron density cube on (lat, lon, alt) with optional
    point-wise parallelism.
    """
    lats = np.asarray(lats, dtype=float)
    lons = np.asarray(lons, dtype=float)
    alts = np.asarray(alts, dtype=float)
    i = np.argmin([np.abs((t - time).total_seconds()) for t in self.dates])
    file = self.files[i]
    df = self.load_data(file)
    df = df[df.alt >= 0]

    out = np.zeros((lats.size, lons.size, alts.size), dtype=float)
    ij = [(ii, jj) for ii in range(lats.size) for jj in range(lons.size)]

    def _job(ii, jj):
        p = self._fetch_profile_from_df(
            df, lats[ii], lons[jj], alts, dlat=dlat, dlon=dlon
        )
        return ii, jj, p

    n_workers = max(1, int(workers))
    if n_workers == 1:
        for ii, jj in ij:
            _, _, p = _job(ii, jj)
            out[ii, jj, :] = p
    else:
        with ThreadPoolExecutor(max_workers=n_workers) as ex:
            for ii, jj, p in ex.map(lambda t: _job(*t), ij):
                out[ii, jj, :] = p

    self.param3d = out
    if to_file:
        savemat(to_file, dict(ne=self.param3d))
    return self.param3d, alts

load_from_file(to_file)

Load a previously saved electron density array from a .mat file.

Parameters
str

Path to the .mat file containing key ne.

Returns

np.ndarray Electron density array as stored.

Source code in hfpytrace/density/gemini.py 
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
def load_from_file(self, to_file: str):
    """Load a previously saved electron density array from a ``.mat`` file.

    Parameters
    ----------
    to_file : str
        Path to the ``.mat`` file containing key ``ne``.

    Returns
    -------
    np.ndarray
        Electron density array as stored.
    """
    logger.info(f"Load from file {to_file.split('/')[-1]}")
    self.param = loadmat(to_file)["ne"]
    return self.param

Source Code

hfpytrace/density/gemini.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
"""GEMINI electron density reader.

Reads GEMINI (Geospace Environment Model of Ion-Neutral Interactions) HDF5
output files and interpolates the electron density onto arbitrary lat/lon/altitude
grids for use in HF ray-tracing.

Requires
--------
h5py, pandas : for HDF5 I/O and grid operations.

Classes
-------
GEMINI2d
    Loads GEMINI HDF5 output files and provides ``fetch_dataset`` /
    ``fetch_dataset_3d`` methods that resample Ne onto user-specified grids.
"""

import datetime as dt
import glob
import os
from concurrent.futures import ThreadPoolExecutor

import h5py
import numpy as np
import pandas as pd
from loguru import logger
from scipy.io import loadmat, savemat

from hfpytrace import utils


class GEMINI2d(object):
    """Electron density from GEMINI simulation output.

    Reads GEMINI HDF5 output and interpolates electron density to requested
    (lat, lon, alt) grids.  The grid coordinate file (``cfg.grid_coordinate_file``)
    is loaded once; individual time-step files are matched by wall-clock time.

    Parameters
    ----------
    cfg : SimpleNamespace
        Config with ``density_file_location``, ``density_file_name`` (glob
        pattern for data files), ``grid_coordinate_file`` (HDF5 coordinate
        file), ``scale``, and ``kind`` (interpolation scale/kind).
    event : datetime.datetime
        Reference event time; the nearest file timestamp is used.
    """

    def __init__(
        self,
        cfg,
        event,
    ):
        self.cfg = cfg
        self.event = event
        self.load_grid()
        self.search_mat_files()
        return

    def search_mat_files(self):
        """
        Search all files
        """
        self.files = glob.glob(os.path.join(self.cfg.density_file_location, "*.mat"))
        self.files.remove(self.ccord_file)
        self.files.sort()
        self.dates = []
        for fname in self.files:
            day = dt.datetime.strptime(fname.split("_")[0].split("/")[-1], "%Y%m%d")
            seconds = int(fname.split("_")[1].split(".")[0])
            day = day + dt.timedelta(seconds=seconds)
            self.dates.append(day)
        return

    def load_grid(self):
        """
        Load grid file
        """
        logger.info(f"Loading GEMINI Grid files")
        self.ccord_file = self.cfg.density_file_location + self.cfg.grid_coordinate_file
        self.grid = pd.DataFrame()
        with h5py.File(self.ccord_file, "r") as fkey:
            self.grid["glat"] = np.array(fkey.get("glat")[0]).tolist()
            self.grid["glon"] = np.mod((fkey.get("glon")[0] + 180), 360) - 180
            self.grid["alt"] = fkey.get("galt")[0, :]
            self.grid["glat"] = fkey.get("glat")[0, :]
        return

    def load_data(self, fname):
        """Load electron density array from a single GEMINI HDF5 time-step file.

        Parameters
        ----------
        fname : str
            Path to the HDF5 data file.

        Returns
        -------
        pandas.DataFrame
            Grid DataFrame extended with a ``nsall`` column (electron density in m⁻³).
        """
        with h5py.File(fname, "r") as fkey:
            o = self.grid.copy()
            o["nsall"] = fkey.get("nsall")[0, :]
            logger.info(f"Load dataset: {fname}")
        return o

    def fetch_dataset(
        self,
        time: dt.datetime,
        lats,
        lons,
        alts,
        to_file: str = None,
        dlat: float = 0.2,
        dlon: float = 0.2,
        intp_edens_xlim_index: float = 0,
    ):
        """Fetch 2D electron density along a route at a given time.

        The nearest GEMINI time-step file is selected.  For each route point
        ``(lat, lon)``, grid cells within ``±dlat`` / ``±dlon`` are averaged
        and then altitude-interpolated.

        Parameters
        ----------
        time : datetime.datetime
            Requested time snapshot.
        lats, lons : array-like, shape (npts,)
            Geographic coordinates of route points [°].
        alts : array-like, shape (nalt,)
            Target altitude levels [km].
        to_file : str, optional
            If given, save the result to a ``.mat`` file at this path.
        dlat, dlon : float
            Half-width of the spatial averaging box [°].  Default ``0.2``.
        intp_edens_xlim_index : int
            If > 0, all columns from this index onward are filled by repeating
            the column at this index (useful for extending the ionosphere past
            the GEMINI domain boundary).

        Returns
        -------
        param : np.ndarray, shape (nalt, npts)
            Electron density in cm⁻³.
        alts : array-like
            The ``alts`` argument (returned for caller convenience).
        """
        self.param = np.zeros((len(alts), len(lats)))
        i = np.argmin([np.abs((t - time).total_seconds()) for t in self.dates])
        file = self.files[i]
        df = self.load_data(file)
        df = df[df.alt >= 0]
        self.param = np.zeros((len(alts), len(lats)))
        for ix, lat, lon in zip(range(len(lats)), lats, lons):
            uf = df[
                (df.glat >= lat - dlat)
                & (df.glat <= lat + dlat)
                & (df.glon >= lon - dlon)
                & (df.glon <= lon + dlon)
            ]
            uf = uf.groupby("alt").mean().reset_index()
            if (len(uf) > 0) and (uf.alt.max() / 1e3 >= max(alts) * 0.7):
                method = "intp" if uf.alt.max() / 1e3 > max(alts) else "extp"
                self.param[:, ix] = (
                    utils.interpolate_by_altitude(
                        np.array(uf.alt) / 1e3,
                        alts,
                        np.array(uf["nsall"]),
                        self.cfg.scale,
                        self.cfg.kind,
                        method=method,
                    )
                    * 1e-6
                )
        logger.info(f"Index: {intp_edens_xlim_index}")
        if intp_edens_xlim_index:
            for j in range(intp_edens_xlim_index, len(lats)):
                self.param[:, j] = self.param[:, intp_edens_xlim_index]
        if to_file:
            savemat(to_file, dict(ne=self.param))
        return self.param, alts

    def _fetch_profile_from_df(self, df, lat, lon, alts, dlat=0.2, dlon=0.2):
        """Return a single altitude-interpolated Ne profile at (lat, lon).

        Parameters
        ----------
        df : pandas.DataFrame
            Loaded grid/data frame with columns ``glat``, ``glon``, ``alt``,
            and ``nsall``.
        lat, lon : float
            Target geographic coordinates [°].
        alts : array-like
            Target altitude levels [km].
        dlat, dlon : float
            Spatial averaging half-width [°].

        Returns
        -------
        np.ndarray, shape (nalt,)
            Electron density in cm⁻³.
        """
        uf = df[
            (df.glat >= lat - dlat)
            & (df.glat <= lat + dlat)
            & (df.glon >= lon - dlon)
            & (df.glon <= lon + dlon)
        ]
        uf = uf.groupby("alt").mean().reset_index()
        p = np.zeros(len(alts), dtype=float)
        if (len(uf) > 0) and (uf.alt.max() / 1e3 >= max(alts) * 0.7):
            method = "intp" if uf.alt.max() / 1e3 > max(alts) else "extp"
            p = (
                utils.interpolate_by_altitude(
                    np.array(uf.alt) / 1e3,
                    alts,
                    np.array(uf["nsall"]),
                    self.cfg.scale,
                    self.cfg.kind,
                    method=method,
                )
                * 1e-6
            )
        return np.asarray(p, dtype=float)

    def fetch_dataset_3d(
        self,
        time: dt.datetime,
        lats,
        lons,
        alts,
        to_file: str = None,
        workers: int = 1,
        dlat: float = 0.2,
        dlon: float = 0.2,
    ):
        """
        Fetch 3D electron density cube on (lat, lon, alt) with optional
        point-wise parallelism.
        """
        lats = np.asarray(lats, dtype=float)
        lons = np.asarray(lons, dtype=float)
        alts = np.asarray(alts, dtype=float)
        i = np.argmin([np.abs((t - time).total_seconds()) for t in self.dates])
        file = self.files[i]
        df = self.load_data(file)
        df = df[df.alt >= 0]

        out = np.zeros((lats.size, lons.size, alts.size), dtype=float)
        ij = [(ii, jj) for ii in range(lats.size) for jj in range(lons.size)]

        def _job(ii, jj):
            p = self._fetch_profile_from_df(
                df, lats[ii], lons[jj], alts, dlat=dlat, dlon=dlon
            )
            return ii, jj, p

        n_workers = max(1, int(workers))
        if n_workers == 1:
            for ii, jj in ij:
                _, _, p = _job(ii, jj)
                out[ii, jj, :] = p
        else:
            with ThreadPoolExecutor(max_workers=n_workers) as ex:
                for ii, jj, p in ex.map(lambda t: _job(*t), ij):
                    out[ii, jj, :] = p

        self.param3d = out
        if to_file:
            savemat(to_file, dict(ne=self.param3d))
        return self.param3d, alts

    def load_from_file(self, to_file: str):
        """Load a previously saved electron density array from a ``.mat`` file.

        Parameters
        ----------
        to_file : str
            Path to the ``.mat`` file containing key ``ne``.

        Returns
        -------
        np.ndarray
            Electron density array as stored.
        """
        logger.info(f"Load from file {to_file.split('/')[-1]}")
        self.param = loadmat(to_file)["ne"]
        return self.param