Python astropy.units.deg() Examples

def test_pixel_resolution_to_nside():

    # Check the different rounding options
    nside = pixel_resolution_to_nside(13 * u.arcmin, round='nearest')
    assert nside == 256

    nside = pixel_resolution_to_nside(13 * u.arcmin, round='up')
    assert nside == 512

    nside = pixel_resolution_to_nside(13 * u.arcmin, round='down')
    assert nside == 256

    # Check that it works with arrays
    nside = pixel_resolution_to_nside([1e3, 10, 1e-3] * u.deg, round='nearest')
    assert_equal(nside, [1, 8, 65536])

    with pytest.raises(ValueError) as exc:
        pixel_resolution_to_nside(13 * u.arcmin, round='peaches')
    assert exc.value.args[0] == "Invalid value for round: 'peaches'"

    with pytest.raises(AttributeError) as exc:
        pixel_resolution_to_nside(13)
    assert exc.value.args[0] == "'int' object has no attribute 'to'" 

def test_equ_random_sample_nodist_vector(self):
        """
        Test that a random sample of the reddening vs. distance curve is drawn
        from the full set of samples. Uses vector of coordinates as input.
        """

        # Prepare coordinates
        l = [d['l']*units.deg for d in self._test_data]
        b = [d['b']*units.deg for d in self._test_data]

        c = coords.SkyCoord(l, b, frame='galactic')

        ebv_data = np.array([d['samples'] for d in self._test_data])
        ebv_calc = self._bayestar(c, mode='random_sample')

        # print 'vector random sample:'
        # print 'ebv_data.shape = {}'.format(ebv_data.shape)
        # print 'ebv_calc.shape = {}'.format(ebv_calc.shape)
        # print ebv_data[0]
        # print ebv_calc[0]

        d_ebv = np.min(np.abs(ebv_data[:,:,:] - ebv_calc[:,None,:]), axis=1)
        np.testing.assert_allclose(d_ebv, 0., atol=0.001, rtol=0.0001) 

def test_equ_med_scalar(self):
        """
        Test that median reddening is correct in at arbitary distances, using
        individual coordinates as input.
        """
        for d in self._test_data:
            l = d['l']*units.deg
            b = d['b']*units.deg

            for reps in range(10):
                dm = 3. + (25.-3.)*np.random.random()
                dist = 10.**(dm/5.-2.)
                c = coords.SkyCoord(l, b, distance=dist*units.kpc, frame='galactic')

                ebv_samples = self._interp_ebv(d, dist)
                ebv_data = np.nanmedian(ebv_samples)
                ebv_calc = self._bayestar(c, mode='median')

                np.testing.assert_allclose(ebv_data, ebv_calc, atol=0.001, rtol=0.0001) 

def test_bounds(self):
        """
        Test that out-of-bounds coordinates return NaN reddening, and that
        in-bounds coordinates do not return NaN reddening.
        """

        for mode in (['random_sample', 'random_sample_per_pix',
                      'median', 'samples', 'mean']):
            # Draw random coordinates, both above and below dec = -30 degree line
            n_pix = 1000
            ra = -180. + 360.*np.random.random(n_pix)
            dec = -75. + 90.*np.random.random(n_pix)    # 45 degrees above/below
            c = coords.SkyCoord(ra, dec, frame='icrs', unit='deg')

            ebv_calc = self._bayestar(c, mode=mode)

            nan_below = np.isnan(ebv_calc[dec < -35.])
            nan_above = np.isnan(ebv_calc[dec > -25.])
            pct_nan_above = np.sum(nan_above) / float(nan_above.size)

            # print r'{:s}: {:.5f}% nan above dec=-25 deg.'.format(mode, 100.*pct_nan_above)

            self.assertTrue(np.all(nan_below))
            self.assertTrue(pct_nan_above < 0.05) 

def test_equ_med_far_vector(self):
        """
        Test that median reddening is correct in the far limit, using a vector
        of coordinates as input.
        """
        l = [d['l']*units.deg for d in self._test_data]
        b = [d['b']*units.deg for d in self._test_data]
        dist = [1.e3*units.kpc for bb in b]
        c = coords.SkyCoord(l, b, distance=dist, frame='galactic')

        ebv_data = np.array([np.nanmedian(d['samples'][:,-1]) for d in self._test_data])
        ebv_calc = self._bayestar(c, mode='median')

        # print 'vector:'
        # print r'% residual:'
        # for ed,ec in zip(ebv_data, ebv_calc):
        #     print '  {: >8.3f}'.format((ec - ed) / (0.02 + 0.02 * ed))

        np.testing.assert_allclose(ebv_data, ebv_calc, atol=0.001, rtol=0.0001) 

def test_fEZ(self):
        """
        Test that fEZ returns correct shape and units.
        """
        exclude_mods=[]

        for mod in self.allmods:
            if mod.__name__ in exclude_mods:
                continue
            if 'fEZ' in mod.__dict__:
                obj = mod()
                # use 3 planets
                d = 10.*np.random.rand(3)*u.AU
                I = np.random.uniform(0.0, 180.0, 3)*u.deg
                fEZs = obj.fEZ(self.TL.MV[0], I, d)
                self.assertEqual(len(fEZs), 3, 'fEZ does not return same number of values as planets tested for {}'.format(mod.__name__))
                self.assertEqual(fEZs.unit, self.unit, 'fEZ does not return 1/arcsec**2 for {}'.format(mod.__name__)) 

def test_shape(self):
        """
        Test that the output shapes are as expected with input coordinate arrays
        of different shapes.
        """

        for mode in ['random_sample', 'median', 'mean', 'samples']:
            for reps in range(5):
                # Draw random coordinates, with different shapes
                n_dim = np.random.randint(1,4)
                shape = np.random.randint(1,7, size=(n_dim,))

                ra = -180. + 360.*np.random.random(shape)
                dec = -90. + 180. * np.random.random(shape)
                c = coords.SkyCoord(ra, dec, frame='icrs', unit='deg')

                ebv_calc = self._bayestar(c, mode=mode)

                np.testing.assert_equal(ebv_calc.shape[:n_dim], shape)

                if mode == 'samples':
                    self.assertEqual(len(ebv_calc.shape), n_dim+2) # sample, distance
                else:
                    self.assertEqual(len(ebv_calc.shape), n_dim+1) # distance 

def test_equ_random_sample_scalar(self):
        """
        Test that random sample of reddening at arbitary distance is actually
        from the set of possible reddening samples at that distance. Uses vector
        of coordinates/distances as input. Uses single set of
        coordinates/distance as input.
        """
        for d in self._test_data:
            # Prepare coordinates (with random distances)
            l = d['l']*units.deg
            b = d['b']*units.deg
            dm = 3. + (25.-3.)*np.random.random()

            dist = 10.**(dm/5.-2.)
            c = coords.SkyCoord(l, b, distance=dist*units.kpc, frame='galactic')

            ebv_data = self._interp_ebv(d, dist)
            ebv_calc = self._bayestar(c, mode='random_sample')

            d_ebv = np.min(np.abs(ebv_data[:] - ebv_calc))

            np.testing.assert_allclose(d_ebv, 0., atol=0.001, rtol=0.0001) 

def test_shape(self):
        """
        Test that the output shapes are as expected with input coordinate arrays
        of different shapes.
        """

        for reps in range(5):
            # Draw random coordinates, with different shapes
            n_dim = np.random.randint(1,4)
            shape = np.random.randint(1,7, size=(n_dim,))

            ra = (-180. + 360.*np.random.random(shape)) * units.deg
            dec = (-90. + 180. * np.random.random(shape)) * units.deg
            c = coords.SkyCoord(ra, dec, frame='icrs')

            E = self._planck(c)

            np.testing.assert_equal(E.shape, shape) 

def test_interpolate_bilinear_invalid():
    values = np.ones(133)
    with pytest.raises(ValueError) as exc:
        interpolate_bilinear_lonlat([1, 3, 4] * u.deg, [3, 2, 6] * u.deg, values)
    assert exc.value.args[0] == 'Number of pixels must be divisible by 12'

    values = np.ones(192)
    with pytest.raises(ValueError) as exc:
        interpolate_bilinear_lonlat([1, 3, 4] * u.deg, [3, 2, 6] * u.deg,
                                    values, order='banana')
    assert exc.value.args[0] == "order must be 'nested' or 'ring'"

    result = interpolate_bilinear_lonlat([0, np.nan] * u.deg,
                                         [0, np.nan] * u.deg, values,
                                         order='nested')
    assert result.shape == (2,)
    assert result[0] == 1
    assert np.isnan(result[1]) 

def test_equ_samples_nodist_vector(self):
        """
        Test that full set of samples of reddening vs. distance curves is
        correct. Uses vector of coordinates as input.
        """

        # Prepare coordinates
        l = [d['l']*units.deg for d in self._test_data]
        b = [d['b']*units.deg for d in self._test_data]

        c = coords.SkyCoord(l, b, frame='galactic')

        ebv_data = np.array([d['samples'] for d in self._test_data])
        ebv_calc = self._bayestar(c, mode='samples')

        # print 'vector random sample:'
        # print 'ebv_data.shape = {}'.format(ebv_data.shape)
        # print 'ebv_calc.shape = {}'.format(ebv_calc.shape)
        # print ebv_data[0]
        # print ebv_calc[0]

        np.testing.assert_allclose(ebv_data, ebv_calc, atol=0.001, rtol=0.0001) 

def test_spiral_rot_rate():
    # Check that a slower rotation rate results in a less wound spiral
    v = 100 * u.km / u.s
    r0 = 0 * u.au
    l0 = 0 * u.deg
    spiral1 = ParkerSpiral(v, r0, l0)
    spiral2 = ParkerSpiral(v, r0, l0, 14 * u.deg / u.day)

    long1 = spiral1.longitude(1 * u.au)
    long2 = spiral1.longitude(2 * u.au)
    assert long1 > long2 

def fetch_alerts(self, parameters):
        """Must return an iterator"""
        response = requests.get(BROKER_URL)
        response.raise_for_status()

        html_data = response.text.split('\n')
        for line in html_data:
            if 'var alerts' in line:
                alerts_data = line.replace('var alerts = ', '')
                alerts_data = alerts_data.replace('\n', '').replace(';', '')

        alert_list = json.loads(alerts_data)

        if parameters['cone'] is not None and len(parameters['cone']) > 0:
            cone_params = parameters['cone'].split(',')
            parameters['cone_ra'] = float(cone_params[0])
            parameters['cone_dec'] = float(cone_params[1])
            parameters['cone_radius'] = float(cone_params[2])*u.deg
            parameters['cone_centre'] = SkyCoord(float(cone_params[0]),
                                                 float(cone_params[1]),
                                                 frame="icrs", unit="deg")

        filtered_alerts = []
        if parameters['target_name'] is not None and \
                len(parameters['target_name']) > 0:
            for alert in alert_list:
                if parameters['target_name'] in alert['name']:
                    filtered_alerts.append(alert)

        elif 'cone_radius' in parameters.keys():
            for alert in alert_list:
                c = SkyCoord(float(alert['ra']), float(alert['dec']),
                             frame="icrs", unit="deg")
                if parameters['cone_centre'].separation(c) <= parameters['cone_radius']:
                    filtered_alerts.append(alert)

        else:
            filtered_alerts = alert_list

        return iter(filtered_alerts) 

def __init__(self, probe):
        self.probe = _check_probe(probe)
        self.units = OrderedDict([
            ('B instrument', u.dimensionless_unscaled),
            ('Bx', u.nT), ('By', u.nT), ('Bz', u.nT),
            ('sigma B', u.nT),
            ('Ion instrument', u.dimensionless_unscaled),
            ('Status', u.dimensionless_unscaled),
            ('Tp_par', u.K), ('Tp_perp', u.K),
            ('carrot', u.dimensionless_unscaled),
            ('r_sun', u.AU), ('clat', u.deg),
            ('clong', u.deg), ('earth_he_angle', u.deg),
            ('n_p', u.cm**-3), ('vp_x', u.km / u.s),
            ('vp_y', u.km / u.s), ('vp_z', u.km / u.s),
            ('vth_p_par', u.km / u.s), ('vth_p_perp', u.km / u.s)]) 

def __init__(self, v, r0, l0, omega_sun=14.713 * (u.deg / u.day)):
        self.v = v
        self.r0 = r0
        self.l0 = l0
        self.omega_sun = omega_sun 

def test_parker_spiral():
    spiral = ParkerSpiral(100 * u.km / u.s, 0 * u.au, 0 * u.deg)
    longs = spiral.longitude([1, 2] * u.au)
    assert_quantity_allclose(longs[0], -254.7492444 * u.deg)
    assert_quantity_allclose(longs[1], -509.4984888 * u.deg) 

def to_target(self):
        target = super().to_target()
        target.type = 'SIDEREAL'
        target.name = (self.catalog_data['name_prefix'] + self.catalog_data['name'])
        c = SkyCoord('{0} {1}'.format(self.catalog_data['ra'], self.catalog_data['dec']), unit=(u.hourangle, u.deg))
        target.ra, target.dec = c.ra.deg, c.dec.deg
        return target 

def cloud_attenuation(lat, lon, el, f, p):
    """
    A method to estimate the attenuation due to clouds along slant paths for
    a given probability.


    Parameters
    ----------
    lat : number, sequence, or numpy.ndarray
        Latitudes of the receiver points
    lon : number, sequence, or numpy.ndarray
        Longitudes of the receiver points
    el : number, sequence, or numpy.ndarray
        Elevation angle of the receiver points (deg)
    f : number
        Frequency (GHz)
    p : number
         Percentage of time exceeded for p% of the average year


    Returns
    -------
    p: numpy.ndarray
        Rainfall rate exceeded for p% of the average year



    References
    ----------
    [1] Attenuation due to clouds and fog:
    https://www.itu.int/rec/R-REC-P.840/en
    """
    global __model
    type_output = type(lat)
    lat = prepare_input_array(lat)
    lon = prepare_input_array(lon)
    lon = np.mod(lon, 360)
    el = prepare_quantity(el, u.deg, 'Elevation angle')
    f = prepare_quantity(f, u.GHz, 'Frequency')
    val = __model.cloud_attenuation(lat, lon, el, f, p)
    return prepare_output_array(val, type_output) * u.dB 

def __init__(self, probe):
        _check_probe(probe, ['8'])
        self.probe = probe

        self.units = OrderedDict(
            [('sw_flag', u.dimensionless_unscaled),
             ('x_gse', u.R_earth), ('y_gse', u.R_earth),
             ('z_gse', u.R_earth), ('y_gsm', u.R_earth),
             ('z_gsm', u.R_earth), ('Nm', u.dimensionless_unscaled),
             ('<|B|>', u.nT), ('|<B>|', u.nT), ('<B_lat>', u.nT),
             ('<B_long>', u.nT), ('Bx_gse', u.nT), ('By_gse', u.nT),
             ('Bz_gse', u.nT), ('By_gsm', u.nT), ('Bz_gsm', u.nT),
             ('sigma|B|', u.nT), ('sigma B', u.nT),
             ('sigma B_x', u.nT), ('sigma B_y', u.nT),
             ('sigma B_z', u.nT),
             ('plas_reg', u.dimensionless_unscaled),
             ('Npp', u.dimensionless_unscaled),
             ('v_fit', u.km / u.s), ('vx_fit_gse', u.km / u.s),
             ('vy_fit_gse', u.km / u.s), ('vz_fit_gse', u.km / u.s),
             ('vlong_fit', u.deg), ('vlat_fit', u.deg),
             ('np_fit', u.cm**-3), ('Tp_fit', u.K),
             ('v_mom', u.km / u.s), ('vx_mom_gse', u.km / u.s),
             ('vy_mom_gse', u.km / u.s), ('vz_mom_gse', u.km / u.s),
             ('vlong_mom', u.deg), ('vlat_mom', u.deg),
             ('np_mom', u.cm**-3), ('Tp_mom', u.K),
             ('FCp', u.dimensionless_unscaled),
             ('DWp', u.dimensionless_unscaled)]) 

def test_bilinear_interpolation_weights(order):

    indices, weights = bilinear_interpolation_weights(100 * u.deg, 10 * u.deg,
                                                      nside=4, order=order)
    if order == 'nested':
        indices = nested_to_ring(indices, nside=4)
    assert_equal(indices, [76, 77, 60, 59])
    assert_allclose(weights, [0.532723, 0.426179, 0.038815, 0.002283], atol=1e-6) 

def test_cone_search_skycoord(self):
        coord = SkyCoord(1 * u.deg, 4 * u.deg, frame='galactic')
        result1 = self.pix.cone_search_skycoord(coord, 1 * u.deg)
        assert len(result1) == 77
        result2 = self.pix.cone_search_skycoord(coord.icrs, 1 * u.deg)
        assert_allclose(result1, result2) 

def test_interpolate_bilinear_skycoord(self):
        values = np.ones(12 * 256 ** 2) * 3
        coord = SkyCoord([1, 2, 3] * u.deg, [4, 3, 1] * u.deg, frame='fk4')
        result = self.pix.interpolate_bilinear_skycoord(coord, values)
        assert_allclose(result, [3, 3, 3])

        # Make sure that coordinate system is correctly taken into account

        values = np.arange(12 * 256 ** 2) * 3
        coord = SkyCoord([1, 2, 3] * u.deg, [4, 3, 1] * u.deg, frame='fk4')

        result1 = self.pix.interpolate_bilinear_skycoord(coord, values)
        result2 = self.pix.interpolate_bilinear_skycoord(coord.icrs, values)

        assert_allclose(result1, result2) 

def test_cone_search_lonlat_invalid(self):
        lon, lat = [1, 2] * u.deg, [3, 4] * u.deg
        with pytest.raises(ValueError) as exc:
            self.pix.cone_search_lonlat(lon, lat, 1 * u.deg)
        assert exc.value.args[0] == ('The longitude, latitude and radius must '
                                     'be scalar Quantity objects') 

def test_cone_search_lonlat(self):
        lon, lat = 1 * u.deg, 4 * u.deg
        result = self.pix.cone_search_lonlat(lon, lat, 1 * u.deg)
        assert len(result) == 77 

def test_interpolate_bilinear_lonlat(self):
        values = np.ones(12 * 256 ** 2) * 3
        result = self.pix.interpolate_bilinear_lonlat([1, 3, 4] * u.deg,
                                                      [3, 2, 6] * u.deg, values)
        assert_allclose(result, [3, 3, 3]) 

def test_bilinear_interpolation_weights(self):
        indices, weights = self.pix.bilinear_interpolation_weights([1, 3, 4] * u.deg,
                                                                   [3, 2, 6] * u.deg)
        assert indices.shape == (4, 3)
        assert weights.shape == (4, 3) 

def healpix_cone_search(lon, lat, radius, nside, order='ring'):
    """
    Find all the HEALPix pixels within a given radius of a longitude/latitude.

    Note that this returns all pixels that overlap, including partially, with
    the search cone. This function can only be used for a single lon/lat pair at
    a time, since different calls to the function may result in a different
    number of matches.

    Parameters
    ----------
    lon, lat : :class:`~astropy.units.Quantity`
        The longitude and latitude to search around
    radius : :class:`~astropy.units.Quantity`
        The search radius
    nside : int
        Number of pixels along the side of each of the 12 top-level HEALPix tiles
    order : { 'nested' | 'ring' }
        Order of HEALPix pixels

    Returns
    -------
    healpix_index : `~numpy.ndarray`
        1-D array with all the matching HEALPix pixel indices.
    """

    lon = lon.to_value(u.deg)
    lat = lat.to_value(u.deg)
    radius = radius.to_value(u.deg)

    _validate_nside(nside)
    order = _validate_order(order)

    return _core.healpix_cone_search(lon, lat, radius, nside, order) 

def nside2pixarea(nside, degrees=False):
    """Drop-in replacement for healpy `~healpy.pixelfunc.nside2pixarea`."""
    area = nside_to_pixel_area(nside)
    if degrees:
        return area.to(u.deg ** 2).value
    else:
        return area.to(u.sr).value 

def _lonlat_to_healpy(lon, lat, lonlat=False):
    # We use in-place operations below to avoid making temporary arrays - this
    # is safe because the lon/lat arrays returned from healpix_to_lonlat are
    # new and not used elsewhere.
    if lonlat:
        return lon.to(u.deg).value, lat.to(u.deg).value
    else:
        lat, lon = lat.to(u.rad).value, lon.to(u.rad).value
        if np.isscalar(lon):
            return PI_2 - lat, lon
        else:
            lat = np.subtract(PI_2, lat, out=lat)
            return lat, lon 

def log_occulterResults(self,DRM,slewTimes,sInd,sd,dV):
        """Updates the given DRM to include occulter values and results
        
        Args:
            DRM (dict):
                Design Reference Mission, contains the results of one complete
                observation (detection and characterization)
            slewTimes (astropy Quantity):
                Time to transfer to new star line of sight in units of days
            sInd (integer):
                Integer index of the star of interest
            sd (astropy Quantity):
                Angular separation between stars in rad
            dV (astropy Quantity):
                Delta-V used to transfer to new star line of sight in units of m/s
                
        Returns:
            dict:
                Design Reference Mission dictionary, contains the results of one complete
                observation (detection and characterization)
        
        """
        
        DRM['slew_time'] = slewTimes.to('day')
        DRM['slew_angle'] = sd.to('deg')
        
        slew_mass_used = slewTimes*self.defburnPortion*self.flowRate
        DRM['slew_dV'] = (slewTimes*self.ao*self.defburnPortion).to('m/s')
        DRM['slew_mass_used'] = slew_mass_used.to('kg')
        self.scMass = self.scMass - slew_mass_used
        DRM['scMass'] = self.scMass.to('kg')
        
        return DRM