Python coordinates.EarthLocation类代码示例

def check_builtin_matches_remote(download_url=True):
    """
    This function checks that the builtin sites registry is consistent with the
    remote registry (or a registry at some other location).

    Note that current this is *not* run by the testing suite (because it
    doesn't start with "test", and is instead meant to be used as a check
    before merging changes in astropy-data)
    """
    builtin_registry = EarthLocation._get_site_registry(force_builtin=True)
    dl_registry = EarthLocation._get_site_registry(force_download=download_url)

    in_dl = {}
    matches = {}
    for name in builtin_registry.names:
        in_dl[name] = name in dl_registry
        if in_dl[name]:
            matches[name] = quantity_allclose(builtin_registry[name], dl_registry[name])
        else:
            matches[name] = False

    if not all(matches.values()):
        # this makes sure we actually see which don't match
        print("In builtin registry but not in download:")
        for name in in_dl:
            if not in_dl[name]:
                print('    ', name)
        print("In both but not the same value:")
        for name in matches:
            if not matches[name] and in_dl[name]:
                print('    ', name, 'builtin:', builtin_registry[name], 'download:', dl_registry[name])
        assert False, "Builtin and download registry aren't consistent - failures printed to stdout"

def test_itrs_vals_5133():
    time = Time('2010-1-1')
    el = EarthLocation.from_geodetic(lon=20*u.deg, lat=45*u.deg, height=0*u.km)

    lons = [20, 30, 20]*u.deg
    lats = [44, 45, 45]*u.deg
    alts = [0, 0, 10]*u.km
    coos = [EarthLocation.from_geodetic(lon, lat, height=alt).get_itrs(time)
            for lon, lat, alt in zip(lons, lats, alts)]

    aaf = AltAz(obstime=time, location=el)
    aacs = [coo.transform_to(aaf) for coo in coos]

    assert all([coo.isscalar for coo in aacs])

    # the ~1 arcsec tolerance is b/c aberration makes it not exact
    assert_quantity_allclose(aacs[0].az, 180*u.deg, atol=1*u.arcsec)
    assert aacs[0].alt < 0*u.deg
    assert aacs[0].distance > 50*u.km

    # it should *not* actually be 90 degrees, b/c constant latitude is not
    # straight east anywhere except the equator... but should be close-ish
    assert_quantity_allclose(aacs[1].az, 90*u.deg, atol=5*u.deg)
    assert aacs[1].alt < 0*u.deg
    assert aacs[1].distance > 50*u.km

    assert_quantity_allclose(aacs[2].alt, 90*u.deg, atol=1*u.arcsec)
    assert_quantity_allclose(aacs[2].distance, 10*u.km)

def airmass_plots(KPNO=False,ING=False,MLO=False):

    observer_site = Observer.at_site("Kitt Peak", timezone="US/Mountain")

    if KPNO:
        print('plotting airmass curves for Kitt Peak')
        observing_location = EarthLocation.of_site('Kitt Peak')
        observer_site = Observer.at_site("Kitt Peak", timezone="US/Mountain")
        start_time = Time('2017-03-12 01:00:00') # UTC time, so 1:00 UTC = 6 pm AZ mountain time
        end_time = Time('2017-03-12 14:00:00')

    elif MLO:
        print('plotting airmass curves for MLO')
        observing_location = EarthLocation.of_site(u'Palomar')
        observer_site = Observer.at_site("Palomar", timezone="US/Pacific")
        # for run starting 2019-Apr-04 at MLO
        start_time = Time('2019-04-03 01:00:00') # need to enter UTC time, MLO UTC+6?
        end_time = Time('2019-04-03 14:00:00')
        
    elif ING:
        print('plotting airmass curves for INT')
        observing_location = EarthLocation.of_site(u'Roque de los Muchachos')
        observer_site = Observer.at_site("Roque de los Muchachos", timezone="GMT")
        # for run starting 2019-Feb-04 at INT
        start_time = Time('2019-02-04 19:00:00') # INT is on UTC
        end_time = Time('2019-02-05 07:00:00')

    #observing_time = Time('2017-05-19 07:00')  # 1am UTC=6pm AZ mountain time
    #observing_time = Time('2018-03-12 07:00')  # 1am UTC=6pm AZ mountain time
    #aa = AltAz(location=observing_location, obstime=observing_time)


    #for i in range(len(pointing_ra)):



    delta_t = end_time - start_time
    observing_time = start_time + delta_t*np.linspace(0, 1, 75)
    nplots = int(sum(obs_mass_flag)/8.)
    print(nplots)
    for j in range(nplots):
        plt.figure()
        legend_list = []
        for i in range(8):
            pointing_center = coords.SkyCoord(pointing_ra[8*j+i]*u.deg, pointing_dec[8*j+i]*u.deg, frame='icrs')
            if i == 3:
                plot_airmass(pointing_center,observer_site,observing_time,brightness_shading=True)
            else:
                plot_airmass(pointing_center,observer_site,observing_time)
            legend_list.append('Pointing %02d'%(8*j+i+1))
    
        plt.legend(legend_list)
        #plt.ylim(0.9,2.5)
        plt.gca().invert_yaxis()
        plt.subplots_adjust(bottom=.15)
        #plt.axvline(x=7*u.hour,ls='--',color='k')
        plt.axhline(y=2,ls='--',color='k')
        plt.savefig(outfile_prefix+'airmass-%02d.png'%(j+1))

def test_regression_simple_5133():
    t = Time('J2010')
    obj = EarthLocation(-1*u.deg, 52*u.deg, height=[100., 0.]*u.km)
    home = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km)
    aa = obj.get_itrs(t).transform_to(AltAz(obstime=t, location=home))

    # az is more-or-less undefined for straight up or down
    assert_quantity_allclose(aa.alt, [90, -90]*u.deg, rtol=1e-5)
    assert_quantity_allclose(aa.distance, [90, 10]*u.km)

def test_EarthLocation_state_online():
    EarthLocation._site_registry = None
    EarthLocation._get_site_registry(force_download=True)
    assert EarthLocation._site_registry is not None

    oldreg = EarthLocation._site_registry
    newreg = EarthLocation._get_site_registry()
    assert oldreg is newreg
    newreg = EarthLocation._get_site_registry(force_download=True)
    assert oldreg is not newreg

def test_is_night():
    lco = Observer(location=EarthLocation.of_site('lco'))  # Las Campanas
    aao = Observer(location=EarthLocation.of_site('aao'))  # Sydney, Australia
    vbo = Observer(location=EarthLocation.of_site('vbo'))  # India

    time1 = Time('2015-07-28 17:00:00')
    nights1 = [observer.is_night(time1) for observer in [lco, aao, vbo]]
    assert np.all(nights1 == [False, True, True])

    time2 = Time('2015-07-28 02:00:00')
    nights2 = [observer.is_night(time2) for observer in [lco, aao, vbo]]
    assert np.all(nights2 == [True, False, False])

    def find_parallax(self, date):
        '''Find the maximum parallax of self.planet on the given date
           from self.observer's location -- in other words, the difference
           in Mars' position between the observer's position and an
           observer at the same latitude but opposite longitude:
           this tells you you how much difference you would see from
           your position if Mars didn't move between your sunrise and sunset.
        '''
        # To calculate from a point on the equator, set lat to 0.
        observer_loc = EarthLocation.from_geodetic(self.location.lon,
                                                   self.location.lat,
                                                   self.location.height)

        # Specify the anti-point.
        # This isn't really an antipode unless lat == 0.
        antipode_loc = EarthLocation.from_geodetic(-observer.lon,
                                                   observer.lat,
                                                   observer.height)

        # XXX Oops, astropy doesn't offer next_rising etc.
        # so we'll need a function to find that before this
        # function can be implemented since it only works
        # when the planet is on the horizon so both the observer
        # and the anti-observer can see it.
        risetime = find_next_rising(planetname, date)

        obs_planet = get_body(planetname, risetime, observer_loc)
        ant_planet = get_body(planetname, risetime, antipode_loc)

        # First, calculate it the straightforward way using the arctan:
        print()
        mars_dist_miles = mars.distance.km / 1.609344
        print("Miles to Mars:", mars_dist_miles)
        earth_mean_radius = 3958.8    # in miles
        half_dist = earth_mean_radius * math.cos(observer_loc.lat)
        print("Distance between observers:", 2. * half_dist)
        par = 2. * math.atan(half_dist / mars_dist_miles) \
              * 180. / math.pi * 3600.
        print("Calculated parallax (arcsec):", par)

        # See what astropy calculates as the difference between observations:
        print()
        print("parallax on %s: RA %f, dec %f" % (antipode.date,
                                                 obs_planet.ra - ant_planet.ra,
                                             obs_planet.dec - ant_planet.dec))
        total_par = (math.sqrt((obs_planet.ra.radians
                                - ant_planet.ra.radians)**2
                               + (obs_planet.dec.radians
                                - ant_planet.dec.radians)**2)
                     * 180. * 3600. / math.pi)
        print("Total parallax (sum of squares): %f arcseconds" % total_par)
        print()

def test_regression_6697():
    """
    Test for regression of a bug in get_gcrs_posvel that introduced errors at the 1m/s level.

    Comparison data is derived from calculation in PINT
    https://github.com/nanograv/PINT/blob/master/pint/erfautils.py
    """
    pint_vels = CartesianRepresentation(*(348.63632871, -212.31704928, -0.60154936), unit=u.m/u.s)
    location = EarthLocation(*(5327448.9957829, -1718665.73869569,  3051566.90295403), unit=u.m)
    t = Time(2458036.161966612, format='jd', scale='utc')
    obsgeopos, obsgeovel = location.get_gcrs_posvel(t)
    delta = (obsgeovel-pint_vels).norm()
    assert delta < 1*u.cm/u.s

def test_EarthLocation_basic():
    greenwichel = EarthLocation.of_site('greenwich')
    lon, lat, el = greenwichel.to_geodetic()
    assert_quantity_allclose(lon, Longitude('0:0:0', unit=u.deg),
                             atol=10*u.arcsec)
    assert_quantity_allclose(lat, Latitude('51:28:40', unit=u.deg),
                             atol=1*u.arcsec)
    assert_quantity_allclose(el, 46*u.m, atol=1*u.m)

    names = EarthLocation.get_site_names()
    assert 'greenwich' in names
    assert 'example_site' in names

    with pytest.raises(KeyError) as exc:
        EarthLocation.of_site('nonexistent site')
    assert exc.value.args[0] == "Site 'nonexistent site' not in database. Use EarthLocation.get_site_names to see available sites."

def test_vega_sirius_transit_seattle():
    """
    Check that time of transit of Vega for an observer in Seattle is
    consistent with PyEphem results (for no atmosphere/pressure=0)
    """
    lat = "47d36m34.92s"
    lon = "122d19m59.16s"
    elevation = 0.0 * u.m
    pressure = 0 * u.bar
    location = EarthLocation.from_geodetic(lon, lat, elevation)
    time = Time("1990-01-01 12:00:00")
    vega = SkyCoord(279.23473479 * u.degree, 38.78368896 * u.degree)
    sirius = SkyCoord(101.28715533 * u.degree, -16.71611586 * u.degree)

    obs = Observer(location=location, pressure=pressure)
    astroplan_vega_transit = obs.target_meridian_transit_time(time, vega, which="next").datetime
    astroplan_sirius_transit = obs.target_meridian_transit_time(time, sirius, which="next").datetime

    astroplan_vector_transit = obs.target_meridian_transit_time(time, [vega, sirius], which="next").datetime

    # Run print_pyephem_vega_sirius_transit() to compute analogous
    # result from PyEphem:
    pyephem_vega_transit = datetime.datetime(1990, 1, 2, 3, 41, 9, 244067)
    pyephem_sirius_transit = datetime.datetime(1990, 1, 1, 15, 51, 15, 135167)

    # Typical difference in this example between PyEphem and astroplan
    # with an atmosphere is <2 min
    threshold_minutes = 8
    assert abs(pyephem_vega_transit - astroplan_vega_transit) < datetime.timedelta(minutes=threshold_minutes)
    assert abs(pyephem_sirius_transit - astroplan_sirius_transit) < datetime.timedelta(minutes=threshold_minutes)

    # Now check vectorized solutions against scalar:
    assert astroplan_vector_transit[0] == astroplan_vega_transit
    assert astroplan_vector_transit[1] == astroplan_sirius_transit

def test_sunrise_sunset_equator_civil_twilight():
    """
    Check that time of sunrise/set for an observer on the equator is
    consistent with PyEphem results (for no atmosphere/pressure=0)
    """
    lat = "00:00:00"
    lon = "00:00:00"
    elevation = 0.0 * u.m
    pressure = 0 * u.bar
    location = EarthLocation.from_geodetic(lon, lat, elevation)
    time = Time("2000-01-01 12:00:00")
    obs = Observer(location=location, pressure=pressure)
    # Manually impose horizon equivalent to civil twilight
    horizon = -6 * u.degree
    astroplan_next_sunrise = obs.sun_rise_time(time, which="next", horizon=horizon).datetime
    astroplan_next_sunset = obs.sun_set_time(time, which="next", horizon=horizon).datetime

    astroplan_prev_sunrise = obs.sun_rise_time(time, which="previous", horizon=horizon).datetime
    astroplan_prev_sunset = obs.sun_set_time(time, which="previous", horizon=horizon).datetime

    # Run print_pyephem_sunrise_sunset_equator_civil_twilight() to compute
    # analogous result from PyEphem:
    pyephem_next_rise = datetime.datetime(2000, 1, 2, 5, 37, 34, 83328)
    pyephem_next_set = datetime.datetime(2000, 1, 1, 18, 29, 29, 195908)
    pyephem_prev_rise = datetime.datetime(2000, 1, 1, 5, 37, 4, 701708)
    pyephem_prev_set = datetime.datetime(1999, 12, 31, 18, 29, 1, 530987)

    threshold_minutes = 8
    assert abs(pyephem_next_rise - astroplan_next_sunrise) < datetime.timedelta(minutes=threshold_minutes)
    assert abs(pyephem_next_set - astroplan_next_sunset) < datetime.timedelta(minutes=threshold_minutes)
    assert abs(pyephem_prev_rise - astroplan_prev_sunrise) < datetime.timedelta(minutes=threshold_minutes)
    assert abs(pyephem_prev_set - astroplan_prev_sunset) < datetime.timedelta(minutes=threshold_minutes)

def test_exceptions():
    lat = "00:00:00"
    lon = "00:00:00"
    elevation = 0.0 * u.m
    location = EarthLocation.from_geodetic(lon, lat, elevation)
    time = Time("2000-01-01 12:00:00")
    vega_coords = SkyCoord("18h36m56.33635s", "+38d47m01.2802s")

    obs = Observer(location=location)

    with pytest.raises(ValueError):
        obs.target_rise_time(time, vega_coords, which="oops").datetime

    with pytest.raises(ValueError):
        obs.target_set_time(time, vega_coords, which="oops").datetime

    with pytest.raises(ValueError):
        obs.target_meridian_transit_time(time, vega_coords, which="oops").datetime

    with pytest.raises(ValueError):
        obs.target_meridian_antitransit_time(time, vega_coords, which="oops").datetime

    with pytest.raises(TypeError):
        FixedTarget(["00:00:00", "00:00:00"], name="VE")

    with pytest.raises(TypeError):
        Observer(location="Greenwich")

    with pytest.raises(TypeError):
        Observer(location=EarthLocation(0, 0, 0), timezone=-6)

    def from_itrf(ant_itrf_xyz, ant_labels):
        """
        Instantiate telescope model from ITRF co-ordinates of antennae

        (aka ECEF-coords, i.e. Earth-Centered-Earth-Fixed reference frame.)

        Takes care of calculating central Latitude and Longitude from the mean
        antenna position, and converts the antenna positions into local-XYZ
        frame.

        Args:
            ant_itrf_xyz (numpy.ndarray): Array co-ordinatates in the ITRF frame
            ant_labels (list[str]): Antennae labels
        Returns:
            Telescope: A telescope class with the given array co-ords.
        """
        mean_posn = np.mean(ant_itrf_xyz, axis=0)
        centre = EarthLocation.from_geocentric(mean_posn[0],
                                               mean_posn[1],
                                               mean_posn[2],
                                               unit=u.m,
                                               )
        lon, lat, height = centre.to_geodetic()

        mean_subbed_itrf = ant_itrf_xyz - mean_posn

        rotation = z_rotation_matrix(lon)
        ant_local_xyz = np.dot(rotation, mean_subbed_itrf.T).T

        return Telescope(
            centre=centre,
            ant_labels=ant_labels,
            ant_itrf_xyz=ant_itrf_xyz,
            ant_local_xyz=ant_local_xyz,
        )

def test_moon_rise_set():
    pyephem_next_rise = datetime.datetime(2017, 10, 7, 23, 50, 24, 407018)
    pyephem_next_set = datetime.datetime(2017, 10, 7, 12, 30, 30, 787116)
    pyephem_prev_rise = datetime.datetime(2017, 10, 6, 23, 13, 43, 644455)
    pyephem_prev_set = datetime.datetime(2017, 10, 6, 11, 20, 9, 340009)

    time = Time('2017-10-07 12:00:00')
    lat = '42:00:00'
    lon = '-70:00:00'
    elevation = 0.0 * u.m
    pressure = 0 * u.bar
    location = EarthLocation.from_geodetic(lon, lat, elevation)

    obs = Observer(location=location)

    astroplan_next_rise = obs.moon_rise_time(time, which='next')
    astroplan_next_set = obs.moon_set_time(time, which='next')
    astroplan_prev_rise = obs.moon_rise_time(time, which='previous')
    astroplan_prev_set = obs.moon_set_time(time, which='previous')

    threshold_minutes = 2
    assert (abs(pyephem_next_rise - astroplan_next_rise.datetime) <
            datetime.timedelta(minutes=threshold_minutes))
    assert (abs(pyephem_next_set - astroplan_next_set.datetime) <
            datetime.timedelta(minutes=threshold_minutes))
    assert (abs(pyephem_prev_rise - astroplan_prev_rise.datetime) <
            datetime.timedelta(minutes=threshold_minutes))
    assert (abs(pyephem_prev_set - astroplan_prev_set.datetime) <
            datetime.timedelta(minutes=threshold_minutes))

def eci2el(x,y,z,dt):
    """
    Convert Earth-Centered Inertial (ECI) cartesian coordinates to ITRS for astropy EarthLocation object.

    Inputs :
    x = ECI X-coordinate 
    y = ECI Y-coordinate 
    z = ECI Z-coordinate 
    dt = UTC time (datetime object)
    """

    from astropy.coordinates import GCRS, ITRS, EarthLocation, CartesianRepresentation
    import astropy.units as u
    
    # convert datetime object to astropy time object
    tt=Time(dt,format='datetime')

    # Read the coordinates in the Geocentric Celestial Reference System
    gcrs = GCRS(CartesianRepresentation(x=x, y=y,z=z), obstime=tt)

    # Convert it to an Earth-fixed frame
    itrs = gcrs.transform_to(ITRS(obstime=tt))

    el = EarthLocation.from_geocentric(itrs.x, itrs.y, itrs.z) 

    return el