Python units.Quantity类代码示例

def make_image():
    table = Table.read('acceptance_curve.fits')
    table.pprint()
    center = SkyCoord(83.63, 22.01, unit='deg').galactic

    counts_image = make_empty_image(nxpix=1000, nypix=1000, binsz=0.01, xref=center.l.deg, yref=center.b.deg,
                                    proj='TAN')
    bkg_image = counts_image.copy()
    data_store = DataStore.from_dir('$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2')

    for events in data_store.load_all("events"):
        center = events.pointing_radec.galactic
        livetime = events.observation_live_time_duration
        solid_angle = Angle(0.01, "deg") ** 2

        counts_image.data += bin_events_in_image(events, counts_image).data

        #interp_param = dict(bounds_error=False, fill_value=None)

        acc_hdu = fill_acceptance_image(bkg_image.header, center, table["offset"], table["Acceptance"])
        acc = Quantity(acc_hdu.data, table["Acceptance"].unit) * solid_angle * livetime
        bkg_image.data += acc.decompose()
        print(acc.decompose().sum())

    counts_image.writeto("counts_image.fits", clobber=True)
    bkg_image.writeto("bkg_image.fits", clobber=True)

class MyQuantity(Quantity):

    def __init__(self,coord,unit):
        self.q = Quantity(coord,unit)

    @property
    def value(self):
        return self.q.value

    @property
    def unit(self):
        return self.q.unit

    def set_value(self,value):
        assert isinstance(value,(float,int))
        _unit = self.q.unit
        self.q = Quantity(value,_unit)

    def change_unit(self,unit):
        assert isinstance(unit,(str,Unit))
        _newQ = self.q.to(unit)
        self.q = _newQ

    def asunit(self,unit):
        _q = self.q.to(unit,equivalencies=spectral())
        return _q

    def evaluate(self, x):
        """Wrapper around the evaluate method on the Astropy model classes.

        Parameters
        ----------
        x : `~gammapy.utils.energy.Energy`
            Evaluation point
        """
        if self.spectral_model == 'PowerLaw':
            flux = models.PowerLaw1D.evaluate(x, self.parameters.norm,
                                              self.parameters.reference,
                                              self.parameters.index)
        elif self.spectral_model == 'LogParabola':
            # LogParabola evaluation does not work with arrays because
            # there is bug when using '**' with Quantities
            # see https://github.com/astropy/astropy/issues/4764
            flux = Quantity([models.LogParabola1D.evaluate(xx,
                                                           self.parameters.norm,
                                                           self.parameters.reference,
                                                           self.parameters.alpha,
                                                           self.parameters.beta)
                             for xx in x])
        else:
            raise NotImplementedError('Not implemented for model {}.'.format(self.spectral_model))

        return flux.to(self.parameters.norm.unit)

    def logspace(cls, vmin, vmax, nbins, unit=None):
        """Create axis with equally log-spaced nodes

        if no unit is given, it will be taken from vmax

        Parameters
        ----------
        vmin : `~astropy.units.Quantity`, float
            Lowest value
        vmax : `~astropy.units.Quantity`, float
            Highest value
        bins : int
            Number of bins
        unit : `~astropy.units.UnitBase`, str
            Unit
        """

        if unit is not None:
            vmin = Quantity(vmin, unit)
            vmax = Quantity(vmax, unit)
        else:
            vmin = Quantity(vmin)
            vmax = Quantity(vmax)
            unit = vmax.unit
            vmin = vmin.to(unit)

        x_min, x_max = np.log10([vmin.value, vmax.value])
        vals = np.logspace(x_min, x_max, nbins)

        return cls(vals * unit, interpolation_mode='log')

def solid_angle(image):
    """Compute the solid angle of each pixel.

    This will only give correct results for CAR maps!

    Parameters
    ----------
    image : `~astropy.io.fits.ImageHDU`
        Input image

    Returns
    -------
    area_image : `~astropy.units.Quantity`
        Solid angle image (matching the input image) in steradians.
    """
    # Area of one pixel at the equator
    cdelt0 = image.header['CDELT1']
    cdelt1 = image.header['CDELT2']
    equator_area = Quantity(abs(cdelt0 * cdelt1), 'deg2')

    # Compute image with fraction of pixel area at equator
    glat = coordinates(image)[1]
    area_fraction = np.cos(np.radians(glat))

    result = area_fraction * equator_area.to('sr')

    return result

def clean_up(table_in):
    """Create a new table with exactly the columns / info we want.
    """
    table = Table()
    
    v = Quantity(table_in['v'].data, table_in['v'].unit)
    energy = v.to('MeV', equivalencies=u.spectral())
    table['energy'] = energy
    
    #vFv = Quantity(table['vFv'].data, table['vFv'].unit)
    #flux = (vFv / (energy ** 2)).to('cm^-2 s^-1 MeV^-1')
    
    table['energy_flux'] = table_in['vFv']
    table['energy_flux_err_lo'] = table_in['ed_vFv'] 
    table['energy_flux_err_hi'] = table_in['eu_vFv']
    # Compute symmetrical error because most chi^2 fitters
    # can't handle asymmetrical Gaussian erros
    table['energy_flux_err'] = 0.5 * (table_in['eu_vFv'] + table_in['ed_vFv'])

    table['component'] = table_in['component']
    table['paper'] = table_in['paper']
    
    mask = table['energy_flux_err_hi'] == 0

    return table

    def make_bkg_cube(self, bkg_norm=True):
        """
        Make the background image for one observation from a bkg model.

        Parameters
        ----------
        bkg_norm : bool
            If true, apply the scaling factor from the number of counts
            outside the exclusion region to the bkg image
        """
        for i_E in range(len(self.energy_reco) - 1):
            energy_band = Energy(
                [self.energy_reco[i_E].value, self.energy_reco[i_E + 1].value],
                self.energy_reco.unit)
            table = self.bkg.acceptance_curve_in_energy_band(
                energy_band=energy_band)
            center = self.obs_center.galactic
            bkg_hdu = fill_acceptance_image(self.header, center,
                                            table["offset"],
                                            table["Acceptance"],
                                            self.offset_band[1])
            bkg_image = Quantity(bkg_hdu.data, table[
                "Acceptance"].unit) * self.bkg_cube.sky_image_ref.solid_angle() * self.livetime
            self.bkg_cube.data[i_E, :, :] = bkg_image.decompose().value

        if bkg_norm:
            scale = self.background_norm_factor()
            self.bkg_cube.data = scale * self.bkg_cube.data
            if self.save_bkg_scale:
                self.table_bkg_scale.add_row([self.obs_id, scale])

def test_saturation_water_pressure():

    args_list = [
        (1.e-30, None, apu.K),
        (1.e-30, None, apu.hPa),
        ]
    check_astro_quantities(atm.saturation_water_pressure, args_list)

    temp = Quantity([100, 200, 300], apu.K)
    press = Quantity([900, 1000, 1100], apu.hPa)

    press_w = Quantity(
        [2.57439748e-17, 3.23857740e-03, 3.55188758e+01], apu.hPa
        )

    assert_quantity_allclose(
        atm.saturation_water_pressure(temp, press),
        press_w
        )
    assert_quantity_allclose(
        atm.saturation_water_pressure(
            temp.to(apu.mK), press.to(apu.Pa)
            ),
        press_w
        )

def test_LogEnergyAxis():
    from scipy.stats import gmean
    energy = Quantity([1, 10, 100], 'TeV')
    energy_axis = LogEnergyAxis(energy)

    energy = Quantity(gmean([1, 10]), 'TeV')
    pix = energy_axis.wcs_world2pix(energy.to('MeV'))
    assert_allclose(pix, 0.5)

    world = energy_axis.wcs_pix2world(pix)
    assert_quantity_allclose(world, energy)

 def yindex(self, index):
     if index is None:
         del self.yindex
         return
     if not isinstance(index, Quantity):
         index = Quantity(index, unit=self._default_yunit, copy=False)
     self.y0 = index[0]
     index.regular = is_regular(index.value)
     if index.regular:
         self.dy = index[1] - index[0]
     else:
         del self.dy
     self._yindex = index

    def _get_range_from_textfields(self, min_text, max_text, linelist_units, plot_units):
        amin = amax = None
        if min_text.hasAcceptableInput() and max_text.hasAcceptableInput():

            amin = float(min_text.text())
            amax = float(max_text.text())

            amin = Quantity(amin, plot_units)
            amax = Quantity(amax, plot_units)

            amin = amin.to(linelist_units, equivalencies=u.spectral())
            amax = amax.to(linelist_units, equivalencies=u.spectral())

        return (amin, amax)

    def evaluate(self, method=None, **kwargs):
        """Evaluate NDData Array

        This function provides a uniform interface to several interpolators.
        The evaluation nodes are given as ``kwargs``.

        Currently available:
        `~scipy.interpolate.RegularGridInterpolator`, methods: linear, nearest

        Parameters
        ----------
        method : str {'linear', 'nearest'}, optional
            Interpolation method
        kwargs : dict
            Keys are the axis names, Values the evaluation points

        Returns
        -------
        array : `~astropy.units.Quantity`
            Interpolated values, axis order is the same as for the NDData array
        """
        values = []
        for axis in self.axes:
            # Extract values for each axis, default: nodes
            temp = Quantity(kwargs.pop(axis.name, axis.nodes))
            # Transform to correct unit
            temp = temp.to(axis.unit).value
            # Transform to match interpolation behaviour of axis
            values.append(np.atleast_1d(axis._interp_values(temp)))

        # This is to catch e.g. typos in axis names
        if kwargs != {}:
            raise ValueError("Input given for unknown axis: {}".format(kwargs))

        if method is None:
            out = self._eval_regular_grid_interp(values)
        elif method == 'linear':
            out = self._eval_regular_grid_interp(values, method='linear')
        elif method == 'nearest':
            out = self._eval_regular_grid_interp(values, method='nearest')
        else:
            raise ValueError('Interpolator {} not available'.format(method))

        # Clip interpolated values to be non-negative
        np.clip(out, 0, None, out=out)
        # Attach units to the output
        out = out * self.data.unit

        return out

    def __init__(self, energy, rad, exposure=None, psf_value=None):

        self.energy = Quantity(energy).to('GeV')
        self.rad = Quantity(rad).to('radian')
        if exposure is None:
            self.exposure = Quantity(np.ones(len(energy)), 'cm^2 s')
        else:
            self.exposure = Quantity(exposure).to('cm^2 s')

        if psf_value is None:
            self.psf_value = Quantity(np.zeros(len(energy), len(rad)), 'sr^-1')
        else:
            self.psf_value = Quantity(psf_value).to('sr^-1')

        # Cache for TablePSF at each energy ... only computed when needed
        self._table_psf_cache = [None] * len(self.energy)

def parse_value(value, default_units, equivalence=None):
    if isinstance(value, string_types):
        v = value.split(",")
        if len(v) == 2:
            value = (float(v[0]), v[1])
        else:
            value = float(v[0])
    if hasattr(value, "to_astropy"):
        value = value.to_astropy()
    if isinstance(value, Quantity):
        q = Quantity(value.value, value.unit)
    elif iterable(value):
        q = Quantity(value[0], value[1])
    else:
        q = Quantity(value, default_units)
    return q.to(default_units, equivalencies=equivalence).value

    def __init__(
        self,
        energy_lo,
        energy_hi,
        theta,
        sigmas,
        norms,
        energy_thresh_lo="0.1 TeV",
        energy_thresh_hi="100 TeV",
    ):
        self.energy_lo = Quantity(energy_lo, "TeV")
        self.energy_hi = Quantity(energy_hi, "TeV")
        ebounds = EnergyBounds.from_lower_and_upper_bounds(
            self.energy_lo, self.energy_hi
        )
        self.energy = ebounds.log_centers
        self.theta = Quantity(theta, "deg")
        sigmas[0][sigmas[0] == 0] = 1
        sigmas[1][sigmas[1] == 0] = 1
        sigmas[2][sigmas[2] == 0] = 1
        self.sigmas = sigmas

        self.norms = norms
        self.energy_thresh_lo = Quantity(energy_thresh_lo, "TeV")
        self.energy_thresh_hi = Quantity(energy_thresh_hi, "TeV")

        self._interp_norms = self._setup_interpolators(self.norms)
        self._interp_sigmas = self._setup_interpolators(self.sigmas)