Module pyaurorax.tools.bounding_box.extract_metric

Extract various metrics from a given bounding box.

# Copyright 2024 University of Calgary
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
"""
Extract various metrics from a given bounding box.
"""

from ._geo import geo
from ._ccd import ccd
from ._mag import mag
from ._elevation import elevation
from ._azimuth import azimuth

__all__ = [
    "geo",
    "ccd",
    "mag",
    "elevation",
    "azimuth",
]

Functions

def azimuth(images: numpy.ndarray, skymap: pyucalgarysrs.data.classes.Skymap, azimuth_bounds: Sequence[Union[float, int]], metric: Literal['mean', 'median', 'sum'] = 'median', n_channels: Optional[int] = None, show_preview: bool = False) ‑> numpy.ndarray

Compute a metric of image data within an azimuthal boundary.

Args

images : numpy.ndarray

A set of images. Normally this would come directly from a data read call, but can also be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images] or [row, cols, channels, num_images].

skymap : Skymap

The skymap corresponding to the image data.

azimuth_bounds (Sequence[int | float]: A 2-element sequence specifying the azimuthal bounds from which to extract the metric. Anticipated order is [az_min, az_max].

metric : str

The name of the metric that is to be computed for the bounded area. Valid metrics are mean, median, sum. Default is median.

n_channels : int

By default, function will assume the type of data passed as input - this argument can be used to manually specify the number of channels contained in image data.

show_preview (bool): Plot a preview of the bounded area.

Returns

A numpy.ndarray containing the metrics computed within azimuth range, for all image frames.

Raises

ValueError

issue encountered with value supplied in parameter

Expand source code

def azimuth(images: np.ndarray,
            skymap: Skymap,
            azimuth_bounds: Sequence[Union[int, float]],
            metric: Literal["mean", "median", "sum"] = "median",
            n_channels: Optional[int] = None,
            show_preview: bool = False) -> np.ndarray:
    """
    Compute a metric of image data within an azimuthal boundary.

    Args:
        images (numpy.ndarray): 
            A set of images. Normally this would come directly from a data `read` call, but can also
            be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images]
            or [row, cols, channels, num_images].
        
        skymap (pyaurorax.data.ucalgary.Skymap): 
            The skymap corresponding to the image data.

        azimuth_bounds (Sequence[int | float]: 
            A 2-element sequence specifying the azimuthal bounds from which to extract the metric. 
            Anticipated order is [az_min, az_max].

        metric (str): 
            The name of the metric that is to be computed for the bounded area. Valid metrics are `mean`,
            `median`, `sum`. Default is `median`.

        n_channels (int): 
            By default, function will assume the type of data passed as input - this argument can be used
            to manually specify the number of channels contained in image data.

        show_preview (bool):
            Plot a preview of the bounded area.

    Returns:
        A numpy.ndarray containing the metrics computed within azimuth range, for all image frames.

    Raises:
        ValueError: issue encountered with value supplied in parameter
    """

    # Select individual azimuths from list
    az_0 = azimuth_bounds[0]
    az_1 = azimuth_bounds[1]

    # Determine if we are using single or 3 channel
    images = np.squeeze(images)
    if n_channels is not None:
        n_channels = n_channels
    else:
        n_channels = 1
        if (len(images.shape) == 3):
            # single channel
            n_channels = 1
        elif (len(images.shape) == 4):
            # three channel
            n_channels = 3

    # Ensure that coordinates are valid
    if az_0 > 360 or az_0 < 0:
        raise ValueError("Invalid Azimuth: " + str(az_0))
    elif az_1 > 360 or az_1 < 0:
        raise ValueError("Invalid Azimuth: " + str(az_1))

    # Ensure that azimuths are properly ordered
    if az_0 > az_1:
        az_0, az_1 = az_1, az_0

    # Ensure that this is a valid polygon
    if (az_0 == az_1):
        raise ValueError("Azimuth bounds defined with zero area.")

    # Obtain azimuth array from skymap
    az = np.squeeze(skymap.full_azimuth)

    # Obtain indices into skymap within azimuth range
    bound_idx = np.where(np.logical_and(az > float(az_0), az < float(az_1)))

    # If boundaries contain no data, raise error
    if len(bound_idx[0]) == 0 or len(bound_idx[1]) == 0:
        raise ValueError("No data within desired bounds. Try a larger area.")

    # Slice out the bounded data, and plot preview if requested
    if n_channels == 1:
        bound_data = images[bound_idx[0], bound_idx[1], :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1]] = 255
            plt.figure()
            plt.imshow(preview_img, cmap="grey", origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    elif n_channels == 3:
        bound_data = images[bound_idx[0], bound_idx[1], :, :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1], 0] = 255
            preview_img[bound_idx[0], bound_idx[1], 1:] = 0
            plt.figure()
            plt.imshow(preview_img, origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    else:
        raise ValueError("Unrecognized image format with shape: " + str(images.shape))

    # Compute metric of interest
    if metric == 'median':
        result = np.median(bound_data, axis=0)
    elif metric == 'mean':
        result = np.mean(bound_data, axis=0)
    elif metric == 'sum':
        result = np.sum(bound_data, axis=0)
    else:
        raise ValueError("Metric " + str(metric) + " is not recognized.")

    return result

def ccd(images: numpy.ndarray, ccd_bounds: Sequence[int], metric: Literal['mean', 'median', 'sum'] = 'median', n_channels: Optional[int] = None, show_preview: bool = False) ‑> numpy.ndarray

Compute a metric of image data within a CCD boundary.

Args

images : numpy.ndarray

A set of images. Normally this would come directly from a data read call, but can also be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images] or [row, cols, channels, num_images].

ccd_bounds : List[int]

A 4-element sequence specifying the (inclusive) CCD bounds from which to extract the metric. Anticipated order is [x_0, x_1, y_0, y_1].

metric : str

The name of the metric that is to be computed for the bounded area. Valid metrics are mean, median, sum. Defaults to median.

n_channels : int

By default, function will assume the type of data passed as input - this argument can be used to manually specify the number of channels contained in image data.

show_preview (bool): Plot a preview of the bounded area.

Returns

A numpy.ndarray containing the metrics computed within CCD bounds, for all image frames.

Raises

ValueError

issue encountered with value supplied in parameter

Expand source code

def ccd(images: np.ndarray,
        ccd_bounds: Sequence[int],
        metric: Literal["mean", "median", "sum"] = "median",
        n_channels: Optional[int] = None,
        show_preview: bool = False) -> np.ndarray:
    """
    Compute a metric of image data within a CCD boundary.

    Args:
        images (numpy.ndarray): 
            A set of images. Normally this would come directly from a data `read` call, but can also
            be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images]
            or [row, cols, channels, num_images].
        
        ccd_bounds (List[int]): 
            A 4-element sequence specifying the (inclusive) CCD bounds from which to extract the metric. 
            Anticipated order is [x_0, x_1, y_0, y_1].

        metric (str): 
            The name of the metric that is to be computed for the bounded area. Valid metrics are `mean`,
            `median`, `sum`. Defaults to `median`.

        n_channels (int): 
            By default, function will assume the type of data passed as input - this argument can be used
            to manually specify the number of channels contained in image data.
        
        show_preview (bool):
            Plot a preview of the bounded area.

    Returns:
        A numpy.ndarray containing the metrics computed within CCD bounds, for all image frames.

    Raises:
        ValueError: issue encountered with value supplied in parameter
    """

    # Select individual ccd x/y from list
    x_0 = ccd_bounds[0]
    x_1 = ccd_bounds[1]
    y_0 = ccd_bounds[2]
    y_1 = ccd_bounds[3]

    # Determine if we are using single or 3 channel
    images = np.squeeze(images)
    if n_channels is not None:
        n_channels = n_channels
    else:
        n_channels = 1
        if (len(images.shape) == 3):
            # single channel
            n_channels = 1
        elif (len(images.shape) == 4):
            # three channel
            n_channels = 3

    # Ensure that coordinates are valid
    max_x = images.shape[0]
    max_y = images.shape[1]

    if y_0 > max_y or y_0 < 0:
        raise ValueError("CCD Y Coordinate " + str(y_0) + " out of range for image of shape " + str((max_y, max_x)) + ".")
    elif y_1 > max_y or y_1 < 0:
        raise ValueError("CCD Y Coordinate " + str(y_1) + " out of range for image of shape " + str((max_y, max_x)) + ".")
    elif x_0 > max_x or x_0 < 0:
        raise ValueError("CCD Y Coordinate " + str(x_0) + " out of range for image of shape " + str((max_y, max_x)) + ".")
    elif x_1 > max_x or x_1 < 0:
        raise ValueError("CCD Y Coordinate " + str(x_1) + " out of range for image of shape " + str((max_y, max_x)) + ".")

    # Ensure that coordinates are properly ordered
    if y_0 > y_1:
        y_0, y_1 = y_1, y_0
    if x_0 > x_1:
        x_0, x_1 = x_1, x_0

    # Ensure that this is a valid polygon
    if (y_0 == y_1) or (x_0 == x_1):
        raise ValueError("Polygon defined with zero area.")

    # Slice out the bounded data
    if n_channels == 1:
        bound_data = images[y_0:y_1, x_0:x_1, :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, 0], top=230)
            preview_img[y_0:y_1, x_0:x_1] = 255
            plt.figure()
            plt.imshow(preview_img, cmap="grey", origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    elif n_channels == 3:
        bound_data = images[y_0:y_1, x_0:x_1, :, :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, :, 0], top=230)
            preview_img[y_0:y_1, x_0:x_1, 0] = 255
            preview_img[y_0:y_1, x_0:x_1, 1:] = 0
            plt.figure()
            plt.imshow(preview_img, origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    else:
        raise ValueError("Unrecognized image format with shape: " + str(images.shape))

    # Compute metric of interest
    if metric == 'median':
        result = np.median(bound_data, axis=(0, 1))
    elif metric == 'mean':
        result = np.mean(bound_data, axis=(0, 1))
    elif metric == 'sum':
        result = np.sum(bound_data, axis=(0, 1))
    else:
        raise ValueError("Metric " + str(metric) + " is not recognized.")

    return result

def elevation(images: numpy.ndarray, skymap: pyucalgarysrs.data.classes.Skymap, elevation_bounds: Sequence[Union[float, int]], metric: Literal['mean', 'median', 'sum'] = 'median', n_channels: Optional[int] = None, show_preview: bool = False) ‑> numpy.ndarray

Compute a metric of image data within an elevation boundary.

Args

images : numpy.ndarray

A set of images. Normally this would come directly from a data read call, but can also be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images] or [row, cols, channels, num_images].

skymap : Skymap

The skymap corresponding to the image data.

elevation_bounds : Sequence

A 2-element sequence specifying the elevation bounds from which to extract the metric. Anticipated order is [el_min, el_max].

metric : str

The name of the metric that is to be computed for the bounded area. Valid metrics are mean, median, sum. Default is median.

n_channels : int

By default, function will assume the type of data passed as input - this argument can be used to manually specify the number of channels contained in image data.

show_preview (bool): Plot a preview of the bounded area.

Returns

A numpy.ndarray containing the metrics computed within elevation range, for all image frames.

Raises

ValueError

issue encountered with value supplied in parameter

Expand source code

def elevation(images: np.ndarray,
              skymap: Skymap,
              elevation_bounds: Sequence[Union[int, float]],
              metric: Literal["mean", "median", "sum"] = "median",
              n_channels: Optional[int] = None,
              show_preview: bool = False) -> np.ndarray:
    """
    Compute a metric of image data within an elevation boundary.

    Args:
        images (numpy.ndarray): 
            A set of images. Normally this would come directly from a data `read` call, but can also
            be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images]
            or [row, cols, channels, num_images].
        
        skymap (pyaurorax.data.ucalgary.Skymap): 
            The skymap corresponding to the image data.

        elevation_bounds (Sequence): 
            A 2-element sequence specifying the elevation bounds from which to extract the metric. 
            Anticipated order is [el_min, el_max].

        metric (str): 
            The name of the metric that is to be computed for the bounded area. Valid metrics are `mean`,
            `median`, `sum`. Default is `median`.

        n_channels (int): 
            By default, function will assume the type of data passed as input - this argument can be used
            to manually specify the number of channels contained in image data.
        
        show_preview (bool):
            Plot a preview of the bounded area.


    Returns:
        A numpy.ndarray containing the metrics computed within elevation range, for all image frames.

    Raises:
        ValueError: issue encountered with value supplied in parameter
    """

    # Select individual elevations from list
    elev_0 = elevation_bounds[0]
    elev_1 = elevation_bounds[1]

    # Determine if we are using single or 3 channel
    images = np.squeeze(images)
    if n_channels is not None:
        n_channels = n_channels
    else:
        n_channels = 1
        if (len(images.shape) == 3):
            # single channel
            n_channels = 1
        elif (len(images.shape) == 4):
            # three channel
            n_channels = 3

    # Ensure that coordinates are valid
    if elev_0 > 90 or elev_0 < 0:
        raise ValueError("Invalid Elevation: " + str(elev_0))
    elif elev_1 > 90 or elev_1 < 0:
        raise ValueError("Invalid Elevation: " + str(elev_1))

    # Ensure that elevations are properly ordered
    if elev_0 > elev_1:
        elev_0, elev_1 = elev_1, elev_0

    # Ensure that this is a valid bounded area
    if (elev_0 == elev_1):
        raise ValueError("Elevation bounds defined with zero area.")

    # Obtain elevation array from skymap
    elev = np.squeeze(skymap.full_elevation)

    # Obtain indices into skymap within elevation range
    bound_idx = np.where(np.logical_and(elev >= float(elev_0), elev <= float(elev_1)))

    # If boundaries contain no data, raise error
    if len(bound_idx[0]) == 0 or len(bound_idx[1]) == 0:
        raise ValueError("No data within desired bounds. Try a larger area.")

    # Slice out the bounded data
    if n_channels == 1:
        bound_data = images[bound_idx[0], bound_idx[1], :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1]] = 255
            plt.figure()
            plt.imshow(preview_img, cmap="grey", origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    elif n_channels == 3:
        bound_data = images[bound_idx[0], bound_idx[1], :, :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1], 0] = 255
            preview_img[bound_idx[0], bound_idx[1], 1:] = 0
            plt.figure()
            plt.imshow(preview_img, origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    else:
        raise ValueError("Unrecognized image format with shape: " + str(images.shape))

    # Compute metric of interest
    if metric == 'median':
        result = np.median(bound_data, axis=0)
    elif metric == 'mean':
        result = np.mean(bound_data, axis=0)
    elif metric == 'sum':
        result = np.sum(bound_data, axis=0)
    else:
        raise ValueError("Metric " + str(metric) + " is not recognized.")

    return result

def geo(images: numpy.ndarray, skymap: pyucalgarysrs.data.classes.Skymap, altitude_km: Union[float, int], lonlat_bounds: Sequence[Union[float, int]], metric: Literal['mean', 'median', 'sum'] = 'median', n_channels: Optional[int] = None, show_preview: bool = False) ‑> numpy.ndarray

Compute a metric of image data within a geographic lat/lon boundary.

Args

images : numpy.ndarray

A set of images. Normally this would come directly from a data read call, but can also be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images] or [row, cols, channels, num_images].

skymap : Skymap

The skymap corresponding to the image data.

altitude_km : int or float

The altitude of the image data in kilometers.

lonlat_bounds : Sequence

A 4-element sequence specifying the lat/lon bounds from which to extract the metric. Anticipated order is [lon_0, lon_1, lat_0, lat_1].

metric : str

The name of the metric that is to be computed for the bounded area. Valid metrics are mean, median, sum. Default is median.

n_channels : int

By default, function will assume the type of data passed as input - this argument can be used to manually specify the number of channels contained in image data.

show_preview (bool): Plot a preview of the bounded area.

Returns

A numpy.ndarray containing the metrics computed within elevation range, for all image frames.

Raises

ValueError

issue encountered with value supplied in parameter

Expand source code

def geo(images: np.ndarray,
        skymap: Skymap,
        altitude_km: Union[int, float],
        lonlat_bounds: Sequence[Union[int, float]],
        metric: Literal["mean", "median", "sum"] = "median",
        n_channels: Optional[int] = None,
        show_preview: bool = False) -> np.ndarray:
    """
    Compute a metric of image data within a geographic lat/lon boundary.

    Args:
        images (numpy.ndarray): 
            A set of images. Normally this would come directly from a data `read` call, but can also
            be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images]
            or [row, cols, channels, num_images].
        
        skymap (pyaurorax.data.ucalgary.Skymap): 
            The skymap corresponding to the image data.
        
        altitude_km (int or float): 
            The altitude of the image data in kilometers.

        lonlat_bounds (Sequence): 
            A 4-element sequence specifying the lat/lon bounds from which to extract the metric. 
            Anticipated order is [lon_0, lon_1, lat_0, lat_1].

        metric (str): 
            The name of the metric that is to be computed for the bounded area. Valid metrics are `mean`,
            `median`, `sum`. Default is `median`.

        n_channels (int): 
            By default, function will assume the type of data passed as input - this argument can be used
            to manually specify the number of channels contained in image data.
        
        show_preview (bool):
            Plot a preview of the bounded area.

    Returns:
        A numpy.ndarray containing the metrics computed within elevation range, for all image frames.

    Raises:
        ValueError: issue encountered with value supplied in parameter
    """

    # Select individual lats/lons from list
    lon_0 = lonlat_bounds[0]
    lon_1 = lonlat_bounds[1]
    lat_0 = lonlat_bounds[2]
    lat_1 = lonlat_bounds[3]

    # Determine if we are using single or 3 channel
    images = np.squeeze(images)
    if n_channels is not None:
        n_channels = n_channels
    else:
        n_channels = 1
        if (len(images.shape) == 3):
            # single channel
            n_channels = 1
        elif (len(images.shape) == 4):
            # three channel
            n_channels = 3

    # Ensure that coordinates are valid
    if lat_0 > 90 or lat_0 < -90:
        raise ValueError("Invalid Latitude: " + str(lat_0))
    elif lat_1 > 90 or lat_1 < -90:
        raise ValueError("Invalid Latitude: " + str(lat_1))
    elif lon_0 > 360 or lon_0 < -180:
        raise ValueError("Invalid Longitude: " + str(lon_0))
    elif lon_1 > 360 or lon_1 < -180:
        raise ValueError("Invalid Longitude: " + str(lon_1))

    # Convert(0,360) longitudes to (-180,180) if entered as such
    if lon_0 > 180:
        lon_0 -= 360.0
    if lon_1 > 180:
        lon_1 -= 360.0

    # Ensure that coordinates are properly ordered
    if lat_0 > lat_1:
        lat_0, lat_1 = lat_1, lat_0
    if lon_0 > lon_1:
        lon_0, lon_1 = lon_1, lon_0

    # Ensure that this is a valid polygon
    if (lat_0 == lat_1) or (lon_0 == lon_1):
        raise ValueError("Polygon defined with zero area.")

    # Obtain lat/lon arrays from skymap
    interp_alts = skymap.full_map_altitude / 1000.0
    if (altitude_km in interp_alts):
        altitude_idx = np.where(altitude_km == interp_alts)

        lats = np.squeeze(skymap.full_map_latitude[altitude_idx, :, :])
        lons = np.squeeze(skymap.full_map_longitude[altitude_idx, :, :])
        lons[np.where(lons > 180)] -= 360.0  # Fix skymap to be in (-180,180) format

    else:
        # Make sure altitude is in range that can be interpolated
        if (altitude_km < interp_alts[0]) or (altitude_km > interp_alts[2]):
            raise ValueError("Altitude " + str(altitude_km) + " outside valid range of " + str((interp_alts[0], interp_alts[2])))

        # Initialze empty lat/lon arrays
        lats = np.full(np.squeeze(skymap.full_map_latitude[0, :, :]).shape, np.nan, dtype=skymap.full_map_latitude[0, :, :].dtype)
        lons = np.full(np.squeeze(skymap.full_map_latitude[0, :, :]).shape, np.nan, dtype=skymap.full_map_latitude[0, :, :].dtype)

        # Interpolate lats and lons at desired altitude
        for i in range(skymap.full_map_latitude.shape[1]):
            for j in range(skymap.full_map_latitude.shape[2]):
                pixel_lats = skymap.full_map_latitude[:, i, j]
                pixel_lons = skymap.full_map_longitude[:, i, j]
                if np.isnan(pixel_lats).any() or np.isnan(pixel_lons).any():
                    continue
                lats[i, j] = np.interp(altitude_km, interp_alts, pixel_lats)
                lons[i, j] = np.interp(altitude_km, interp_alts, pixel_lons)

        lons[np.where(lons > 180)] -= 360.0  # Fix skymap to be in (-180,180) format

    # Check that lat/lon range is reasonable
    min_skymap_lat = np.nanmin(lats)
    max_skymap_lat = np.nanmax(lats)
    min_skymap_lon = np.nanmin(lons)
    max_skymap_lon = np.nanmax(lons)
    if (lat_0 <= min_skymap_lat) or (lat_1 >= max_skymap_lat):
        raise ValueError(f"Latitude range supplied is outside the valid range for this skymap {(min_skymap_lat,max_skymap_lat)}.")
    if (lon_0 <= min_skymap_lon) or (lon_1 >= max_skymap_lon):
        raise ValueError(f"Longitude range supplied is outside the valid range for this skymap {(min_skymap_lon,max_skymap_lon)}.")

    # Obtain indices into skymap within lat/lon range
    bound_idx = np.where(np.logical_and.reduce((lats >= float(lat_0), lats <= float(lat_1), lons >= float(lon_0), lons <= float(lon_1))))

    # If boundaries contain no data, raise error
    if len(bound_idx[0]) == 0 or len(bound_idx[1]) == 0:
        raise ValueError("No data within desired bounds. Try a larger area.")

    # Convert from skymap coords to image coords
    bound_idx = tuple(i - 1 for i in bound_idx)
    bound_idx = tuple(np.maximum(idx, 0) for idx in bound_idx)

    # Slice out the bounded data
    if n_channels == 1:
        bound_data = images[bound_idx[0], bound_idx[1], :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1]] = 255
            plt.figure()
            plt.imshow(preview_img, cmap="grey", origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    elif n_channels == 3:
        bound_data = images[bound_idx[0], bound_idx[1], :, :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1], 0] = 255
            preview_img[bound_idx[0], bound_idx[1], 1:] = 0
            plt.figure()
            plt.imshow(preview_img, origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    else:
        raise ValueError("Unrecognized image format with shape: " + str(images.shape))

    # Compute metric of interest
    if metric == 'median':
        result = np.median(bound_data, axis=0)
    elif metric == 'mean':
        result = np.mean(bound_data, axis=0)
    elif metric == 'sum':
        result = np.sum(bound_data, axis=0)
    else:
        raise ValueError("Metric " + str(metric) + " is not recognized.")

    return result

def mag(images: numpy.ndarray, timestamp: datetime.datetime, skymap: pyucalgarysrs.data.classes.Skymap, altitude_km: Union[float, int], lonlat_bounds: Sequence[Union[float, int]], metric: Literal['mean', 'median', 'sum'] = 'median', n_channels: Optional[int] = None, show_preview: bool = False) ‑> numpy.ndarray

Compute a metric of image data within a magnetic lat/lon boundary.

Args

images : numpy.ndarray

A set of images. Normally this would come directly from a data read call, but can also be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images] or [row, cols, channels, num_images].

skymap : Skymap

The skymap corresponding to the image data.

altitude_km : int or float

The altitude of the image data in kilometers.

lonlat_bounds : Sequence

A 4-element sequence specifying the magnetic lat/lon bounds from which to extract the metric. Anticipated order is [lon_0, lon_1, lat_0, lat_1].

metric : str

The name of the metric that is to be computed for the bounded area. Valid metrics are mean, median, sum. Default is median.

n_channels : int

By default, function will assume the type of data passed as input - this argument can be used to manually specify the number of channels contained in image data.

show_preview (bool): Plot a preview of the bounded area.

Returns

A numpy.ndarray containing the metrics computed within elevation range, for all image frames.

Raises

ValueError

issue encountered with value supplied in parameter

Expand source code

def mag(images: np.ndarray,
        timestamp: datetime.datetime,
        skymap: Skymap,
        altitude_km: Union[int, float],
        lonlat_bounds: Sequence[Union[int, float]],
        metric: Literal["mean", "median", "sum"] = "median",
        n_channels: Optional[int] = None,
        show_preview: bool = False) -> np.ndarray:
    """
    Compute a metric of image data within a magnetic lat/lon boundary.

    Args:
        images (numpy.ndarray): 
            A set of images. Normally this would come directly from a data `read` call, but can also
            be any arbitrary set of images. It is anticipated that the order of axes is [rows, cols, num_images]
            or [row, cols, channels, num_images].
        
        skymap (pyaurorax.data.ucalgary.Skymap): 
            The skymap corresponding to the image data.
        
        altitude_km (int or float): 
            The altitude of the image data in kilometers.

        lonlat_bounds (Sequence): 
            A 4-element sequence specifying the magnetic lat/lon bounds from which to extract the metric. 
            Anticipated order is [lon_0, lon_1, lat_0, lat_1].

        metric (str): 
            The name of the metric that is to be computed for the bounded area. Valid metrics are `mean`,
            `median`, `sum`. Default is `median`.

        n_channels (int): 
            By default, function will assume the type of data passed as input - this argument can be used
            to manually specify the number of channels contained in image data.

        show_preview (bool):
            Plot a preview of the bounded area.

    Returns:
        A numpy.ndarray containing the metrics computed within elevation range, for all image frames.

    Raises:
        ValueError: issue encountered with value supplied in parameter
    """

    # Select individual lats/lons from list
    lon_0 = lonlat_bounds[0]
    lon_1 = lonlat_bounds[1]
    lat_0 = lonlat_bounds[2]
    lat_1 = lonlat_bounds[3]

    # Determine if we are using single or 3 channel
    images = np.squeeze(images)
    if n_channels is not None:
        n_channels = n_channels
    else:
        n_channels = 1
        if (len(images.shape) == 3):
            # single channel
            n_channels = 1
        elif (len(images.shape) == 4):
            # three channel
            n_channels = 3

    # Ensure that coordinates are valid
    if lat_0 > 90 or lat_0 < -90:
        raise ValueError("Invalid Latitude: " + str(lat_0))
    elif lat_1 > 90 or lat_1 < -90:
        raise ValueError("Invalid Latitude: " + str(lat_1))
    elif lon_0 > 360 or lon_0 < -180:
        raise ValueError("Invalid Longitude: " + str(lon_0))
    elif lon_1 > 360 or lon_1 < -180:
        raise ValueError("Invalid Longitude: " + str(lon_1))

    # Convert(0,360) longitudes to (-180,180) if entered as such
    if lon_0 > 180:
        lon_0 -= 360.0
    if lon_1 > 180:
        lon_1 -= 360.0

    # Ensure that coordinates are properly ordered
    if lat_0 > lat_1:
        lat_0, lat_1 = lat_1, lat_0
    if lon_0 > lon_1:
        lon_0, lon_1 = lon_1, lon_0

    # Ensure that this is a valid polygon
    if (lat_0 == lat_1) or (lon_0 == lon_1):
        raise ValueError("Polygon defined with zero area.")

    # Obtain lat/lon arrays from skymap
    interp_alts = skymap.full_map_altitude / 1000.0
    if (altitude_km in interp_alts):
        altitude_idx = np.where(altitude_km == interp_alts)

        lats = np.squeeze(skymap.full_map_latitude[altitude_idx, :, :])
        lons = np.squeeze(skymap.full_map_longitude[altitude_idx, :, :])
        lons[np.where(lons > 180)] -= 360.0  # Fix skymap to be in (-180,180) format

    else:
        # Make sure altitude is in range that can be interpolated
        if (altitude_km < interp_alts[0]) or (altitude_km > interp_alts[2]):
            raise ValueError("Altitude " + str(altitude_km) + " outside valid range of " + str((interp_alts[0], interp_alts[2])))

        # Initialze empty lat/lon arrays
        lats = np.full(np.squeeze(skymap.full_map_latitude[0, :, :]).shape, np.nan, dtype=skymap.full_map_latitude[0, :, :].dtype)
        lons = np.full(np.squeeze(skymap.full_map_latitude[0, :, :]).shape, np.nan, dtype=skymap.full_map_latitude[0, :, :].dtype)

        # Interpolate lats and lons at desired altitude
        for i in range(skymap.full_map_latitude.shape[1]):
            for j in range(skymap.full_map_latitude.shape[2]):
                pixel_lats = skymap.full_map_latitude[:, i, j]
                pixel_lons = skymap.full_map_longitude[:, i, j]
                if np.isnan(pixel_lats).any() or np.isnan(pixel_lons).any():
                    continue
                lats[i, j] = np.interp(altitude_km, interp_alts, pixel_lats)
                lons[i, j] = np.interp(altitude_km, interp_alts, pixel_lons)

        lons[np.where(lons > 180)] -= 360.0  # Fix skymap to be in (-180,180) format

    # Convert skymap to magnetic coords
    mag_lats, mag_lons, mag_alts = aacgmv2.convert_latlon_arr(lats.flatten(), lons.flatten(), (lons * 0.0).flatten(), timestamp, method_code='G2A')
    mag_lats = np.reshape(mag_lats, lats.shape)
    mag_lons = np.reshape(mag_lons, lons.shape)

    # Check that lat/lon range is reasonable
    min_skymap_lat = np.nanmin(mag_lats)
    max_skymap_lat = np.nanmax(mag_lats)
    min_skymap_lon = np.nanmin(mag_lons)
    max_skymap_lon = np.nanmax(mag_lons)
    if (lat_0 <= min_skymap_lat) or (lat_1 >= max_skymap_lat):
        raise ValueError(f"Latitude range supplied is outside the valid range for this skymap {(min_skymap_lat,max_skymap_lat)}.")
    if (lon_0 <= min_skymap_lon) or (lon_1 >= max_skymap_lon):
        raise ValueError(f"Latitude range supplied is outside the valid range for this skymap {(min_skymap_lon,max_skymap_lon)}.")

    # Obtain indices into skymap within lat/lon range
    bound_idx = np.where(
        np.logical_and.reduce((
            mag_lats >= float(lat_0),
            mag_lats <= float(lat_1),
            mag_lons >= float(lon_0),
            mag_lons <= float(lon_1),
        )))

    # If boundaries contain no data, raise error
    if len(bound_idx[0]) == 0 or len(bound_idx[1]) == 0:
        raise ValueError("No data within desired bounds. Try a larger area.")

    # Convert from skymap coords to image coords
    bound_idx = tuple(i - 1 for i in bound_idx)
    bound_idx = tuple(np.maximum(idx, 0) for idx in bound_idx)

    # Slice out the bounded data
    if n_channels == 1:
        bound_data = images[bound_idx[0], bound_idx[1], :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1]] = 255
            plt.figure()
            plt.imshow(preview_img, cmap="grey", origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    elif n_channels == 3:
        bound_data = images[bound_idx[0], bound_idx[1], :, :]
        if show_preview:
            preview_img = scale_intensity(images[:, :, :, 0], top=230)
            preview_img[bound_idx[0], bound_idx[1], 0] = 255
            preview_img[bound_idx[0], bound_idx[1], 1:] = 0
            plt.figure()
            plt.imshow(preview_img, origin="lower")
            plt.title("Bounded Area Preview")
            plt.axis("off")
            plt.show()
    else:
        raise ValueError("Unrecognized image format with shape: " + str(images.shape))

    # Compute metric of interest
    if metric == 'median':
        result = np.median(bound_data, axis=0)
    elif metric == 'mean':
        result = np.mean(bound_data, axis=0)
    elif metric == 'sum':
        result = np.sum(bound_data, axis=0)
    else:
        raise ValueError("Metric " + str(metric) + " is not recognized.")

    return result