Python astropy.units.Quantity() Examples

def prepare_quantity(value, units=None, name_val=None):
    """ The function verifys that a
    """
    if value is None:
        return None

    if isinstance(value, u.Quantity):
        if units in [u.K, u.deg_C, u.Kelvin, u.Celsius, u.imperial.deg_F]:
            return value.to(units, equivalencies=u.temperature()).value
        else:
            return value.to(units).value

    elif isinstance(value, numbers.Number) and units is not None:
        return value
    elif isinstance(value, np.ndarray) and units is not None:
        return value
    elif isinstance(value, list) and units is not None:
        return np.array([prepare_quantity(v, units, name_val) for v in value])
    elif isinstance(value, tuple) and units is not None:
        return np.array([prepare_quantity(v, units, name_val) for v in value])
    else:
        raise ValueError('%s has not the correct format. It must be a value,'
                         'sequence, array, or a Quantity with %s units' %
                         (name_val, str(units))) 

def absmag_to_pc(absmag, mag):
    """
    To convert absolute magnitude to parsec, Magic Number will be preserved

    :param absmag: absolute magnitude
    :type absmag: Union[float, ndarray]
    :param mag: apparent magnitude
    :type mag: Union[float, ndarray]
    :return: parsec
    :rtype: astropy Quantity
    :History: 2017-Nov-16 - Written - Henry Leung (University of Toronto)
    """
    absmag = np.array(absmag)
    mag = np.array(mag)
    magic_idx = ((absmag == MAGIC_NUMBER) | (mag == MAGIC_NUMBER))  # check for magic number

    with warnings.catch_warnings():  # suppress numpy Runtime warning caused by MAGIC_NUMBER
        warnings.simplefilter("ignore")
        pc = (10. ** (((mag - absmag) / 5.) + 1.))

    if absmag.shape != ():  # check if its only 1 element
        pc[magic_idx] = MAGIC_NUMBER
        return pc * u.parsec
    else:
        return (MAGIC_NUMBER if magic_idx == [1] else pc) * u.parsec 

def eclip2equat(self, r_eclip, currentTime):
        """Rotates heliocentric coordinates from ecliptic to equatorial frame.
        
        Args:
            r_eclip (astropy Quantity nx3 array):
                Positions vector in heliocentric ecliptic frame in units of AU
            currentTime (astropy Time array):
                Current absolute mission time in MJD
        
        Returns:
            r_equat (astropy Quantity nx3 array):
                Positions vector in heliocentric equatorial frame in units of AU
        
        """
        
        r_equat = self.equat2eclip(r_eclip, currentTime, rotsign=-1)
        
        return r_equat 

def mean_absolute_percentage_error(x, y, axis=None):
    """
    | NumPy implementation of tf.keras.metrics.mean_absolute_percentage_error with capability to deal with ``magicnumber``
      and astropy Quantity
    | Either both x and y are ndarray or both x and y are astropy.Quatity, return has no astropy units in all case

    :param x: prediction
    :type x: Union[ndarray, float, astropy.Quatity]
    :param y: ground truth
    :type y: Union[ndarray, float, astropy.Quatity]
    :param axis: NumPy axis
    :type axis: Union[NoneType, int]
    :raise: TypeError when only either x or y contains astropy units. Both x, y should carry/not carry astropy units at the same time
    :return: Mean Absolute Percentage Error
    :rtype: Union[ndarray, float]
    :History: 2018-Apr-11 - Written - Henry Leung (University of Toronto)
    """
    return mape_core(x, y, axis=axis, mode='mean') 

def get_frame_rate(self):
        """Determine the number of frames per second.

        The frame rate is calculated from the time elapsed between the
        first two frames, as inferred from their time stamps.

        Returns
        -------
        frame_rate : `~astropy.units.Quantity`
            Frames per second.
        """
        with self.temporary_offset(0):
            header0 = self.find_header()
            self.seek(header0.frame_nbytes, 1)
            header1 = self.read_header()

        # Mark 4 specification states frames-lengths range from 1.25 ms
        # to 160 ms.
        tdelta = (header1.fraction[0] - header0.fraction[0]) % 1.
        return u.Quantity(1 / tdelta, u.Hz).round() 

def test_construct_optics():
    """ create an OpticsDescription and make sure it
    fails if units are missing """
    OpticsDescription(
        name="test",
        num_mirrors=1,
        num_mirror_tiles=100,
        mirror_area=u.Quantity(550, u.m ** 2),
        equivalent_focal_length=u.Quantity(10, u.m),
    )

    with pytest.raises(TypeError):
        OpticsDescription(
            name="test",
            num_mirrors=1,
            num_mirror_tiles=100,
            mirror_area=550,
            equivalent_focal_length=10,
        ) 

def mape_core(x, y, axis=None, mode=None):
    if isinstance(x, list):
        x = np.array(x)
    if isinstance(y, list):
        y = np.array(y)
    if isinstance(x, u.Quantity) and isinstance(y, u.Quantity):
        percentage = ((x - y) / y).value
        # still need to take the value for creating mask
        x = x.value
        y = y.value
    elif (isinstance(x, u.Quantity) and not isinstance(y, u.Quantity)) or \
            (not isinstance(x, u.Quantity) and isinstance(y, u.Quantity)):
        raise TypeError("Only one of your data provided has astropy units \n"
                        "Either both x and y are ndarray or both x and y are astropy.Quatity, "
                        "return without astropy units in all case")
    else:
        percentage = (x - y) / y
    if mode == 'mean':
        return np.ma.mean(np.ma.array(np.abs(percentage) * 100., mask=((x == MAGIC_NUMBER) | (y == MAGIC_NUMBER))),
                          axis=axis)
    elif mode == 'median':
        return np.ma.median(np.ma.array(np.abs(percentage) * 100., mask=[(x == MAGIC_NUMBER) | (y == MAGIC_NUMBER)]),
                            axis=axis) 

def test_construct():
    """ Check we can make a CameraReadout from scratch """
    camera_name = "Unknown"
    sampling_rate = u.Quantity(2, u.GHz)
    reference_pulse_shape = np.ones((2, 20)).astype(np.float)
    reference_pulse_sample_width = u.Quantity(0.5, u.ns)
    readout = CameraReadout(
        camera_name=camera_name,
        sampling_rate=sampling_rate,
        reference_pulse_shape=reference_pulse_shape,
        reference_pulse_sample_width=reference_pulse_sample_width,
    )

    assert readout.camera_name == camera_name
    assert readout.sampling_rate == sampling_rate
    assert (readout.reference_pulse_shape == reference_pulse_shape).all()
    assert readout.reference_pulse_sample_width == reference_pulse_sample_width 

def get_time(self, frame_rate=None):
        """Converts ref_epoch, seconds, and frame_nr to Time object.

        Uses 'ref_epoch', which stores the number of half-years from 2000,
        and 'seconds'.  By default, it also calculates the offset using
        the current frame number.  For non-zero 'frame_nr', this requires the
        frame rate, which is calculated from the sample rate in the header.
        The latter can also be explicitly passed on.

        Parameters
        ----------
        frame_rate : `~astropy.units.Quantity`, optional
            For non-zero 'frame_nr', this is used to calculate the
            corresponding offset. If not given, the frame rate from the
            header is used (if it is non-zero).

        Returns
        -------
        time : `~astropy.time.Time`
        """
        if frame_rate is None and self['sampling_rate'] != 0:
            frame_rate = self.frame_rate
        return super().get_time(frame_rate=frame_rate) 

def stellarTeff(self, sInds):
        """Calculate the effective stellar temperature based on B-V color.
        
        This method uses the empirical fit from Ballesteros (2012) doi:10.1209/0295-5075/97/34008
        
        Args:
            sInds (integer ndarray):
                Indices of the stars of interest
        
        Returns:
            Quantity array:
                Stellar effective temperatures in degrees K
        
        """
        
        # cast sInds to array
        sInds = np.array(sInds, ndmin=1, copy=False)
        
        Teff = 4600.0*u.K * (1.0/(0.92*self.BV[sInds] + 1.7) + 1.0/(0.92*self.BV[sInds] + 0.62))
        
        return Teff 

def test_guess_area():
    x = u.Quantity([0, 1, 2], u.cm)
    y = u.Quantity([0, 0, 0], u.cm)
    n_pixels = len(x)

    geom = CameraGeometry(
        "test",
        pix_id=np.arange(n_pixels),
        pix_area=None,
        pix_x=x,
        pix_y=y,
        pix_type="rect",
    )

    assert np.all(geom.pix_area == 1 * u.cm ** 2)

    geom = CameraGeometry(
        "test",
        pix_id=np.arange(n_pixels),
        pix_area=None,
        pix_x=x,
        pix_y=y,
        pix_type="hexagonal",
    )
    assert u.allclose(geom.pix_area, 2 * np.sqrt(3) * (0.5 * u.cm) ** 2) 

def test_calc_pixel_neighbors_square_diagonal():
    """
    check that neighbors for square-pixel cameras are what we expect,
    namely that the diagonals are included if requested.
    """
    x, y = np.meshgrid(np.arange(20), np.arange(20))

    cam = CameraGeometry(
        camera_name="test",
        pix_id=np.arange(400),
        pix_type="rectangular",
        pix_x=u.Quantity(x.ravel(), u.cm),
        pix_y=u.Quantity(y.ravel(), u.cm),
        pix_area=u.Quantity(np.ones(400), u.cm ** 2),
    )

    cam._neighbors = cam.calc_pixel_neighbors(diagonal=True)
    assert set(cam.neighbors[21]) == {0, 1, 2, 20, 22, 40, 41, 42} 

def test_find_neighbor_pixels():
    """ test basic neighbor functionality """
    n_pixels = 5
    x, y = u.Quantity(
        np.meshgrid(np.linspace(-5, 5, n_pixels), np.linspace(-5, 5, n_pixels)), u.cm
    )

    geom = CameraGeometry(
        "test",
        pix_id=np.arange(n_pixels),
        pix_area=u.Quantity(4, u.cm ** 2),
        pix_x=x.ravel(),
        pix_y=y.ravel(),
        pix_type="rectangular",
    )

    neigh = geom.neighbors
    assert set(neigh[11]) == {16, 6, 10, 12} 

def mean_absolute_error(x, y, axis=None):
    """
    NumPy implementation of tf.keras.metrics.mean_absolute_error  with capability to deal with ``magicnumber``
    and astropy Quantity

    Either both x and y are ndarray or both x and y are astropy.Quatity, return without astropy units in all case

    :param x: prediction
    :type x: Union[ndarray, float, astropy.Quatity]
    :param y: ground truth
    :type y: Union[ndarray, float, astropy.Quatity]
    :param axis: NumPy axis
    :type axis: Union[NoneType, int]
    :raise: TypeError when only either x or y contains astropy units. Both x, y should carry/not carry astropy units at the same time
    :return: Mean Absolute Error
    :rtype: Union[ndarray, float]
    :History: 2018-Apr-11 - Written - Henry Leung (University of Toronto)
    """
    return mae_core(x, y, axis=axis, mode='mean') 

def mae_core(x, y, axis=None, mode=None):
    if isinstance(x, list):
        x = np.array(x)
    if isinstance(y, list):
        y = np.array(y)
    if isinstance(x, u.Quantity) and isinstance(y, u.Quantity):
        diff = (x - y).value
        # still need to take the value for creating mask
        x = x.value
        y = y.value
    elif (isinstance(x, u.Quantity) and not isinstance(y, u.Quantity)) or \
            (not isinstance(x, u.Quantity) and isinstance(y, u.Quantity)):
        raise TypeError("Only one of your data provided has astropy units \n"
                        "Either both x and y are ndarray or both x and y are astropy.Quatity, "
                        "return without astropy units in all case")
    else:
        diff = (x - y)
    if mode == 'mean':
        return np.ma.mean(np.ma.array(np.abs(diff), mask=((x == MAGIC_NUMBER) | (y == MAGIC_NUMBER))), axis=axis)
    elif mode == 'median':
        return np.ma.median(np.ma.array(np.abs(diff), mask=[(x == MAGIC_NUMBER) | (y == MAGIC_NUMBER)]), axis=axis) 

def query_equ(self, ra, dec, d=None, frame='icrs', **kwargs):
        """
        A web API version of :obj:`DustMap.query_equ()`. See the documentation for
        the corresponding local query object. Queries using Equatorial
        coordinates. By default, the ICRS frame is used, although other frames
        implemented by :obj:`astropy.coordinates` may also be specified.

        Args:
            ra (:obj:`float`, scalar or array-like): Galactic longitude, in degrees,
                or as an :obj:`astropy.unit.Quantity`.
            dec (:obj:`float`, scalar or array-like): Galactic latitude, in degrees,
                or as an :obj:`astropy.unit.Quantity`.
            d (Optional[:obj:`float`, scalar or array-like]): Distance from the Solar
                System, in kpc, or as an :obj:`astropy.unit.Quantity`. Defaults to
                ``None``, meaning no distance is specified.
            frame (Optional[icrs]): The coordinate system. Can be 'icrs' (the
                default), 'fk5', 'fk4' or 'fk4noeterms'.
            **kwargs: Any additional keyword arguments accepted by derived
                classes.

        Returns:
            The results of the query.
        """
        pass 

def query_gal(self, l, b, d=None, **kwargs):
        """
        A web API version of :obj:`DustMap.query_gal()`. See the documentation for
        the corresponding local query object. Queries using Galactic
        coordinates.

        Args:
            l (:obj:`float`, scalar or array-like): Galactic longitude, in degrees,
                or as an :obj:`astropy.unit.Quantity`.
            b (:obj:`float`, scalar or array-like): Galactic latitude, in degrees,
                or as an :obj:`astropy.unit.Quantity`.
            d (Optional[:obj:`float`, scalar or array-like]): Distance from the Solar
                System, in kpc, or as an :obj:`astropy.unit.Quantity`. Defaults to
                ``None``, meaning no distance is specified.
            **kwargs: Any additional keyword arguments accepted by derived
                classes.

        Returns:
            The results of the query.
        """
        pass 

def median_absolute_percentage_error(x, y, axis=None):
    """
    | NumPy implementation of a median version of tf.keras.metrics.mean_absolute_percentage_error with capability to
    | deal with ``magicnumber`` and astropy Quantity
    | Either both x and y are ndarray or both x and y are astropy.Quatity, return has no astropy units in all case

    :param x: prediction
    :type x: Union[ndarray, float, astropy.Quatity]
    :param y: ground truth
    :type y: Union[ndarray, float, astropy.Quatity]
    :param axis: NumPy axis
    :type axis: Union[NoneType, int]
    :raise: TypeError when only either x or y contains astropy units. Both x, y should carry/not carry astropy units at the same time
    :return: Median Absolute Percentage Error
    :rtype: Union[ndarray, float]
    :History: 2018-May-13 - Written - Henry Leung (University of Toronto)
    """
    return mape_core(x, y, axis=axis, mode='median') 

def prepare_output_array(array, type_input=None):
    """ Formats the output to have the same shape and type as the input
    """
    global output_quantity

    if isinstance(array, u.Quantity):
        value = array.value
        unit = array.unit
    else:
        value = array
        unit = None

    # Squeeze output array to remove singleton dimensions
    if type(value) in [np.ndarray, list]:
        value = np.array(value).squeeze()

    if (type_input in __NUMERIC_TYPES__ and
        (type(array) in __NUMERIC_TYPES__) or
        ((isinstance(array, np.ndarray) and array.size == 1) or
         (not type(array) not in __NUMERIC_TYPES__ and len(array) == 1))):
        value = float(value)
    elif type_input is list:
        if isinstance(value, np.ndarray):
            value = value.tolist()
        else:
            value = list(value)
    else:
        value = value

    if unit is not None:
        return value * unit
    else:
        return value 

def __init__(
        self,
        camera_name,
        sampling_rate,
        reference_pulse_shape,
        reference_pulse_sample_width,
    ):
        """Stores properties related to the readout of a Cherenkov Camera.

        Parameters
        ----------
        camera_name: str
             Camera name (e.g. NectarCam, LSTCam, ...)
        sampling_rate : u.Quantity[frequency]
            Sampling rate of the waveform
        reference_pulse_shape : ndarray
            Expected pulse shape for a signal in the waveform. 2 dimensional,
            first dimension is gain channel.
        reference_pulse_sample_width : u.Quantity[time]
            The amount of time corresponding to each sample in the 2nd
            dimension of reference_pulse_shape

        """
        self.camera_name = camera_name
        self.sampling_rate = sampling_rate
        self.reference_pulse_shape = reference_pulse_shape
        self.reference_pulse_sample_width = reference_pulse_sample_width 

def reference_pulse_sample_time(self):
        """
        Time axis for the reference pulse
        """
        _, n_samples = self.reference_pulse_shape.shape
        sample_width_ns = self.reference_pulse_sample_width.to_value(u.ns)
        pulse_max_sample = n_samples * sample_width_ns
        sample_time = np.arange(0, pulse_max_sample, sample_width_ns)
        return u.Quantity(sample_time, u.ns) 

def from_table(cls, url_or_table, **kwargs):
        """Load a CameraReadout from an `astropy.table.Table` instance or a
        file that is readable by `astropy.table.Table.read()`.

        Parameters
        ----------
        url_or_table: string or astropy.table.Table
            either input filename/url or a Table instance
        kwargs: extra keyword arguments
            extra arguments passed to `astropy.table.read()`, depending on
            file type (e.g. format, hdu, path)

        """
        tab = url_or_table
        if not isinstance(url_or_table, Table):
            tab = Table.read(url_or_table, **kwargs)

        camera_name = tab.meta.get("CAM_ID", "Unknown")
        n_channels = tab.meta["NCHAN"]
        sampling_rate = u.Quantity(tab.meta["SAMPFREQ"], u.GHz)
        reference_pulse_sample_width = u.Quantity(tab.meta["REF_WIDTH"], u.ns)
        reference_pulse_shape = np.array(
            [tab[f"reference_pulse_shape_channel{i}"] for i in range(n_channels)]
        )

        return cls(
            camera_name=camera_name,
            sampling_rate=sampling_rate,
            reference_pulse_shape=reference_pulse_shape,
            reference_pulse_sample_width=reference_pulse_sample_width,
        ) 

def deserialize(self, spectrum):
        """
        Constructs a Spectrum1D from the spectrum value stored in a ReducedDatum

        :param spectrum: JSON representation used to construct the Spectrum1D
        :type spectrum: str

        :returns: Spectrum1D representing the spectrum information
        :rtype: specutil.Spectrum1D
        """
        data = json.loads(spectrum)
        flux = Quantity(value=data['photon_flux'], unit=data['photon_flux_units'])
        wavelength = Quantity(value=data['wavelength'], unit=data['wavelength_units'])
        spectrum = Spectrum1D(flux=flux, spectral_axis=wavelength)
        return spectrum 

def set_time(self, time, frame_rate=None):
        """Converts Time object to ref_epoch, seconds, and frame_nr.

        For non-integer seconds, a frame rate is needed to calculate the
        'frame_nr'.

        Parameters
        ----------
        time : `~astropy.time.Time`
            The time to use for this header.
        frame_rate : `~astropy.units.Quantity`, optional
            For calculating 'frame_nr' from the fractional seconds.
        """
        seconds = time - self.ref_time
        int_sec = int(seconds.sec)

        # Round to nearest ns to handle timestamp difference errors.
        frac_sec = seconds - int_sec * u.s
        if abs(frac_sec) < 1. * u.ns:
            frame_nr = 0
        elif abs(1. * u.s - frac_sec) < 1. * u.ns:
            int_sec += 1
            frame_nr = 0
        else:
            if frame_rate is None:
                raise ValueError("this header does not provide a frame "
                                 "rate. Pass it in explicitly.")

            frame_nr = int((frac_sec * frame_rate).to(u.one).round())
            if abs(frame_nr / frame_rate - 1. * u.s) < 1. * u.ns:
                frame_nr = 0
                int_sec += 1

        self['seconds'] = int_sec
        self['frame_nr'] = frame_nr 

def _healpy_to_lonlat(theta, phi, lonlat=False):
    # Unlike in _lonlat_to_healpy, we don't use in-place operations since we
    # don't want to modify theta and phi since the user may be using them
    # elsewhere.
    if lonlat:
        lon = np.asarray(theta) / RAD2DEG
        lat = np.asarray(phi) / RAD2DEG
    else:
        lat = PI_2 - np.asarray(theta)
        lon = np.asarray(phi)
    return u.Quantity(lon, u.rad, copy=False), u.Quantity(lat, u.rad, copy=False) 

def get_frame_rate(self):
        """Determine the number of frames per second.

        This method first tries to determine the frame rate by looking for
        the highest frame number in the first second of data.  If that fails,
        it uses the time difference between two consecutive frames. This can
        fail if the headers do not store fractional seconds, or if the data
        rate is above 512 Mbps.

        Returns
        -------
        frame_rate : `~astropy.units.Quantity`
            Frames per second.
        """
        try:
            return super().get_frame_rate()
        except Exception as exc:
            with self.temporary_offset(0):
                try:
                    header0 = self.read_header()
                    self.seek(header0.payload_nbytes, 1)
                    header1 = self.read_header()
                    tdelta = header1.fraction - header0.fraction
                    if tdelta == 0.:
                        exc.args += ("frame rate can also not be determined "
                                     "from the first two headers, as they "
                                     "have identical fractional seconds.",)
                    return u.Quantity(1 / tdelta, u.Hz).round()
                except Exception:
                    pass
            raise exc 

def readout():
    camera_name = "Unknown"
    sampling_rate = u.Quantity(2, u.GHz)
    reference_pulse_shape = np.ones((2, 20)).astype(np.float)
    reference_pulse_sample_width = u.Quantity(0.5, u.ns)
    return CameraReadout(
        camera_name=camera_name,
        sampling_rate=sampling_rate,
        reference_pulse_shape=reference_pulse_shape,
        reference_pulse_sample_width=reference_pulse_sample_width,
    ) 

def sample_rate(self):
        """Number of complete samples per second.

        Assumes the 'sampling_rate' header field represents a per-channel
        sample rate for complex samples, or half the sample rate for real ones.
        """
        # Interprets sample rate correctly for EDV=3, but may not for EDV=1.
        return u.Quantity(self['sampling_rate']
                          * (1 if self['complex_data'] else 2),
                          u.MHz if self['sampling_unit'] else u.kHz) 

def update(self, *, time=None, frame_rate=None, sample_rate=None,
               verify=True, **kwargs):
        """Update the header by setting keywords or properties.

        Here, any keywords matching header keys are applied first, and any
        remaining ones are used to set header properties, in the order set
        by the class (in ``_properties``).

        Parameters
        ----------
        time : `~astropy.time.Time`, optional
            A possible new time.  This is updated last.
        frame_rate : `~astropy.units.Quantity`, optional
            The frame rate to use in calculating the frame number from the
            time.  Needed for times at non-integer seconds.  Ignored if no
            ``time`` is given.
        sample_rate : `~astropy.units.Quantity`, optional
            Alternative to ``frame_rate``.  Ignored if no ``time`` is given,
            or if ``frame_rate`` is given.
        verify : bool, optional
            If `True` (default), verify integrity after updating.
        **kwargs
            Arguments used to set keywords and properties.
        """
        super().update(verify=False, **kwargs)

        if time is not None:
            if frame_rate is None and sample_rate is not None:
                frame_rate = sample_rate / self.samples_per_frame

            self.set_time(time, frame_rate=frame_rate)

        if verify:
            self.verify() 

def calcfZmax(self, sInds, Obs, TL, TK, mode, hashname):
        """Finds the maximum zodiacal light values for each star over an entire orbit of the sun not including keeoput angles.
        
        Note:
            Prototype includes keepout angles because the values are all the same

        Args:
            sInds[sInds] (integer array):
                the star indicies we would like fZmax and fZmaxInds returned for
            Obs (module):
                Observatory module
            TL (module):
                Target List Module
            TK (TimeKeeping object):
                TimeKeeping object
            mode (dict):
                Selected observing mode
            hashname (string):
                hashname describing the files specific to the current json script
                
        Returns:
            tuple:
            valfZmax[sInds] (astropy Quantity array):
                the maximum fZ (for the prototype, these all have the same value) with units 1/arcsec**2
            absTimefZmax[sInds] (astropy Time array):
                returns the absolute Time the maximum fZ occurs (for the prototype, these all have the same value)
        """
        # cast sInds to array
        sInds = np.array(sInds, ndmin=1, copy=False)
        # get all array sizes
        nStars = sInds.size

        nZ = np.ones(nStars)
        valfZmax = nZ*10**(-0.4*self.magZ)/u.arcsec**2

        absTimefZmax = nZ*u.d + TK.currentTimeAbs

        return valfZmax[sInds], absTimefZmax[sInds]