Python astropy.coordinates.Angle() Examples

def circle2circle(line):
    """
    Parse a string that describes a circle in ds9 format.

    Parameters
    ----------
    line : str
        A string containing a DS9 region command for a circle.

    Returns
    -------
    circle : [ra, dec, radius]
        The center and radius of the circle.
    """
    words = re.split('[(,\s)]', line)
    ra = words[1]
    dec = words[2]
    radius = words[3][:-1]  # strip the "
    if ":" in ra:
        ra = Angle(ra, unit=u.hour)
    else:
        ra = Angle(ra, unit=u.degree)
    dec = Angle(dec, unit=u.degree)
    radius = Angle(radius, unit=u.arcsecond)
    return [ra.degree, dec.degree, radius.degree] 

def __init__(self, x, y, length, width, psi):
        """Create 2D Gaussian model for a shower image in a camera.

        Parameters
        ----------
        centroid : u.Quantity[length, shape=(2, )]
            position of the centroid of the shower in camera coordinates
        width: u.Quantity[length]
            width of shower (minor axis)
        length: u.Quantity[length]
            length of shower (major axis)
        psi : convertable to `astropy.coordinates.Angle`
            rotation angle about the centroid (0=x-axis)

        Returns
        -------
        a `scipy.stats` object

        """
        self.x = x
        self.y = y
        self.width = width
        self.length = length
        self.psi = psi 

def test_against_jpl_horizons():
    """Check that Astropy gives consistent results with the JPL Horizons example.

    The input parameters and reference results are taken from this page:
    (from the first row of the Results table at the bottom of that page)
    http://ssd.jpl.nasa.gov/?horizons_tutorial
    """
    obstime = Time('1998-07-28 03:00')
    location = EarthLocation(lon=Angle('248.405300d'),
                             lat=Angle('31.9585d'),
                             height=2.06 * u.km)
    # No atmosphere
    altaz_frame = AltAz(obstime=obstime, location=location)

    altaz = SkyCoord('143.2970d 2.6223d', frame=altaz_frame)
    radec_actual = altaz.transform_to('icrs')
    radec_expected = SkyCoord('19h24m55.01s -40d56m28.9s', frame='icrs')
    distance = radec_actual.separation(radec_expected).to('arcsec')
    # Current value: 0.238111 arcsec
    assert distance < 1 * u.arcsec 

def test_nested():

    with quantity_support():

        with quantity_support():

            fig = plt.figure()
            ax = fig.add_subplot(1, 1, 1)
            ax.scatter(Angle([1, 2, 3], u.deg), [3, 4, 5] * u.kg)

            assert ax.xaxis.get_units() == u.deg
            assert ax.yaxis.get_units() == u.kg

        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1)
        ax.scatter(Angle([1, 2, 3], u.arcsec), [3, 4, 5] * u.pc)

        assert ax.xaxis.get_units() == u.arcsec
        assert ax.yaxis.get_units() == u.pc 

def test_quantity_subclass():
    """Check that subclasses are recognized.

    This sadly is not done by matplotlib.units itself, though
    there is a PR to change it:
    https://github.com/matplotlib/matplotlib/pull/13536
    """
    plt.figure()

    with quantity_support():
        plt.scatter(Angle([1, 2, 3], u.deg), [3, 4, 5] * u.kg)
        plt.scatter([105, 210, 315] * u.arcsec, [3050, 3025, 3010] * u.g)
        plt.plot(Angle([105, 210, 315], u.arcsec), [3050, 3025, 3010] * u.g)

        assert plt.gca().xaxis.get_units() == u.deg
        assert plt.gca().yaxis.get_units() == u.kg 

def test_Ellipse2D():
    """Test Ellipse2D model."""
    amplitude = 7.5
    x0, y0 = 15, 15
    theta = Angle(45, 'deg')
    em = models.Ellipse2D(amplitude, x0, y0, 7, 3, theta.radian)
    y, x = np.mgrid[0:30, 0:30]
    e = em(x, y)
    assert np.all(e[e > 0] == amplitude)
    assert e[y0, x0] == amplitude

    rotation = models.Rotation2D(angle=theta.degree)
    point1 = [2, 0]      # Rotation2D center is (0, 0)
    point2 = rotation(*point1)
    point1 = np.array(point1) + [x0, y0]
    point2 = np.array(point2) + [x0, y0]
    e1 = models.Ellipse2D(amplitude, x0, y0, 7, 3, theta=0.)
    e2 = models.Ellipse2D(amplitude, x0, y0, 7, 3, theta=theta.radian)
    assert e1(*point1) == e2(*point2) 

def setup_class(self):
        self.impact_reco = ImPACTReconstructor(root_dir=".")
        self.horizon_frame = AltAz()

        self.h1 = HillasParametersContainer(
            x=1 * u.deg,
            y=1 * u.deg,
            r=1 * u.deg,
            phi=Angle(0 * u.rad),
            intensity=100,
            length=0.4 * u.deg,
            width=0.4 * u.deg,
            psi=Angle(0 * u.rad),
            skewness=0,
            kurtosis=0,
        )

    # @pytest.mark.skip('need a dataset for this to work') 

def generate_IRAF_input(writefn=None):
    dt = test_input_time.utc.datetime

    coos = _get_test_input_radecs()

    lines = []
    for ra, dec in zip(coos.ra, coos.dec):
        rastr = Angle(ra).to_string(u.hour, sep=':')
        decstr = Angle(dec).to_string(u.deg, sep=':')

        msg = '{yr} {mo} {day} {uth}:{utmin} {ra} {dec}'
        lines.append(msg.format(yr=dt.year, mo=dt.month, day=dt.day,
                                uth=dt.hour, utmin=dt.minute,
                                ra=rastr, dec=decstr))
    if writefn:
        with open(writefn, 'w') as f:
            for l in lines:
                f.write(l)
    else:
        for l in lines:
            print(l)
    print('Run IRAF as:\nastutil\nrvcorrect f=<filename> observatory=Paranal') 

def test_search_around_scalar():
    from astropy.coordinates import SkyCoord, Angle

    cat = SkyCoord([1, 2, 3], [-30, 45, 8], unit="deg")
    target = SkyCoord('1.1 -30.1', unit="deg")

    with pytest.raises(ValueError) as excinfo:
        cat.search_around_sky(target, Angle('2d'))

    # make sure the error message is *specific* to search_around_sky rather than
    # generic as reported in #3359
    assert 'search_around_sky' in str(excinfo.value)

    with pytest.raises(ValueError) as excinfo:
        cat.search_around_3d(target, Angle('2d'))
    assert 'search_around_3d' in str(excinfo.value) 

def get_time(self, coords):
        """Gets time, including light travel time correction to solar system barycenter for object given location"""
        t0 = self.post_obj.time - self.post_obj.barycorr
        
        ra = Angle(coords[0], u.deg)
        dec = Angle(coords[1], u.deg)
        
        greenwich = coord.EarthLocation.of_site('greenwich')
        times = time.Time(t0+2457000, format='jd',
                   scale='utc', location=greenwich)
        ltt_bary = times.light_travel_time(SkyCoord(ra, dec)).value
        
        self.time = t0 + ltt_bary
        self.barycorr = ltt_bary 

def formatCoord(x, y):
    """Custom coordinate format"""
    global midRA, midDec

    RA, Dec = factor.directions.xy2radec([x], [y], midRA, midDec)
    RA_str = Angle(RA[0], unit='deg').to_string('hour')
    Dec_str = Angle(Dec[0], unit='deg').to_string('deg')
    return 'RA = {0} Dec = {1}'.format(RA_str, Dec_str) 

def coords_from_name(name):
    """ Finds the RA, Dec for a given target name."""
    from astroquery.simbad import Simbad
    result_table = Simbad.query_object(name)
    coords = SkyCoord(Angle(result_table['RA'][0] + ' hours'),
                      Angle(result_table['DEC'][0] + ' degrees'))
    return (coords.ra.deg, coords.dec.deg) 

def _add_gaia_figure_elements(self, tpf, fig, magnitude_limit=18):
        """Make the Gaia Figure Elements"""
        # Get the positions of the Gaia sources
        c1 = SkyCoord(tpf.ra, tpf.dec, frame='icrs', unit='deg')
        # Use pixel scale for query size
        pix_scale = 21.0
        # We are querying with a diameter as the radius, overfilling by 2x.
        from astroquery.vizier import Vizier
        Vizier.ROW_LIMIT = -1
        result = Vizier.query_region(c1, catalog=["I/345/gaia2"],
                                     radius=Angle(np.max(tpf.shape[1:]) * pix_scale, "arcsec"))
        no_targets_found_message = ValueError('Either no sources were found in the query region '
                                              'or Vizier is unavailable')
        too_few_found_message = ValueError('No sources found brighter than {:0.1f}'.format(magnitude_limit))
        if result is None:
            raise no_targets_found_message
        elif len(result) == 0:
            raise too_few_found_message
        result = result["I/345/gaia2"].to_pandas()
        result = result[result.Gmag < magnitude_limit]
        if len(result) == 0:
            raise no_targets_found_message
        radecs = np.vstack([result['RA_ICRS'], result['DE_ICRS']]).T
        coords = tpf.wcs.all_world2pix(radecs, 1) ## TODO, is origin supposed to be zero or one?
        year = ((tpf.astropy_time[0].jd - 2457206.375) * u.day).to(u.year)
        pmra = ((np.nan_to_num(np.asarray(result.pmRA)) * u.milliarcsecond/u.year) * year).to(u.arcsec).value
        pmdec = ((np.nan_to_num(np.asarray(result.pmDE)) * u.milliarcsecond/u.year) * year).to(u.arcsec).value
        result.RA_ICRS += pmra
        result.DE_ICRS += pmdec

        # Gently size the points by their Gaia magnitude
        sizes = 10000.0 / 2**(result['Gmag']/2)

        plt.scatter(coords[:, 0]+tpf.column, coords[:, 1]+tpf.row, c='firebrick', alpha=0.5, edgecolors='r', s=sizes)
        plt.scatter(coords[:, 0]+tpf.column, coords[:, 1]+tpf.row, c='None', edgecolors='r', s=sizes)
        plt.xlim([tpf.column, tpf.column+tpf.shape[1]])
        plt.ylim([tpf.row, tpf.row+tpf.shape[2]])

        return fig 

def test_ra2dec():
    """Test ra2dec against astropy conversion"""
    # Test against the astropy calculations
    for ra in ['14:21:45.003', '-12 04 22', '-00 01 12.003']:
        ans = at.ra2dec(ra)
        desired = Angle(ra, unit=u.hourangle).hour * 15
        assert_approx_equal(ans, desired, "{0} != {1}".format(ans, desired)) 

def radec2azel(
    ra_deg: "ndarray", dec_deg: "ndarray", lat_deg: "ndarray", lon_deg: "ndarray", time: datetime, *, use_astropy: bool = False
) -> typing.Tuple["ndarray", "ndarray"]:
    """
    sky coordinates (ra, dec) to viewing angle (az, el)

    Parameters
    ----------
    ra_deg : "ndarray"
         ecliptic right ascension (degress)
    dec_deg : "ndarray"
         ecliptic declination (degrees)
    lat_deg : "ndarray"
              observer latitude [-90, 90]
    lon_deg : "ndarray"
              observer longitude [-180, 180] (degrees)
    time : datetime.datetime or str
           time of observation
    use_astropy : bool, optional
                 default use astropy.

    Returns
    -------
    az_deg : "ndarray"
             azimuth [degrees clockwize from North]
    el_deg : "ndarray"
             elevation [degrees above horizon (neglecting aberration)]
    """

    if use_astropy and Time is not None:
        obs = EarthLocation(lat=lat_deg * u.deg, lon=lon_deg * u.deg)
        points = SkyCoord(Angle(ra_deg, unit=u.deg), Angle(dec_deg, unit=u.deg), equinox="J2000.0")
        altaz = points.transform_to(AltAz(location=obs, obstime=Time(str2dt(time))))

        return altaz.az.degree, altaz.alt.degree

    return vradec2azel(ra_deg, dec_deg, lat_deg, lon_deg, time) 

def test_dec2dec():
    """Test dec2dec against astropy conversion"""
    # Test against the astropy calculations
    for dec in ['+14:21:45.003', '-99 04 22', '-00 01 23.456', '00 01']:
        ans = at.dec2dec(dec)
        desired = Angle(dec, unit=u.degree).degree
        assert_approx_equal(ans, desired, err_msg="{0} != {1}".format(ans, desired)) 

def from_tree(cls, node, ctx):
        angle = Angle(QuantityType.from_tree(node['ra']['wrap_angle'], ctx))
        wrap_angle = Angle(angle)
        ra = Longitude(
            node['ra']['value'],
            unit=node['ra']['unit'],
            wrap_angle=wrap_angle)
        dec = Latitude(node['dec']['value'], unit=node['dec']['unit'])

        return ICRS(ra=ra, dec=dec) 

def calculateSeparation(ra1, dec1, ra2, dec2):
    """
    Returns angular separation between two coordinates (all in degrees).

    Parameters
    ----------
    ra1 : float or numpy array
        RA of coordinate 1 in degrees
    dec1 : float or numpy array
        Dec of coordinate 1 in degrees
    ra2 : float
        RA of coordinate 2 in degrees
    dec2 : float
        Dec of coordinate 2 in degrees

    Returns
    -------
    separation : astropy Angle or numpy array
        Angular separation in degrees

    """
    from astropy.coordinates import SkyCoord
    import astropy.units as u

    coord1 = SkyCoord(ra1, dec1, unit=(u.degree, u.degree), frame='fk5')
    coord2 = SkyCoord(ra2, dec2, unit=(u.degree, u.degree), frame='fk5')

    return coord1.separation(coord2) 

def test_against_pyephem():
    """Check that Astropy gives consistent results with one PyEphem example.

    PyEphem: http://rhodesmill.org/pyephem/

    See example input and output here:
    https://gist.github.com/zonca/1672906
    https://github.com/phn/pytpm/issues/2#issuecomment-3698679
    """
    obstime = Time('2011-09-18 08:50:00')
    location = EarthLocation(lon=Angle('-109d24m53.1s'),
                             lat=Angle('33d41m46.0s'),
                             height=30000. * u.m)
    # We are using the default pressure and temperature in PyEphem
    # relative_humidity = ?
    # obswl = ?
    altaz_frame = AltAz(obstime=obstime, location=location,
                        temperature=15 * u.deg_C, pressure=1.010 * u.bar)

    altaz = SkyCoord('6.8927d -60.7665d', frame=altaz_frame)
    radec_actual = altaz.transform_to('icrs')

    radec_expected = SkyCoord('196.497518d -4.569323d', frame='icrs')  # EPHEM
    # radec_expected = SkyCoord('196.496220d -4.569390d', frame='icrs')  # HORIZON
    distance = radec_actual.separation(radec_expected).to('arcsec')
    # TODO: why is this difference so large?
    # It currently is: 31.45187984720655 arcsec
    assert distance < 1e3 * u.arcsec

    # Add assert on current Astropy result so that we notice if something changes
    radec_expected = SkyCoord('196.495372d -4.560694d', frame='icrs')
    distance = radec_actual.separation(radec_expected).to('arcsec')
    # Current value: 0.0031402822944751997 arcsec
    assert distance < 1 * u.arcsec 

def test_fk5_equinox_and_epoch_j2000_0_to_topocentric_observed():
    """
    http://phn.github.io/pytpm/conversions.html#fk5-equinox-and-epoch-j2000-0-to-topocentric-observed
    """
    # Observatory position for `kpno` from here:
    # http://idlastro.gsfc.nasa.gov/ftp/pro/astro/observatory.pro
    location = EarthLocation(lon=Angle('-111.598333d'),
                             lat=Angle('31.956389d'),
                             height=2093.093 * u.m)  # TODO: height correct?

    obstime = Time('2010-01-01 12:00:00')
    # relative_humidity = ?
    # obswl = ?
    altaz_frame = AltAz(obstime=obstime, location=location,
                        temperature=0 * u.deg_C, pressure=0.781 * u.bar)

    radec = SkyCoord('12h22m54.899s 15d49m20.57s', frame='fk5')

    altaz_actual = radec.transform_to(altaz_frame)

    altaz_expected = SkyCoord('264d55m06s 37d54m41s', frame='altaz')
    # altaz_expected = SkyCoord('343.586827647d 15.7683070508d', frame='altaz')
    # altaz_expected = SkyCoord('133.498195532d 22.0162383595d', frame='altaz')
    distance = altaz_actual.separation(altaz_expected)
    # print(altaz_actual)
    # print(altaz_expected)
    # print(distance)
    """TODO: Current output is completely incorrect ... xfailing this test for now.

    <SkyCoord (AltAz: obstime=2010-01-01 12:00:00.000, location=(-1994497.7199061865, -5037954.447348028, 3357437.2294832403) m, pressure=781.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron):00:00.000, location=(-1994497.7199061865, -5037954.447348028, 3357437.2294832403) m, pressure=781.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron): az=133.4869896371561 deg, alt=67.97857990957701 deg>
    <SkyCoord (AltAz: obstime=None, location=None, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0, obswl=1.0 micron): az=264.91833333333335 deg, alt=37.91138888888889 deg>
    68d02m45.732s
    """

    assert distance < 1 * u.arcsec 

def test_angsep():
    """
    Tests that the angular separation object also behaves correctly.
    """
    from astropy.coordinates.angle_utilities import angular_separation

    # check it both works with floats in radians, Quantities, or Angles
    for conv in (np.deg2rad,
                 lambda x: u.Quantity(x, "deg"),
                 lambda x: Angle(x, "deg")):
        for (lon1, lat1, lon2, lat2), corrsep in zip(coords, correct_seps):
            angsep = angular_separation(conv(lon1), conv(lat1),
                                        conv(lon2), conv(lat2))
            assert np.fabs(angsep - conv(corrsep)) < conv(correctness_margin) 

def test_Gaussian2DRotation():
    amplitude = 42
    x_mean, y_mean = 0, 0
    x_stddev, y_stddev = 2, 3
    theta = Angle(10, 'deg')
    pars = dict(amplitude=amplitude, x_mean=x_mean, y_mean=y_mean,
                x_stddev=x_stddev, y_stddev=y_stddev)
    rotation = models.Rotation2D(angle=theta.degree)
    point1 = (x_mean + 2 * x_stddev, y_mean + 2 * y_stddev)
    point2 = rotation(*point1)
    g1 = models.Gaussian2D(theta=0, **pars)
    g2 = models.Gaussian2D(theta=theta.radian, **pars)
    value1 = g1(*point1)
    value2 = g2(*point2)
    assert_allclose(value1, value2) 

def test_icrs_basic(tmpdir):
    wrap_angle = Angle(1.5, unit=units.rad)
    ra = Longitude(25, unit=units.deg, wrap_angle=wrap_angle)
    dec = Latitude(45, unit=units.deg)

    tree = {'coord': ICRS(ra=ra, dec=dec)}

    assert_roundtrip_tree(tree, tmpdir) 

def test_angle(tmpdir):
    tree = {'angle': Angle(100, u.deg)}
    assert_roundtrip_tree(tree, tmpdir) 

def from_tree(cls, node, ctx):
        return Angle(super().from_tree(node, ctx)) 

def get_lst_for_time(jd_array, latitude, longitude, altitude):
    """
    Get the lsts for a set of jd times at an earth location.

    Parameters
    ----------
    jd_array : ndarray of float
        JD times to get lsts for.
    latitude : float
        Latitude of location to get lst for in degrees.
    longitude : float
        Longitude of location to get lst for in degrees.
    altitude : float
        Altitude of location to get lst for in meters.

    Returns
    -------
    ndarray of float
        LSTs in radians corresponding to the jd_array.

    """
    lst_array = np.zeros_like(jd_array)
    jd, reverse_inds = np.unique(jd_array, return_inverse=True)
    times = Time(
        jd,
        format="jd",
        location=(Angle(longitude, unit="deg"), Angle(latitude, unit="deg")),
    )
    if iers.conf.auto_max_age is None:  # pragma: no cover
        delta, status = times.get_delta_ut1_utc(return_status=True)
        if np.any(
            np.isin(status, (iers.TIME_BEFORE_IERS_RANGE, iers.TIME_BEYOND_IERS_RANGE))
        ):
            warnings.warn(
                "time is out of IERS range, setting delta ut1 utc to "
                "extrapolated value"
            )
            times.delta_ut1_utc = delta
    lst_array = times.sidereal_time("apparent").radian[reverse_inds]

    return lst_array 

def set_vector_rho_phi(self, rho, phi, c=None, **kwargs):
        """sets the vector field using R, Phi for each telescope

        Parameters
        ----------
        rho: float or array[float]
            vector magnitude for each telescope
        phi: array[Angle]
            vector angle for each telescope
        c: color or list of colors
            vector color for each telescope (or one for all)
        """
        phi = Angle(phi).rad
        uu, vv = polar_to_cart(rho, phi)
        self.set_vector_uv(uu, vv, c=c, **kwargs) 

def hours_min_to_decimal(val):
    hours, minutes = val.split(':')
    angle = Angle('{0}h{1}m'.format(hours, minutes))
    return angle.to(u.degree).value 

def deg_to_sexigesimal(value, fmt):
    """
    Displays a degree coordinate value in sexigesimal, given a format of hms or dms.
    """
    a = Angle(value, unit=u.degree)
    if fmt == 'hms':
        return '{0:02.0f}:{1:02.0f}:{2:05.3f}'.format(a.hms.h, a.hms.m, a.hms.s)
    elif fmt == 'dms':
        rep = a.signed_dms
        sign = '-' if rep.sign < 0 else '+'
        return '{0}{1:02.0f}:{2:02.0f}:{3:05.3f}'.format(sign, rep.d, rep.m, rep.s)
    else:
        return 'fmt must be "hms" or "dms"' 

def to_python(self, value):
        try:
            a = float(value)
            return a
        except ValueError:
            try:
                if self.c_type == 'ra':
                    a = Angle(value, unit=u.hourangle)
                else:
                    a = Angle(value, unit=u.degree)
                return a.to(u.degree).value
            except Exception:
                raise ValidationError('Invalid format. Please use sexigesimal or degrees')