GLOW 2D Local Polar Coordinate Transformation

glow2d.polar_model

Run GLOW model looking along heading from the current location and return the model output in (ZA, R) local coordinates where ZA is zenith angle in radians and R is distance in kilometers.

Parameters:

Name	Type	Description	Default
time	datetime	Datetime of GLOW calculation.	required
lat	Numeric	Latitude of starting location.	required
lon	Numeric	Longitude of starting location.	required
heading	Numeric	Heading (look direction).	required
altitude	Numeric	Altitude of local polar coordinate system origin in km above ASL. Must be < 100 km. Defaults to 0.	0
max_alt	Numeric	Maximum altitude where intersection is considered (km). Defaults to 1000, i.e. exobase.	1000
n_pts	int	Number of GEO coordinate angular grid points (i.e. number of GLOW runs), must be even and > 20. Defaults to 50.	50
n_bins	int	Number of energy bins. Defaults to 100.	100
n_alt	int	Number of altitude bins, must be > 100. Defaults to None, i.e. uses same number of bins as a single GLOW run.	None
with_prodloss	bool	Calculate production and loss parameters in local coordinates. Defaults to False.	False
n_threads	int	Number of threads for parallel GLOW runs. Set to None to use all system threads. Defaults to None.	None
full_output	bool	Returns only local coordinate GLOW output if False, and a tuple of local and GEO outputs if True. Defaults to False.	False
resamp	Numeric	Number of R and ZA points in local coordinate output. len(R) = len(alt_km) * resamp and len(ZA) = n_pts * resamp. Must be > 0.5. Defaults to 1.5.	1.5
show_progress	bool	Use TQDM to show progress of GLOW model calculations. Defaults to True.	True
kwargs	dict	Passed to glowpython.generic.	{}

Returns:

Name	Type	Description
iono	xarray.Dataset	Ionospheric parameters and brightnesses (with or without production and loss) in local coordinates. This is a reference and should not be modified.
	xarray.Dataset | tuple(xarray.Dataset, xarray.Dataset)	iono, bds (xarray.Dataset, xarray.Dataset): These values are returned only if full_output == True. Both are references and should not be modified.
	xarray.Dataset | tuple(xarray.Dataset, xarray.Dataset)	Ionospheric parameters and brightnesses (with or without production and loss) in local coordinates.
	xarray.Dataset | tuple(xarray.Dataset, xarray.Dataset)	Ionospheric parameters and brightnesses (with production and loss) in GEO coordinates.

Raises:

Type	Description
ValueError	Number of position bins can not be odd.
ValueError	Number of position bins can not be < 20.
ValueError	n_alt can not be < 100.
ValueError	Resampling can not be < 0.5.
ValueError	altitude must be in the range [0, 100].

Warns

ResourceWarning: Number of threads requested is more than available system threads.

Source code in glow2d/_glow2d.py

def polar_model(time: datetime, lat: Numeric, lon: Numeric, heading: Numeric, altitude: Numeric = 0, max_alt: Numeric = 1000, n_pts: int = 50, n_bins: int = 100, *, n_alt: int = None, with_prodloss: bool = False, n_threads: int = None, full_output: bool = False, resamp: Numeric = 1.5, show_progress: bool = True, **kwargs) -> xarray.Dataset | tuple(xarray.Dataset, xarray.Dataset):
    """Run GLOW model looking along heading from the current location and return the model output in
    (ZA, R) local coordinates where ZA is zenith angle in radians and R is distance in kilometers.

    Args:
        time (datetime): Datetime of GLOW calculation.
        lat (Numeric): Latitude of starting location.
        lon (Numeric): Longitude of starting location.
        heading (Numeric): Heading (look direction).
        altitude (Numeric, optional): Altitude of local polar coordinate system origin in km above ASL. Must be < 100 km. Defaults to 0.
        max_alt (Numeric, optional): Maximum altitude where intersection is considered (km). Defaults to 1000, i.e. exobase.
        n_pts (int, optional): Number of GEO coordinate angular grid points (i.e. number of GLOW runs), must be even and > 20. Defaults to 50.
        n_bins (int, optional): Number of energy bins. Defaults to 100.
        n_alt (int, optional): Number of altitude bins, must be > 100. Defaults to None, i.e. uses same number of bins as a single GLOW run.
        with_prodloss (bool, optional): Calculate production and loss parameters in local coordinates. Defaults to False.
        n_threads (int, optional):  Number of threads for parallel GLOW runs. Set to None to use all system threads. Defaults to None.
        full_output (bool, optional): Returns only local coordinate GLOW output if False, and a tuple of local and GEO outputs if True. Defaults to False.
        resamp (Numeric, optional): Number of R and ZA points in local coordinate output. ``len(R) = len(alt_km) * resamp`` and ``len(ZA) = n_pts * resamp``. Must be > 0.5. Defaults to 1.5.
        show_progress (bool, optional): Use TQDM to show progress of GLOW model calculations. Defaults to True.
        kwargs (dict, optional): Passed to `glowpython.generic`.

    Returns:
        iono (xarray.Dataset): Ionospheric parameters and brightnesses (with or without production and loss) in local coordinates. This is a reference and should not be modified.

        iono, bds (xarray.Dataset, xarray.Dataset): These values are returned only if ``full_output == True``. Both are references and should not be modified.

        - Ionospheric parameters and brightnesses (with or without production and loss) in local coordinates.
        - Ionospheric parameters and brightnesses (with production and loss) in GEO coordinates.


    Raises:
        ValueError: Number of position bins can not be odd.
        ValueError: Number of position bins can not be < 20.
        ValueError: n_alt can not be < 100.
        ValueError: Resampling can not be < 0.5.
        ValueError: altitude must be in the range [0, 100].

    Warns:
        ResourceWarning: Number of threads requested is more than available system threads.
    """
    grobj = glow2d_geo(time, lat, lon, heading, max_alt, n_pts, n_bins, n_alt=n_alt, uniformize_glow=True,
                       n_threads=n_threads, show_progress=show_progress, **kwargs)
    bds = grobj.run_model()
    grobj = glow2d_polar(bds, altitude, with_prodloss=with_prodloss, resamp=resamp)
    iono = grobj.transform_coord()
    if not full_output:
        return iono
    else:
        return (iono, bds)

Use GLOW model output evaluated on a 2D grid using glow2d_geo and convert it to a local ZA, R coordinate system at the origin location.

Source code in glow2d/_glow2d.py

class glow2d_polar:
    """Use GLOW model output evaluated on a 2D grid using `glow2d_geo` and convert it to a local ZA, R coordinate system at the origin location.
    """

    def __init__(self, bds: xarray.Dataset, altitude: Numeric = 0, *, with_prodloss: bool = False, resamp: Numeric = 1.5):
        """Create a GLOWRaycast object.

        Args:
            bds (xarray.Dataset): GLOW model evaluated on a 2D grid using `GLOW2D`.
            altitude (Numeric, optional): Altitude of local polar coordinate system origin in km above ASL. Must be < 100 km. Defaults to 0.
            with_prodloss (bool, optional): Calculate production and loss parameters in local coordinates. Defaults to False.
            resamp (Numeric, optional): Number of R and ZA points in local coordinate output. ``len(R) = len(alt_km) * resamp`` and ``len(ZA) = n_pts * resamp``. Must be > 0.5. Defaults to 1.5.

        Raises:
            ValueError: Resampling can not be < 0.5.
        """
        if resamp < 0.5:
            raise ValueError('Resampling can not be < 0.5.')
        if not (0 <= altitude <= 100):
            raise ValueError('Altitude can not be > 100 km.')
        self._resamp = resamp
        self._wprodloss = with_prodloss
        self._bds = bds.copy()
        self._iono = None
        self._r0 = altitude + EARTH_RADIUS

    def transform_coord(self) -> xarray.Dataset:
        """Run the coordinate transform to convert GLOW output from GEO to local coordinate system.

        Returns:
            xarray.Dataset: GLOW output in (ZA, r) coordinates. This is a reference and should not be modified.
        """
        if self._iono is not None:
            return self._iono
        tt, rr = self.get_local_coords(
            self._bds.angle.values, self._bds.alt_km.values + EARTH_RADIUS, r0=self._r0)  # get local coords from geocentric coords

        self._rmin, self._rmax = self._bds.alt_km.values.min(), rr.max()  # nearest and farthest local pts
        # highest and lowest za
        self._tmin, self._tmax = tt.min(), np.pi / 2  # 0, tt.max()
        self._nr_num = round(len(self._bds.alt_km.values) * self._resamp)  # resample
        self._nt_num = round(len(self._bds.angle.values) * self._resamp)   # resample
        outp_shape = (self._nt_num, self._nr_num)

        # ttidx = np.where(tt < 0)  # angle below horizon (LA < 0)
        # # get distribution of global -> local points in local grid
        # res = np.histogram2d(rr.flatten(), tt.flatten(), range=([rr.min(), rr.max()], [0, tt.max()]))
        # gd = resize(res[0], outp_shape, mode='edge')  # remap to right size
        # gd *= res[0].sum() / gd.sum()  # conserve sum of points
        # window_length = int(25 * self._resamp)  # smoothing window
        # window_length = window_length if window_length % 2 else window_length + 1  # must be odd
        # gd = savgol_filter(gd, window_length=window_length, polyorder=5, mode='nearest')  # smooth the distribution

        self._altkm = altkm = self._bds.alt_km.values  # store the altkm
        self._theta = theta = self._bds.angle.values  # store the angles
        rmin, rmax = self._rmin, self._rmax  # local names
        tmin, tmax = self._tmin, self._tmax
        self._nr = nr = np.linspace(
            rmin, rmax, self._nr_num, endpoint=True)  # local r
        self._nt = nt = np.linspace(
            tmin, tmax, self._nt_num, endpoint=True)  # local look angle
        # get meshgrid of the R, T coord system from regular r, za grid
        self._ntt, self._nrr = self.get_global_coords(nt, nr, r0=self._r0)
        # calculate jacobian
        jacobian = self.get_jacobian_glob2loc_glob(self._nrr, self._ntt, r0=self._r0)
        # convert to pixel coordinates
        self._ntt = self._ntt.flatten()  # flatten T, works as _global_from_local LUT
        self._nrr = self._nrr.flatten()  # flatten R, works as _global_from_local LUT
        self._ntt = (self._ntt - self._theta.min()) / \
            (self._theta.max() - self._theta.min()) * \
            len(self._theta)  # calculate equivalent index (pixel coord) from original T grid
        self._nrr = (self._nrr - self._altkm.min() - self._r0) / \
            (self._altkm.max() - self._altkm.min()) * \
            len(self._altkm)  # calculate equivalent index (pixel coord) from original R (alt_km) grid
        # start transformation
        data_vars = {}
        bds = self._bds
        coord_wavelength = bds.wavelength.values  # wl axis
        coord_state = bds.state.values  # state axis
        coord_energy = bds.energy.values  # egrid
        bds_attr = bds.attrs  # attrs
        single_keys = ['Tn',
                       'pedersen',
                       'hall',
                       'Te',
                       'Ti']  # (angle, alt_km) vars
        density_keys = [
            'O',
            'N2',
            'O2',
            'NO',
            'NeIn',
            'NeOut',
            'ionrate',
            'O+',
            'O2+',
            'NO+',
            'N2D',
        ]
        state_keys = [
            'production',
            'loss',
            'excitedDensity'
        ]  # (angle, alt_km, state) vars
        # start = perf_counter_ns()
        # map all the single key types from (angle, alt_km) -> (la, r)
        for key in single_keys:
            inp = self._bds[key].values
            inp[np.where(np.isnan(inp))] = 0
            out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape, mode='nearest')
            out[np.where(out < 0)] = 0
            # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape)
            data_vars[key] = (('za', 'r'), out)
        for key in density_keys:
            inp = self._bds[key].values
            inp[np.where(np.isnan(inp))] = 0
            out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape, mode='nearest') / jacobian
            out[np.where(out < 0)] = 0
            # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape)
            data_vars[key] = (('za', 'r'), out)
        # end = perf_counter_ns()
        # print('Single_key conversion: %.3f us'%((end - start)*1e-3))
        # start = perf_counter_ns()
        # dataset of (angle, alt_km) vars
        iono = xarray.Dataset(data_vars=data_vars, coords={
            'za': nt, 'r': nr})
        # end = perf_counter_ns()
        # print('Single_key dataset: %.3f us'%((end - start)*1e-3))
        # start = perf_counter_ns()
        ver = []
        # map all the wavelength data from (angle, alt_km, wavelength) -> (la, r, wavelength)
        for key in coord_wavelength:
            inp = bds['ver'].loc[dict(wavelength=key)].values
            inp[np.where(np.isnan(inp))] = 0
            # scaled by point distribution because flux is conserved, not brightness
            # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape, mode='nearest') * gd
            out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape, mode='nearest') / jacobian
            out[np.where(out < 0)] = 0
            # inp[ttidx] = 0
            # inpsum = inp.sum()  # sum of input for valid angles
            # outpsum = out.sum()  # sum of output
            # out = out * (inpsum / outpsum)  # scale the sum to conserve total flux
            ver.append(out.T)
        # end = perf_counter_ns()
        # print('VER eval: %.3f us'%((end - start)*1e-3))
        # start = perf_counter_ns()
        ver = np.asarray(ver).T
        ver = xr.DataArray(
            ver,
            coords={'za': nt, 'r': nr,
                    'wavelength': coord_wavelength},
            dims=['za', 'r', 'wavelength'],
            name='ver'
        )  # create wl dataset
        # end = perf_counter_ns()
        # print('VER dataset: %.3f us'%((end - start)*1e-3))
        # start = perf_counter_ns()
        if self._wprodloss:
            d = {}
            for key in state_keys:  # for each var with (angle, alt_km, state)
                res = []

                def convert_state_stuff(st):
                    inp = bds[key].loc[dict(state=st)].values
                    inp[np.where(np.isnan(inp))] = 0
                    out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape)
                    out[np.where(out < 0)] = 0
                    if key in ('production', 'excitedDensity'):
                        out /= jacobian
                    # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape)
                    return out.T
                res = list(map(convert_state_stuff, coord_state))
                res = np.asarray(res).T
                d[key] = (('za', 'r', 'state'), res)
            # end = perf_counter_ns()
            # print('Prod_Loss Eval: %.3f us'%((end - start)*1e-3))
            # start = perf_counter_ns()
            prodloss = xarray.Dataset(
                data_vars=d,
                coords={'za': nt, 'r': nr, 'state': coord_state}
            )  # calculate (angle, alt_km, state) -> (la, r, state) dataset
        else:
            prodloss = xarray.Dataset()
        # end = perf_counter_ns()
        # print('Prod_Loss DS: %.3f us'%((end - start)*1e-3))
        ## EGrid conversion (angle, energy) -> (r, energy) ##
        # EGrid is avaliable really at (angle, alt_km = 0, energy)
        # So we get local coords for (angle, R=R0)
        # we discard the angle information because it is meaningless, EGrid is spatial
        # start = perf_counter_ns()
        _rr, _ = self.get_local_coords(
            bds.angle.values, np.ones(bds.angle.values.shape)*self._r0, r0=self._r0)
        _rr = rr[:, 0]  # spatial EGrid
        d = []
        for en in coord_energy:  # for each energy
            inp = bds['precip'].loc[dict(energy=en)].values
            # interpolate to appropriate energy grid
            out = np.interp(nr, _rr, inp)
            d.append(out)
        d = np.asarray(d).T
        precip = xarray.Dataset({'precip': (('r', 'energy'), d)}, coords={
            'r': nr, 'energy': coord_energy})
        # end = perf_counter_ns()
        # print('Precip interp and ds: %.3f us'%((end - start)*1e-3))

        # start = perf_counter_ns()
        iono = xr.merge((iono, ver, prodloss, precip))  # merge all datasets
        bds_attr['altitude'] = {'values': self._r0 - EARTH_RADIUS, 'units': 'km',
                                'description': 'Altitude of local polar coordinate origin ASL'}
        iono.attrs.update(bds_attr)  # copy original attrs

        _ = list(map(lambda x: iono[x].attrs.update(bds[x].attrs), tuple(iono.data_vars.keys())))  # update attrs from bds

        unit_desc_dict = {
            'za': ('radians', 'Zenith angle'),
            'r': ('km', 'Radial distance in km')
        }
        _ = list(map(lambda x: iono[x].attrs.update(
            {'units': unit_desc_dict[x][0], 'description': unit_desc_dict[x][1]}), unit_desc_dict.keys()))
        # end = perf_counter_ns()
        # print('Merging: %.3f us'%((end - start)*1e-3))
        self._iono = iono
        return iono

    @staticmethod
    def get_emission(iono: xarray.Dataset, feature: str = '5577', za_min: Numeric | Iterable = np.deg2rad(20), za_max: Numeric | Iterable = np.deg2rad(25), num_zapts: int = 10, *, rmin: Numeric = None, rmax: Numeric = None, num_rpts: int = 100) -> float | np.ndarray:
        """Calculate number of photons per azimuth angle (radians) per unit area per second coming from a region of (`rmin`, `rmax`, `za_min`, `za_max`).

        Args:
            iono (xarray.Dataset, optional): GLOW model output in local polar coordinates calculated using `glow2d.glow2d_polar.transform_coord`.
            feature (str, optional):GLOW emission feature. Defaults to '5577'.
            za_min (Numeric | Iterable, optional): Minimum zenith angle. Defaults to np.deg2rad(20).
            za_max (Numeric | Iterable, optional): Maximum zenith angle. Defaults to np.deg2rad(25).
            num_zapts (int, optional): Number of points to interpolate to. Defaults to 10.
            rmin (Numeric, optional): Minimum distance. Defaults to None.
            rmax (Numeric, optional): Maximum distance. Defaults to None.
            num_rpts (int, optional): Number of distance points. The default is used only if minimum or maximum distance is not None. Defaults to 100.

        Raises:
            ValueError: iono is not an xarray.Dataset.
            ValueError: ZA min and max arrays must be of the same dimension.
            ValueError: ZA min not between 0 deg and 90 deg.
            ValueError: ZA max is not between 0 deg and 90 deg.
            ValueError: ZA min > ZA max.
            ValueError: Selected feature is invalid.

        Returns:
            float | np.ndarray: Number of photons/rad/cm^2/s
        """
        if iono is None or not isinstance(iono, xarray.Dataset):
            raise ValueError('iono is not an xarray.Dataset.')
        if isinstance(za_min, Iterable) or isinstance(za_max, Iterable):
            if len(za_min) != len(za_max):
                raise ValueError('ZA min and max arrays must be of the same dimension.')
            callable = partial(glow2d_polar.get_emission, iono=iono, feature=feature,
                               num_zapts=num_zapts, rmin=rmin, rmax=rmax, num_rpts=num_rpts)
            out = list(map(lambda idx: callable(za_min=za_min[idx], za_max=za_max[idx]), range(len(za_min))))
            return np.asarray(out, dtype=float)
        if not (0 <= za_min <= np.deg2rad(90)):
            raise ValueError('ZA must be between 0 deg and 90 deg')
        if not (0 <= za_max <= np.deg2rad(90)):
            raise ValueError('ZA must be between 0 deg and 90 deg')
        if za_min > za_max:
            raise ValueError('ZA min > ZA max')
        if feature not in iono.wavelength.values:
            raise ValueError('Feature %s is invalid. Valid features: %s' % (feature, str(iono.wavelength.values)))
        za: np.ndarray = iono.za.values
        zaxis = iono.za.values
        r: np.ndarray = iono.r.values
        rr = iono.r.values

        if za_min is not None or za_max is not None:
            if (za_min == 0) and (za_max == np.deg2rad(90)):
                pass
            else:
                za_min = za.min() if za_min is None else za_min
                za_max = za.max() if za_max is None else za_max
                zaxis = np.linspace(za_min, za_max, num_zapts, endpoint=True)

        if rmin is not None or rmax is not None:
            rmin = r.min() if rmin is None else rmin
            rmax = r.max() if rmax is None else rmax
            rr = np.linspace(rmin, rmax, num_rpts, endpoint=True)

        ver = iono.ver.loc[dict(wavelength=feature)].values
        ver = interp2d(r, za, ver)(rr, zaxis)  # interpolate to integration axes

        ver = ver*np.sin(zaxis[:, None])  # integration is VER * sin(phi) * d(phi) * d(r)
        return simps(simps(ver.T, zaxis), rr * 1e5)  # do the double integral

    # get global coord index from local coord index, implemented as LUT
    def _global_from_local(self, pt: tuple(int, int)) -> tuple(float, float):
        # if not self.firstrun % 10000:
        #     print('Input:', pt)
        tl, rl = pt  # pixel coord
        # rl = (rl * (self._rmax - self._rmin) / self._nr_num) + self._rmin # real coord
        # tl = (tl * (self._tmax - self._tmin) / self._nt_num) + self._tmin
        # rl = self._nr[rl]
        # tl = self._nt[tl]
        # t, r = self._get_global_coords(tl, rl)
        # # if not self.firstrun % 10000:
        # #     print((rl, tl), '->', (r, t), ':', (self._altkm.min(), self._altkm.max()))
        # r = (r - self._altkm.min() - EARTH_RADIUS) / (self._altkm.max() - self._altkm.min()) * len(self._altkm)
        # t = (t - self._theta.min()) / (self._theta.max() - self._theta.min()) * len(self._theta)
        # if not self.firstrun % 10000:
        #     print((float(r), float(t)))
        return (float(self._ntt[tl*self._nr_num + rl]), float(self._nrr[tl*self._nr_num + rl]))

    @staticmethod
    def get_global_coords(t: np.ndarray | Numeric, r: np.ndarray | Numeric, r0: Numeric = EARTH_RADIUS, meshgrid: bool = True) -> tuple(np.ndarray, np.ndarray):
        """Get GEO coordinates from local coordinates.

        $$
            R = \\sqrt{\\left\\{ (r\\cos{\\phi} + R_0)^2 + r^2\\sin{\\phi}^2 \\right\\}}, \\\\
            \\theta = \\arctan{\\frac{r\\sin{\\phi}}{r\\cos{\\phi} + R_0}}
        $$

        Args:
            t (np.ndarray | Numeric): Angles in radians.
            r (np.ndarray | Numeric): Distance in km.
            r0 (Numeric, optional): Distance to origin. Defaults to geopy.distance.EARTH_RADIUS.
            meshgrid (bool, optional): Optionally convert 1-D inputs to a meshgrid. Defaults to True.

        Raises:
            ValueError: ``r`` and ``t`` does not have the same dimensions
            TypeError: ``r`` and ``t`` are not ``numpy.ndarray``.

        Returns:
            (np.ndarray, np.ndarray): (angles, distances) in GEO coordinates.
        """
        if isinstance(r, np.ndarray) and isinstance(t, np.ndarray):  # if array
            if r.ndim != t.ndim:  # if dims don't match get out
                raise ValueError('r and t does not have the same dimensions')
            if r.ndim == 1 and meshgrid:
                _r, _t = np.meshgrid(r, t)
            elif r.ndim == 1 and not meshgrid:
                _r, _t = r, t
            else:
                _r, _t = r.copy(), t.copy()  # already a meshgrid?
                r = _r[0]
                t = _t[:, 0]
        elif isinstance(r, Numeric) and isinstance(t, Numeric):  # floats
            _r = np.atleast_1d(r)
            _t = np.atleast_1d(t)
        else:
            raise TypeError('r and t must be np.ndarray.')
        # _t = np.pi/2 - _t
        rr = np.sqrt((_r*np.cos(_t) + r0)**2 +
                     (_r*np.sin(_t))**2)  # r, la to R, T
        tt = np.arctan2(_r*np.sin(_t), _r*np.cos(_t) + r0)
        return tt, rr

    @staticmethod
    def get_local_coords(t: np.ndarray | Numeric, r: np.ndarray | Numeric, r0: Numeric = EARTH_RADIUS, meshgrid: bool = True) -> tuple(np.ndarray, np.ndarray):
        """Get local coordinates from GEO coordinates.

        $$
            r = \\sqrt{\\left\\{ (R\\cos{\\theta} - R_0)^2 + R^2\\sin{\\theta}^2 \\right\\}}, \\\\
            \\phi = \\arctan{\\frac{R\\sin{\\theta}}{R\\cos{\\theta} - R_0}}
        $$

        Args:
            t (np.ndarray | Numeric): Angles in radians.
            r (np.ndarray | Numeric): Distance in km.
            r0 (Numeric, optional): Distance to origin. Defaults to geopy.distance.EARTH_RADIUS.
            meshgrid (bool, optional): Optionally convert 1-D inputs to a meshgrid. Defaults to True.

        Raises:
            ValueError: ``r`` and ``t`` does not have the same dimensions
            TypeError: ``r`` and ``t`` are not ``numpy.ndarray``.

        Returns:
            (np.ndarray, np.ndarray): (angles, distances) in local coordinates.
        """
        if isinstance(r, np.ndarray) and isinstance(t, np.ndarray):
            if r.ndim != t.ndim:  # if dims don't match get out
                raise ValueError('r and t does not have the same dimensions')
            if r.ndim == 1 and meshgrid:
                _r, _t = np.meshgrid(r, t)
            elif r.ndim == 1 and not meshgrid:
                _r, _t = r, t
            else:
                _r, _t = r.copy(), t.copy()  # already a meshgrid?
                r = _r[0]
                t = _t[:, 0]
        elif isinstance(r, Numeric) and isinstance(t, Numeric):
            _r = np.atleast_1d(r)
            _t = np.atleast_1d(t)
        else:
            raise TypeError('r and t must be np.ndarray.')
        rr = np.sqrt((_r*np.cos(_t) - r0)**2 +
                     (_r*np.sin(_t))**2)  # R, T to r, la
        tt = np.arctan2(_r*np.sin(_t), _r*np.cos(_t) - r0)
        return tt, rr

    @staticmethod
    def get_jacobian_glob2loc_loc(r: np.ndarray, t: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
        """Jacobian \\(|J_{R\\rightarrow r}|\\) for global to local coordinate transform, evaluated at points in local coordinate.

        $$
            |J_{R\\rightarrow r}| = \\frac{R}{r^3}\\left(R^2 + R_0^2 - 2 R R_0 \\cos{\\theta}\\right)
        $$

        Args:
            r (np.ndarray): 2-dimensional array of r.
            t (np.ndarray): 2-dimensional array of phi.
            r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

        Raises:
            ValueError: Dimension of inputs must be 2.

        Returns:
            np.ndarray: Jacobian evaluated at points.
        """
        if r.ndim != 2 or t.ndim != 2:
            raise ValueError('Dimension of inputs must be 2.')
        gt, gr = glow2d_polar.get_global_coords(t, r, r0=r0)
        jac = (gr / (r**3)) * ((gr**2) + (r0**2) - (2*gr*r0*np.cos(gt)))
        return jac

    @staticmethod
    def get_jacobian_loc2glob_loc(r: np.ndarray, t: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
        """Jacobian \\(|J_{r\\rightarrow R}|\\) for local to global coordinate transform, evaluated at points in local coordinate.

        $$
            |J_{r\\rightarrow R}| = \\frac{r}{R^3}\\left(r^2 + R_0^2 + 2 r R_0 \\cos{\\phi}\\right)
        $$

        Args:
            r (np.ndarray): 2-dimensional array of r.
            t (np.ndarray): 2-dimensional array of phi.
            r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

        Raises:
            ValueError: Dimension of inputs must be 2.

        Returns:
            np.ndarray: Jacobian evaluated at points.
        """
        if r.ndim != 2 or t.ndim != 2:
            raise ValueError('Dimension of inputs must be 2')
        gt, gr = glow2d_polar.get_global_coords(t, r, r0=r0)
        jac = (r/(gr**3))*((r**2) + (r0**2) + (2*r*r0*np.cos(t)))
        return jac

    @staticmethod
    def get_jacobian_glob2loc_glob(gr: np.ndarray, gt: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
        """Jacobian determinant \\(|J_{R\\rightarrow r}|\\) for global to local coordinate transform, evaluated at points in global coordinate.

        $$
            |J_{R\\rightarrow r}| = \\frac{R}{r^3}\\left(R^2 + R_0^2 - 2 R R_0 \\cos{\\theta}\\right)
        $$

        Args:
            gr (np.ndarray): 2-dimensional array of R.
            gt (np.ndarray): 2-dimensional array of Theta.
            r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

        Raises:
            ValueError: Dimension of inputs must be 2.

        Returns:
            np.ndarray: Jacobian evaluated at points.
        """
        if gr.ndim != 2 or gt.ndim != 2:
            raise ValueError('Dimension of inputs must be 2')
        t, r = glow2d_polar.get_local_coords(gt, gr, r0=r0)
        jac = (gr / (r**3)) * ((gr**2) + (r0**2) - (2*gr*r0*np.cos(gt)))
        return jac

    @staticmethod
    def get_jacobian_loc2glob_glob(gr: np.ndarray, gt: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
        """Jacobian \\(|J_{r\\rightarrow R}|\\) for global to local coordinate transform, evaluated at points in local coordinate.

        $$
            |J_{r\\rightarrow R}| = \\frac{r}{R^3}\\left(r^2 + R_0^2 + 2 r R_0 \\cos{\\phi}\\right)
        $$

        Args:
            gr (np.ndarray): 2-dimensional array of R.
            gt (np.ndarray): 2-dimensional array of Theta.
            r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

        Raises:
            ValueError: Dimension of inputs must be 2.

        Returns:
            np.ndarray: Jacobian evaluated at points.
        """
        if gr.ndim != 2 or gt.ndim != 2:
            raise ValueError('Dimension of inputs must be 2')
        t, r = glow2d_polar.get_local_coords(gt, gr, r0=EARTH_RADIUS)
        jac = (r/(gr**3))*((r**2) + (r0**2) + (2*r*r0*np.cos(t)))
        return jac

__init__(bds, altitude=0, *, with_prodloss=False, resamp=1.5)

Create a GLOWRaycast object.

Parameters:

Name	Type	Description	Default
bds	xarray.Dataset	GLOW model evaluated on a 2D grid using GLOW2D.	required
altitude	Numeric	Altitude of local polar coordinate system origin in km above ASL. Must be < 100 km. Defaults to 0.	0
with_prodloss	bool	Calculate production and loss parameters in local coordinates. Defaults to False.	False
resamp	Numeric	Number of R and ZA points in local coordinate output. len(R) = len(alt_km) * resamp and len(ZA) = n_pts * resamp. Must be > 0.5. Defaults to 1.5.	1.5

Raises:

Type	Description
ValueError	Resampling can not be < 0.5.

Source code in glow2d/_glow2d.py

def __init__(self, bds: xarray.Dataset, altitude: Numeric = 0, *, with_prodloss: bool = False, resamp: Numeric = 1.5):
    """Create a GLOWRaycast object.

    Args:
        bds (xarray.Dataset): GLOW model evaluated on a 2D grid using `GLOW2D`.
        altitude (Numeric, optional): Altitude of local polar coordinate system origin in km above ASL. Must be < 100 km. Defaults to 0.
        with_prodloss (bool, optional): Calculate production and loss parameters in local coordinates. Defaults to False.
        resamp (Numeric, optional): Number of R and ZA points in local coordinate output. ``len(R) = len(alt_km) * resamp`` and ``len(ZA) = n_pts * resamp``. Must be > 0.5. Defaults to 1.5.

    Raises:
        ValueError: Resampling can not be < 0.5.
    """
    if resamp < 0.5:
        raise ValueError('Resampling can not be < 0.5.')
    if not (0 <= altitude <= 100):
        raise ValueError('Altitude can not be > 100 km.')
    self._resamp = resamp
    self._wprodloss = with_prodloss
    self._bds = bds.copy()
    self._iono = None
    self._r0 = altitude + EARTH_RADIUS

get_emission(iono, feature='5577', za_min=np.deg2rad(20), za_max=np.deg2rad(25), num_zapts=10, *, rmin=None, rmax=None, num_rpts=100) staticmethod

Calculate number of photons per azimuth angle (radians) per unit area per second coming from a region of (rmin, rmax, za_min, za_max).

Parameters:

Name	Type	Description	Default
iono	xarray.Dataset	GLOW model output in local polar coordinates calculated using glow2d.glow2d_polar.transform_coord.	required
feature	str	GLOW emission feature. Defaults to '5577'.	'5577'
za_min	Numeric | Iterable	Minimum zenith angle. Defaults to np.deg2rad(20).	np.deg2rad(20)
za_max	Numeric | Iterable	Maximum zenith angle. Defaults to np.deg2rad(25).	np.deg2rad(25)
num_zapts	int	Number of points to interpolate to. Defaults to 10.	10
rmin	Numeric	Minimum distance. Defaults to None.	None
rmax	Numeric	Maximum distance. Defaults to None.	None
num_rpts	int	Number of distance points. The default is used only if minimum or maximum distance is not None. Defaults to 100.	100

Raises:

Type	Description
ValueError	iono is not an xarray.Dataset.
ValueError	ZA min and max arrays must be of the same dimension.
ValueError	ZA min not between 0 deg and 90 deg.
ValueError	ZA max is not between 0 deg and 90 deg.
ValueError	ZA min > ZA max.
ValueError	Selected feature is invalid.

Returns:

Type	Description
float | np.ndarray	float | np.ndarray: Number of photons/rad/cm^2/s

Source code in glow2d/_glow2d.py

@staticmethod
def get_emission(iono: xarray.Dataset, feature: str = '5577', za_min: Numeric | Iterable = np.deg2rad(20), za_max: Numeric | Iterable = np.deg2rad(25), num_zapts: int = 10, *, rmin: Numeric = None, rmax: Numeric = None, num_rpts: int = 100) -> float | np.ndarray:
    """Calculate number of photons per azimuth angle (radians) per unit area per second coming from a region of (`rmin`, `rmax`, `za_min`, `za_max`).

    Args:
        iono (xarray.Dataset, optional): GLOW model output in local polar coordinates calculated using `glow2d.glow2d_polar.transform_coord`.
        feature (str, optional):GLOW emission feature. Defaults to '5577'.
        za_min (Numeric | Iterable, optional): Minimum zenith angle. Defaults to np.deg2rad(20).
        za_max (Numeric | Iterable, optional): Maximum zenith angle. Defaults to np.deg2rad(25).
        num_zapts (int, optional): Number of points to interpolate to. Defaults to 10.
        rmin (Numeric, optional): Minimum distance. Defaults to None.
        rmax (Numeric, optional): Maximum distance. Defaults to None.
        num_rpts (int, optional): Number of distance points. The default is used only if minimum or maximum distance is not None. Defaults to 100.

    Raises:
        ValueError: iono is not an xarray.Dataset.
        ValueError: ZA min and max arrays must be of the same dimension.
        ValueError: ZA min not between 0 deg and 90 deg.
        ValueError: ZA max is not between 0 deg and 90 deg.
        ValueError: ZA min > ZA max.
        ValueError: Selected feature is invalid.

    Returns:
        float | np.ndarray: Number of photons/rad/cm^2/s
    """
    if iono is None or not isinstance(iono, xarray.Dataset):
        raise ValueError('iono is not an xarray.Dataset.')
    if isinstance(za_min, Iterable) or isinstance(za_max, Iterable):
        if len(za_min) != len(za_max):
            raise ValueError('ZA min and max arrays must be of the same dimension.')
        callable = partial(glow2d_polar.get_emission, iono=iono, feature=feature,
                           num_zapts=num_zapts, rmin=rmin, rmax=rmax, num_rpts=num_rpts)
        out = list(map(lambda idx: callable(za_min=za_min[idx], za_max=za_max[idx]), range(len(za_min))))
        return np.asarray(out, dtype=float)
    if not (0 <= za_min <= np.deg2rad(90)):
        raise ValueError('ZA must be between 0 deg and 90 deg')
    if not (0 <= za_max <= np.deg2rad(90)):
        raise ValueError('ZA must be between 0 deg and 90 deg')
    if za_min > za_max:
        raise ValueError('ZA min > ZA max')
    if feature not in iono.wavelength.values:
        raise ValueError('Feature %s is invalid. Valid features: %s' % (feature, str(iono.wavelength.values)))
    za: np.ndarray = iono.za.values
    zaxis = iono.za.values
    r: np.ndarray = iono.r.values
    rr = iono.r.values

    if za_min is not None or za_max is not None:
        if (za_min == 0) and (za_max == np.deg2rad(90)):
            pass
        else:
            za_min = za.min() if za_min is None else za_min
            za_max = za.max() if za_max is None else za_max
            zaxis = np.linspace(za_min, za_max, num_zapts, endpoint=True)

    if rmin is not None or rmax is not None:
        rmin = r.min() if rmin is None else rmin
        rmax = r.max() if rmax is None else rmax
        rr = np.linspace(rmin, rmax, num_rpts, endpoint=True)

    ver = iono.ver.loc[dict(wavelength=feature)].values
    ver = interp2d(r, za, ver)(rr, zaxis)  # interpolate to integration axes

    ver = ver*np.sin(zaxis[:, None])  # integration is VER * sin(phi) * d(phi) * d(r)
    return simps(simps(ver.T, zaxis), rr * 1e5)  # do the double integral

get_global_coords(t, r, r0=EARTH_RADIUS, meshgrid=True) staticmethod

Get GEO coordinates from local coordinates.

$$R = \sqrt{\left\{ (r\cos{\phi} + R_0)^2 + r^2\sin{\phi}^2 \right\}}, \\ \theta = \arctan{\frac{r\sin{\phi}}{r\cos{\phi} + R_0}}$$

Parameters:

Name	Type	Description	Default
t	np.ndarray | Numeric	Angles in radians.	required
r	np.ndarray | Numeric	Distance in km.	required
r0	Numeric	Distance to origin. Defaults to geopy.distance.EARTH_RADIUS.	EARTH_RADIUS
meshgrid	bool	Optionally convert 1-D inputs to a meshgrid. Defaults to True.	True

Raises:

Type	Description
ValueError	r and t does not have the same dimensions
TypeError	r and t are not numpy.ndarray.

Returns:

Type	Description
np.ndarray, np.ndarray	(angles, distances) in GEO coordinates.

Source code in glow2d/_glow2d.py

@staticmethod
def get_global_coords(t: np.ndarray | Numeric, r: np.ndarray | Numeric, r0: Numeric = EARTH_RADIUS, meshgrid: bool = True) -> tuple(np.ndarray, np.ndarray):
    """Get GEO coordinates from local coordinates.

    $$
        R = \\sqrt{\\left\\{ (r\\cos{\\phi} + R_0)^2 + r^2\\sin{\\phi}^2 \\right\\}}, \\\\
        \\theta = \\arctan{\\frac{r\\sin{\\phi}}{r\\cos{\\phi} + R_0}}
    $$

    Args:
        t (np.ndarray | Numeric): Angles in radians.
        r (np.ndarray | Numeric): Distance in km.
        r0 (Numeric, optional): Distance to origin. Defaults to geopy.distance.EARTH_RADIUS.
        meshgrid (bool, optional): Optionally convert 1-D inputs to a meshgrid. Defaults to True.

    Raises:
        ValueError: ``r`` and ``t`` does not have the same dimensions
        TypeError: ``r`` and ``t`` are not ``numpy.ndarray``.

    Returns:
        (np.ndarray, np.ndarray): (angles, distances) in GEO coordinates.
    """
    if isinstance(r, np.ndarray) and isinstance(t, np.ndarray):  # if array
        if r.ndim != t.ndim:  # if dims don't match get out
            raise ValueError('r and t does not have the same dimensions')
        if r.ndim == 1 and meshgrid:
            _r, _t = np.meshgrid(r, t)
        elif r.ndim == 1 and not meshgrid:
            _r, _t = r, t
        else:
            _r, _t = r.copy(), t.copy()  # already a meshgrid?
            r = _r[0]
            t = _t[:, 0]
    elif isinstance(r, Numeric) and isinstance(t, Numeric):  # floats
        _r = np.atleast_1d(r)
        _t = np.atleast_1d(t)
    else:
        raise TypeError('r and t must be np.ndarray.')
    # _t = np.pi/2 - _t
    rr = np.sqrt((_r*np.cos(_t) + r0)**2 +
                 (_r*np.sin(_t))**2)  # r, la to R, T
    tt = np.arctan2(_r*np.sin(_t), _r*np.cos(_t) + r0)
    return tt, rr

get_jacobian_glob2loc_glob(gr, gt, r0=EARTH_RADIUS) staticmethod

Jacobian determinant $|J_{R\rightarrow r}|$ for global to local coordinate transform, evaluated at points in global coordinate.

$$|J_{R\rightarrow r}| = \frac{R}{r^3}\left(R^2 + R_0^2 - 2 R R_0 \cos{\theta}\right)$$

Parameters:

Name	Type	Description	Default
gr	np.ndarray	2-dimensional array of R.	required
gt	np.ndarray	2-dimensional array of Theta.	required
r0	Numeric	Coordinate transform offset. Defaults to EARTH_RADIUS.	EARTH_RADIUS

Raises:

Type	Description
ValueError	Dimension of inputs must be 2.

Returns:

Type	Description
np.ndarray	np.ndarray: Jacobian evaluated at points.

Source code in glow2d/_glow2d.py

@staticmethod
def get_jacobian_glob2loc_glob(gr: np.ndarray, gt: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
    """Jacobian determinant \\(|J_{R\\rightarrow r}|\\) for global to local coordinate transform, evaluated at points in global coordinate.

    $$
        |J_{R\\rightarrow r}| = \\frac{R}{r^3}\\left(R^2 + R_0^2 - 2 R R_0 \\cos{\\theta}\\right)
    $$

    Args:
        gr (np.ndarray): 2-dimensional array of R.
        gt (np.ndarray): 2-dimensional array of Theta.
        r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

    Raises:
        ValueError: Dimension of inputs must be 2.

    Returns:
        np.ndarray: Jacobian evaluated at points.
    """
    if gr.ndim != 2 or gt.ndim != 2:
        raise ValueError('Dimension of inputs must be 2')
    t, r = glow2d_polar.get_local_coords(gt, gr, r0=r0)
    jac = (gr / (r**3)) * ((gr**2) + (r0**2) - (2*gr*r0*np.cos(gt)))
    return jac

get_jacobian_glob2loc_loc(r, t, r0=EARTH_RADIUS) staticmethod

Jacobian $|J_{R\rightarrow r}|$ for global to local coordinate transform, evaluated at points in local coordinate.

$$|J_{R\rightarrow r}| = \frac{R}{r^3}\left(R^2 + R_0^2 - 2 R R_0 \cos{\theta}\right)$$

Parameters:

Name	Type	Description	Default
r	np.ndarray	2-dimensional array of r.	required
t	np.ndarray	2-dimensional array of phi.	required
r0	Numeric	Coordinate transform offset. Defaults to EARTH_RADIUS.	EARTH_RADIUS

Raises:

Type	Description
ValueError	Dimension of inputs must be 2.

Returns:

Type	Description
np.ndarray	np.ndarray: Jacobian evaluated at points.

Source code in glow2d/_glow2d.py

@staticmethod
def get_jacobian_glob2loc_loc(r: np.ndarray, t: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
    """Jacobian \\(|J_{R\\rightarrow r}|\\) for global to local coordinate transform, evaluated at points in local coordinate.

    $$
        |J_{R\\rightarrow r}| = \\frac{R}{r^3}\\left(R^2 + R_0^2 - 2 R R_0 \\cos{\\theta}\\right)
    $$

    Args:
        r (np.ndarray): 2-dimensional array of r.
        t (np.ndarray): 2-dimensional array of phi.
        r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

    Raises:
        ValueError: Dimension of inputs must be 2.

    Returns:
        np.ndarray: Jacobian evaluated at points.
    """
    if r.ndim != 2 or t.ndim != 2:
        raise ValueError('Dimension of inputs must be 2.')
    gt, gr = glow2d_polar.get_global_coords(t, r, r0=r0)
    jac = (gr / (r**3)) * ((gr**2) + (r0**2) - (2*gr*r0*np.cos(gt)))
    return jac

get_jacobian_loc2glob_glob(gr, gt, r0=EARTH_RADIUS) staticmethod

Jacobian $|J_{r\rightarrow R}|$ for global to local coordinate transform, evaluated at points in local coordinate.

$$|J_{r\rightarrow R}| = \frac{r}{R^3}\left(r^2 + R_0^2 + 2 r R_0 \cos{\phi}\right)$$

Parameters:

Name	Type	Description	Default
gr	np.ndarray	2-dimensional array of R.	required
gt	np.ndarray	2-dimensional array of Theta.	required
r0	Numeric	Coordinate transform offset. Defaults to EARTH_RADIUS.	EARTH_RADIUS

Raises:

Type	Description
ValueError	Dimension of inputs must be 2.

Returns:

Type	Description
np.ndarray	np.ndarray: Jacobian evaluated at points.

Source code in glow2d/_glow2d.py

@staticmethod
def get_jacobian_loc2glob_glob(gr: np.ndarray, gt: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
    """Jacobian \\(|J_{r\\rightarrow R}|\\) for global to local coordinate transform, evaluated at points in local coordinate.

    $$
        |J_{r\\rightarrow R}| = \\frac{r}{R^3}\\left(r^2 + R_0^2 + 2 r R_0 \\cos{\\phi}\\right)
    $$

    Args:
        gr (np.ndarray): 2-dimensional array of R.
        gt (np.ndarray): 2-dimensional array of Theta.
        r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

    Raises:
        ValueError: Dimension of inputs must be 2.

    Returns:
        np.ndarray: Jacobian evaluated at points.
    """
    if gr.ndim != 2 or gt.ndim != 2:
        raise ValueError('Dimension of inputs must be 2')
    t, r = glow2d_polar.get_local_coords(gt, gr, r0=EARTH_RADIUS)
    jac = (r/(gr**3))*((r**2) + (r0**2) + (2*r*r0*np.cos(t)))
    return jac

get_jacobian_loc2glob_loc(r, t, r0=EARTH_RADIUS) staticmethod

Jacobian $|J_{r\rightarrow R}|$ for local to global coordinate transform, evaluated at points in local coordinate.

$$|J_{r\rightarrow R}| = \frac{r}{R^3}\left(r^2 + R_0^2 + 2 r R_0 \cos{\phi}\right)$$

Parameters:

Name	Type	Description	Default
r	np.ndarray	2-dimensional array of r.	required
t	np.ndarray	2-dimensional array of phi.	required
r0	Numeric	Coordinate transform offset. Defaults to EARTH_RADIUS.	EARTH_RADIUS

Raises:

Type	Description
ValueError	Dimension of inputs must be 2.

Returns:

Type	Description
np.ndarray	np.ndarray: Jacobian evaluated at points.

Source code in glow2d/_glow2d.py

@staticmethod
def get_jacobian_loc2glob_loc(r: np.ndarray, t: np.ndarray, r0: Numeric = EARTH_RADIUS) -> np.ndarray:
    """Jacobian \\(|J_{r\\rightarrow R}|\\) for local to global coordinate transform, evaluated at points in local coordinate.

    $$
        |J_{r\\rightarrow R}| = \\frac{r}{R^3}\\left(r^2 + R_0^2 + 2 r R_0 \\cos{\\phi}\\right)
    $$

    Args:
        r (np.ndarray): 2-dimensional array of r.
        t (np.ndarray): 2-dimensional array of phi.
        r0 (Numeric): Coordinate transform offset. Defaults to EARTH_RADIUS.

    Raises:
        ValueError: Dimension of inputs must be 2.

    Returns:
        np.ndarray: Jacobian evaluated at points.
    """
    if r.ndim != 2 or t.ndim != 2:
        raise ValueError('Dimension of inputs must be 2')
    gt, gr = glow2d_polar.get_global_coords(t, r, r0=r0)
    jac = (r/(gr**3))*((r**2) + (r0**2) + (2*r*r0*np.cos(t)))
    return jac

get_local_coords(t, r, r0=EARTH_RADIUS, meshgrid=True) staticmethod

Get local coordinates from GEO coordinates.

$$r = \sqrt{\left\{ (R\cos{\theta} - R_0)^2 + R^2\sin{\theta}^2 \right\}}, \\ \phi = \arctan{\frac{R\sin{\theta}}{R\cos{\theta} - R_0}}$$

Parameters:

Name	Type	Description	Default
t	np.ndarray | Numeric	Angles in radians.	required
r	np.ndarray | Numeric	Distance in km.	required
r0	Numeric	Distance to origin. Defaults to geopy.distance.EARTH_RADIUS.	EARTH_RADIUS
meshgrid	bool	Optionally convert 1-D inputs to a meshgrid. Defaults to True.	True

Raises:

Type	Description
ValueError	r and t does not have the same dimensions
TypeError	r and t are not numpy.ndarray.

Returns:

Type	Description
np.ndarray, np.ndarray	(angles, distances) in local coordinates.

Source code in glow2d/_glow2d.py

@staticmethod
def get_local_coords(t: np.ndarray | Numeric, r: np.ndarray | Numeric, r0: Numeric = EARTH_RADIUS, meshgrid: bool = True) -> tuple(np.ndarray, np.ndarray):
    """Get local coordinates from GEO coordinates.

    $$
        r = \\sqrt{\\left\\{ (R\\cos{\\theta} - R_0)^2 + R^2\\sin{\\theta}^2 \\right\\}}, \\\\
        \\phi = \\arctan{\\frac{R\\sin{\\theta}}{R\\cos{\\theta} - R_0}}
    $$

    Args:
        t (np.ndarray | Numeric): Angles in radians.
        r (np.ndarray | Numeric): Distance in km.
        r0 (Numeric, optional): Distance to origin. Defaults to geopy.distance.EARTH_RADIUS.
        meshgrid (bool, optional): Optionally convert 1-D inputs to a meshgrid. Defaults to True.

    Raises:
        ValueError: ``r`` and ``t`` does not have the same dimensions
        TypeError: ``r`` and ``t`` are not ``numpy.ndarray``.

    Returns:
        (np.ndarray, np.ndarray): (angles, distances) in local coordinates.
    """
    if isinstance(r, np.ndarray) and isinstance(t, np.ndarray):
        if r.ndim != t.ndim:  # if dims don't match get out
            raise ValueError('r and t does not have the same dimensions')
        if r.ndim == 1 and meshgrid:
            _r, _t = np.meshgrid(r, t)
        elif r.ndim == 1 and not meshgrid:
            _r, _t = r, t
        else:
            _r, _t = r.copy(), t.copy()  # already a meshgrid?
            r = _r[0]
            t = _t[:, 0]
    elif isinstance(r, Numeric) and isinstance(t, Numeric):
        _r = np.atleast_1d(r)
        _t = np.atleast_1d(t)
    else:
        raise TypeError('r and t must be np.ndarray.')
    rr = np.sqrt((_r*np.cos(_t) - r0)**2 +
                 (_r*np.sin(_t))**2)  # R, T to r, la
    tt = np.arctan2(_r*np.sin(_t), _r*np.cos(_t) - r0)
    return tt, rr

transform_coord()

Run the coordinate transform to convert GLOW output from GEO to local coordinate system.

Returns:

Type	Description
xarray.Dataset	xarray.Dataset: GLOW output in (ZA, r) coordinates. This is a reference and should not be modified.

Source code in glow2d/_glow2d.py

def transform_coord(self) -> xarray.Dataset:
    """Run the coordinate transform to convert GLOW output from GEO to local coordinate system.

    Returns:
        xarray.Dataset: GLOW output in (ZA, r) coordinates. This is a reference and should not be modified.
    """
    if self._iono is not None:
        return self._iono
    tt, rr = self.get_local_coords(
        self._bds.angle.values, self._bds.alt_km.values + EARTH_RADIUS, r0=self._r0)  # get local coords from geocentric coords

    self._rmin, self._rmax = self._bds.alt_km.values.min(), rr.max()  # nearest and farthest local pts
    # highest and lowest za
    self._tmin, self._tmax = tt.min(), np.pi / 2  # 0, tt.max()
    self._nr_num = round(len(self._bds.alt_km.values) * self._resamp)  # resample
    self._nt_num = round(len(self._bds.angle.values) * self._resamp)   # resample
    outp_shape = (self._nt_num, self._nr_num)

    # ttidx = np.where(tt < 0)  # angle below horizon (LA < 0)
    # # get distribution of global -> local points in local grid
    # res = np.histogram2d(rr.flatten(), tt.flatten(), range=([rr.min(), rr.max()], [0, tt.max()]))
    # gd = resize(res[0], outp_shape, mode='edge')  # remap to right size
    # gd *= res[0].sum() / gd.sum()  # conserve sum of points
    # window_length = int(25 * self._resamp)  # smoothing window
    # window_length = window_length if window_length % 2 else window_length + 1  # must be odd
    # gd = savgol_filter(gd, window_length=window_length, polyorder=5, mode='nearest')  # smooth the distribution

    self._altkm = altkm = self._bds.alt_km.values  # store the altkm
    self._theta = theta = self._bds.angle.values  # store the angles
    rmin, rmax = self._rmin, self._rmax  # local names
    tmin, tmax = self._tmin, self._tmax
    self._nr = nr = np.linspace(
        rmin, rmax, self._nr_num, endpoint=True)  # local r
    self._nt = nt = np.linspace(
        tmin, tmax, self._nt_num, endpoint=True)  # local look angle
    # get meshgrid of the R, T coord system from regular r, za grid
    self._ntt, self._nrr = self.get_global_coords(nt, nr, r0=self._r0)
    # calculate jacobian
    jacobian = self.get_jacobian_glob2loc_glob(self._nrr, self._ntt, r0=self._r0)
    # convert to pixel coordinates
    self._ntt = self._ntt.flatten()  # flatten T, works as _global_from_local LUT
    self._nrr = self._nrr.flatten()  # flatten R, works as _global_from_local LUT
    self._ntt = (self._ntt - self._theta.min()) / \
        (self._theta.max() - self._theta.min()) * \
        len(self._theta)  # calculate equivalent index (pixel coord) from original T grid
    self._nrr = (self._nrr - self._altkm.min() - self._r0) / \
        (self._altkm.max() - self._altkm.min()) * \
        len(self._altkm)  # calculate equivalent index (pixel coord) from original R (alt_km) grid
    # start transformation
    data_vars = {}
    bds = self._bds
    coord_wavelength = bds.wavelength.values  # wl axis
    coord_state = bds.state.values  # state axis
    coord_energy = bds.energy.values  # egrid
    bds_attr = bds.attrs  # attrs
    single_keys = ['Tn',
                   'pedersen',
                   'hall',
                   'Te',
                   'Ti']  # (angle, alt_km) vars
    density_keys = [
        'O',
        'N2',
        'O2',
        'NO',
        'NeIn',
        'NeOut',
        'ionrate',
        'O+',
        'O2+',
        'NO+',
        'N2D',
    ]
    state_keys = [
        'production',
        'loss',
        'excitedDensity'
    ]  # (angle, alt_km, state) vars
    # start = perf_counter_ns()
    # map all the single key types from (angle, alt_km) -> (la, r)
    for key in single_keys:
        inp = self._bds[key].values
        inp[np.where(np.isnan(inp))] = 0
        out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape, mode='nearest')
        out[np.where(out < 0)] = 0
        # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape)
        data_vars[key] = (('za', 'r'), out)
    for key in density_keys:
        inp = self._bds[key].values
        inp[np.where(np.isnan(inp))] = 0
        out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape, mode='nearest') / jacobian
        out[np.where(out < 0)] = 0
        # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape)
        data_vars[key] = (('za', 'r'), out)
    # end = perf_counter_ns()
    # print('Single_key conversion: %.3f us'%((end - start)*1e-3))
    # start = perf_counter_ns()
    # dataset of (angle, alt_km) vars
    iono = xarray.Dataset(data_vars=data_vars, coords={
        'za': nt, 'r': nr})
    # end = perf_counter_ns()
    # print('Single_key dataset: %.3f us'%((end - start)*1e-3))
    # start = perf_counter_ns()
    ver = []
    # map all the wavelength data from (angle, alt_km, wavelength) -> (la, r, wavelength)
    for key in coord_wavelength:
        inp = bds['ver'].loc[dict(wavelength=key)].values
        inp[np.where(np.isnan(inp))] = 0
        # scaled by point distribution because flux is conserved, not brightness
        # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape, mode='nearest') * gd
        out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape, mode='nearest') / jacobian
        out[np.where(out < 0)] = 0
        # inp[ttidx] = 0
        # inpsum = inp.sum()  # sum of input for valid angles
        # outpsum = out.sum()  # sum of output
        # out = out * (inpsum / outpsum)  # scale the sum to conserve total flux
        ver.append(out.T)
    # end = perf_counter_ns()
    # print('VER eval: %.3f us'%((end - start)*1e-3))
    # start = perf_counter_ns()
    ver = np.asarray(ver).T
    ver = xr.DataArray(
        ver,
        coords={'za': nt, 'r': nr,
                'wavelength': coord_wavelength},
        dims=['za', 'r', 'wavelength'],
        name='ver'
    )  # create wl dataset
    # end = perf_counter_ns()
    # print('VER dataset: %.3f us'%((end - start)*1e-3))
    # start = perf_counter_ns()
    if self._wprodloss:
        d = {}
        for key in state_keys:  # for each var with (angle, alt_km, state)
            res = []

            def convert_state_stuff(st):
                inp = bds[key].loc[dict(state=st)].values
                inp[np.where(np.isnan(inp))] = 0
                out = geometric_transform(inp, mapping=self._global_from_local, output_shape=outp_shape)
                out[np.where(out < 0)] = 0
                if key in ('production', 'excitedDensity'):
                    out /= jacobian
                # out = warp(inp, inverse_map=(2, self._ntt, self._nrr), output_shape=outp_shape)
                return out.T
            res = list(map(convert_state_stuff, coord_state))
            res = np.asarray(res).T
            d[key] = (('za', 'r', 'state'), res)
        # end = perf_counter_ns()
        # print('Prod_Loss Eval: %.3f us'%((end - start)*1e-3))
        # start = perf_counter_ns()
        prodloss = xarray.Dataset(
            data_vars=d,
            coords={'za': nt, 'r': nr, 'state': coord_state}
        )  # calculate (angle, alt_km, state) -> (la, r, state) dataset
    else:
        prodloss = xarray.Dataset()
    # end = perf_counter_ns()
    # print('Prod_Loss DS: %.3f us'%((end - start)*1e-3))
    ## EGrid conversion (angle, energy) -> (r, energy) ##
    # EGrid is avaliable really at (angle, alt_km = 0, energy)
    # So we get local coords for (angle, R=R0)
    # we discard the angle information because it is meaningless, EGrid is spatial
    # start = perf_counter_ns()
    _rr, _ = self.get_local_coords(
        bds.angle.values, np.ones(bds.angle.values.shape)*self._r0, r0=self._r0)
    _rr = rr[:, 0]  # spatial EGrid
    d = []
    for en in coord_energy:  # for each energy
        inp = bds['precip'].loc[dict(energy=en)].values
        # interpolate to appropriate energy grid
        out = np.interp(nr, _rr, inp)
        d.append(out)
    d = np.asarray(d).T
    precip = xarray.Dataset({'precip': (('r', 'energy'), d)}, coords={
        'r': nr, 'energy': coord_energy})
    # end = perf_counter_ns()
    # print('Precip interp and ds: %.3f us'%((end - start)*1e-3))

    # start = perf_counter_ns()
    iono = xr.merge((iono, ver, prodloss, precip))  # merge all datasets
    bds_attr['altitude'] = {'values': self._r0 - EARTH_RADIUS, 'units': 'km',
                            'description': 'Altitude of local polar coordinate origin ASL'}
    iono.attrs.update(bds_attr)  # copy original attrs

    _ = list(map(lambda x: iono[x].attrs.update(bds[x].attrs), tuple(iono.data_vars.keys())))  # update attrs from bds

    unit_desc_dict = {
        'za': ('radians', 'Zenith angle'),
        'r': ('km', 'Radial distance in km')
    }
    _ = list(map(lambda x: iono[x].attrs.update(
        {'units': unit_desc_dict[x][0], 'description': unit_desc_dict[x][1]}), unit_desc_dict.keys()))
    # end = perf_counter_ns()
    # print('Merging: %.3f us'%((end - start)*1e-3))
    self._iono = iono
    return iono