Python astropy.coordinates.SkyCoord() Examples

def radec_obs_vec_mpc(inds, mpc_object_data):
    """Compute vector of observed ra,dec values for MPC tracking data.

       Args:
           inds (int array): line numbers of data in file
           mpc_object_data (ndarray): MPC observation data for object

       Returns:
           rov (1xlen(inds) array): vector of ra/dec observed values
    """
    rov = np.zeros((2*len(inds)))
    for i in range(0,len(inds)):
        indm1 = inds[i]-1
        # extract observations data
        timeobs = Time( datetime(mpc_object_data['yr'][indm1],
                                 mpc_object_data['month'][indm1],
                                 mpc_object_data['day'][indm1]) + timedelta(days=mpc_object_data['utc'][indm1]) )
        obs_t_ra_dec = SkyCoord(mpc_object_data['radec'][indm1], unit=(uts.hourangle, uts.deg), obstime=timeobs)
        rov[2*i-2], rov[2*i-1] = obs_t_ra_dec.ra.rad, obs_t_ra_dec.dec.rad
    return rov

# compute residuals vector for ra/dec observations with pre-computed observed radec values vector 

def query(self, coords, **kwargs):
        """
        Returns E(B-V) (or a different Planck dust inference, depending on how
        the class was intialized) at the specified location(s) on the sky.

        Args:
            coords (:obj:`astropy.coordinates.SkyCoord`): The coordinates to query.

        Returns:
            A float array of the selected Planck component, at the given
            coordinates. The shape of the output is the same as the shape of the
            coordinates stored by ``coords``. If extragalactic E(B-V), tau_353
            or radiance was chosen, then the output has units of magnitudes of
            E(B-V). If the selected Planck component is temperature (or
            temperature error), then an :obj:`astropy.Quantity` is returned, with
            units of Kelvin. If beta (or beta error) was chosen, then the output
            is unitless.
        """
        return self._scale * super(PlanckQuery, self).query(coords, **kwargs) 

def deserialize_skycoord(d):
    """
    Deserializes a JSONified :obj:`astropy.coordinates.SkyCoord`.

    Args:
        d (:obj:`dict`): A dictionary representation of a :obj:`SkyCoord` object.

    Returns:
        A :obj:`SkyCoord` object.
    """
    if 'distance' in d:
        args = (d['lon'], d['lat'], d['distance'])
    else:
        args = (d['lon'], d['lat'])

    return coords.SkyCoord(
        *args,
        frame=d['frame'],
        representation='spherical') 

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

        for reps in range(10):
            # 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._sfd(c)

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

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_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_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_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_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_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_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_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 inverse_method(self,N,d):
        
        t = np.linspace(1e-3,0.999,N)
        f = np.log( t / (1 - t) )
        f = f/f[0]
        
        psi= np.pi*f
        cosPsi = np.cos(psi)
        sinTheta = ( np.abs(cosPsi) + (1-np.abs(cosPsi))*np.random.rand(len(cosPsi)))
        
        theta = np.arcsin(sinTheta)
        theta = np.pi-theta + (2*theta - np.pi)*np.round(np.random.rand(len(t)))
        cosPhi = cosPsi/sinTheta
        phi = np.arccos(cosPhi)*(-1)**np.round(np.random.rand(len(t)))
        
        coords = SkyCoord(phi*u.rad,(np.pi/2-theta)*u.rad,d*np.ones(len(phi))*u.pc)

        return coords 

def radec_residual_mpc(x, t_ra_dec_datapoint, long, parallax_s, parallax_c):
    """Compute observed minus computed (O-C) residual for a given ra/dec
    datapoint, represented as a SkyCoord object, for MPC observation data.

       Args:
           x (1x6 array): set of Keplerian elements
           t_ra_dec_datapoint (SkyCoord): ra/dec datapoint
           long (float): longitude of observing site
           parallax_s (float): parallax constant S of observing site
           parallax_c (float): parallax constant C of observing site

       Returns:
           (1x2 array): right ascension difference, declination difference
    """
    ra_comp, dec_comp = rhovec2radec(long, parallax_s, parallax_c, t_ra_dec_datapoint.obstime,
                                     x[0], x[1], x[2], x[3], x[4], x[5])
    ra_obs, dec_obs = t_ra_dec_datapoint.ra.rad, t_ra_dec_datapoint.dec.rad
    #"unsigned" distance between points in torus
    diff_ra = angle_diff_rad(ra_obs, ra_comp)
    diff_dec = angle_diff_rad(dec_obs, dec_comp)
    return np.array((diff_ra,diff_dec)) 

def get_observer_pos_wrt_earth(sat_observatories_data, obs_radec, site_codes):
    """Compute position of observer at Earth's surface, with respect
    to the Earth, in equatorial frame, during 3 distinct instants.

       Args:
           sat_observatories_data (string): path to file containing COSPAR satellite tracking stations data.
           obs_radec (1x3 SkyCoord array): three rad/dec observations
           site_codes (1x3 int array): COSPAR codes of observation sites

       Returns:
           R (1x3 array): cartesian position vectors (observer wrt Earth)
    """
    R = np.array((np.zeros((3,)),np.zeros((3,)),np.zeros((3,))))
    # load MPC observatory data
    obsite1 = get_station_data(site_codes[0], sat_observatories_data)
    obsite2 = get_station_data(site_codes[1], sat_observatories_data)
    obsite3 = get_station_data(site_codes[2], sat_observatories_data)

    R[0] = observerpos_sat(obsite1['Latitude'], obsite1['Longitude'], obsite1['Elev'], obs_radec[0].obstime)
    R[1] = observerpos_sat(obsite2['Latitude'], obsite2['Longitude'], obsite2['Elev'], obs_radec[1].obstime)
    R[2] = observerpos_sat(obsite3['Latitude'], obsite3['Longitude'], obsite3['Elev'], obs_radec[2].obstime)

    return R 

def _setup_skyonly(self,min_nobs,max_nobs_std,max_plxerr,max_scd,min_lat):
        self._tgas_sid= (self._full_tgas['source_id']/2**(35.\
                               +2*(12.-numpy.log2(_BASE_NSIDE)))).astype('int')
        self._tgas_sid_skyonlyindx= (self._full_jk > 0.)\
            *(self._full_jk < 0.8)\
            *(self._full_twomass['j_mag'] > 6.)\
            *(self._full_twomass['j_mag'] < 10.)
        nstar, e= numpy.histogram(self._tgas_sid[self._tgas_sid_skyonlyindx],
                                  range=[-0.5,_BASE_NPIX-0.5],bins=_BASE_NPIX)
        self._nstar_tgas_skyonly= nstar
        self._nobs_tgas_skyonly= self._compute_mean_quantity_tgas(\
            'astrometric_n_good_obs_al',lambda x: x/9.)
        self._nobsstd_tgas_skyonly= numpy.sqrt(\
            self._compute_mean_quantity_tgas(\
                'astrometric_n_good_obs_al',lambda x: (x/9.)**2.)
            -self._nobs_tgas_skyonly**2.)
        self._scank4_tgas_skyonly= self._compute_mean_quantity_tgas(\
            'scan_direction_strength_k4')
        self._plxerr_tgas_skyonly= self._compute_mean_quantity_tgas(\
            'parallax_error')
        tmp_decs, ras= healpy.pix2ang(_BASE_NSIDE,numpy.arange(_BASE_NPIX),
                                      nest=True)
        coos= apco.SkyCoord(ras,numpy.pi/2.-tmp_decs,unit="rad")
        coos= coos.transform_to(apco.GeocentricTrueEcliptic)
        self._eclat_skyonly= coos.lat.to('deg').value
        self._exclude_mask_skyonly= \
            (self._nobs_tgas_skyonly < min_nobs)\
            +(self._nobsstd_tgas_skyonly > max_nobs_std)\
            +(numpy.fabs(self._eclat_skyonly) < min_lat)\
            +(self._plxerr_tgas_skyonly > max_plxerr)\
            +(self._scank4_tgas_skyonly > max_scd)
        return None 

def get_observations_data(mpc_object_data, inds):
    """Extract three ra/dec observations from MPC observation data file.

       Args:
           mpc_object_data (string): file path to MPC observation data of object
           inds (int array): indices of requested data

       Returns:
           obs_radec (1x3 SkyCoord array): ra/dec observation data
           obs_t (1x3 Time array): time observation data
           site_codes (1x3 int array): corresponding codes of observation sites
    """
    # construct SkyCoord 3-element array with observational information
    timeobs = np.zeros((3,), dtype=Time)
    obs_radec = np.zeros((3,), dtype=SkyCoord)
    obs_t = np.zeros((3,))

    timeobs[0] = Time( datetime(mpc_object_data['yr'][inds[0]],
                                mpc_object_data['month'][inds[0]],
                                mpc_object_data['day'][inds[0]]) + timedelta(days=mpc_object_data['utc'][inds[0]]) )
    timeobs[1] = Time( datetime(mpc_object_data['yr'][inds[1]],
                                mpc_object_data['month'][inds[1]],
                                mpc_object_data['day'][inds[1]]) + timedelta(days=mpc_object_data['utc'][inds[1]]) )
    timeobs[2] = Time( datetime(mpc_object_data['yr'][inds[2]],
                                mpc_object_data['month'][inds[2]],
                                mpc_object_data['day'][inds[2]]) + timedelta(days=mpc_object_data['utc'][inds[2]]) )

    obs_radec[0] = SkyCoord(mpc_object_data['radec'][inds[0]], unit=(uts.hourangle, uts.deg), obstime=timeobs[0])
    obs_radec[1] = SkyCoord(mpc_object_data['radec'][inds[1]], unit=(uts.hourangle, uts.deg), obstime=timeobs[1])
    obs_radec[2] = SkyCoord(mpc_object_data['radec'][inds[2]], unit=(uts.hourangle, uts.deg), obstime=timeobs[2])

    # construct vector of observation time (continous variable)
    obs_t[0] = obs_radec[0].obstime.tdb.jd
    obs_t[1] = obs_radec[1].obstime.tdb.jd
    obs_t[2] = obs_radec[2].obstime.tdb.jd

    site_codes = [mpc_object_data['observatory'][inds[0]],
                  mpc_object_data['observatory'][inds[1]],
                  mpc_object_data['observatory'][inds[2]]]

    return obs_radec, obs_t, site_codes 

def coords(self):
        """
        A :class:`~astropy.coordinates.SkyCoord` object.
        """
        if self._frame not in spice_astropy_frame_mapping:
            raise ValueError(f'Current frame "{self._frame}" not in list of '
                             f'known coordinate frames implemented in astropy '
                             f'or sunpy ({spice_astropy_frame_mapping})')
        if self._abcorr.lower() != 'none':
            raise NotImplementedError(
                'Can only convert to astropy coordinates if the aberration '
                'correction is set to "none" '
                f'(currently set to {self._abcorr})')

        frame = spice_astropy_frame_mapping[self._frame][0]

        # Error if sunpy < 2 due to changes in Heliographic coordinates then
        if (frame == suncoords.HeliographicCarrington and
                int(sunpy.__version__[0]) < 2):
            raise NotImplementedError(
                'Converting to Carrington coordinates only works for '
                f'sunpy versions >= 2.0 (found sunpy {sunpy.__version__} '
                'installed)'
            )

        kwargs = spice_astropy_frame_mapping[self._frame][1]
        coords = astrocoords.SkyCoord(
            self.x, self.y, self.z,
            frame=frame, representation_type='cartesian',
            obstime=self.times,
            **kwargs)
        coords.representation_type = frame().default_representation
        return coords 

def boundaries_skycoord(self, healpix_index, step):
        """
        Return the celestial coordinates of the edges of HEALPix pixels

        This returns the celestial coordinates of points along the edge of each
        HEALPIX pixel. The number of points returned for each pixel is ``4 * step``,
        so setting ``step`` to 1 returns just the corners.

        This method requires that a celestial frame was specified when
        initializing HEALPix.  If you don't know or need the celestial frame,
        you can instead use :meth:`~astropy_healpix.HEALPix.boundaries_lonlat`.

        Parameters
        ----------
        healpix_index : `~numpy.ndarray`
            1-D array of HEALPix pixels
        step : int
            The number of steps to take along each edge.

        Returns
        -------
        skycoord : :class:`~astropy.coordinates.SkyCoord`
            The celestial coordinates of the HEALPix pixel boundaries
        """
        if self.frame is None:
            raise NoFrameError("boundaries_skycoord")
        lon, lat = self.boundaries_lonlat(healpix_index, step)
        representation = UnitSphericalRepresentation(lon, lat, copy=False)
        return SkyCoord(self.frame.realize_frame(representation)) 

def cone_search_skycoord(self, skycoord, radius):
        """
        Find all the HEALPix pixels within a given radius of a celestial position.

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

        This method requires that a celestial frame was specified when
        initializing HEALPix.  If you don't know or need the celestial frame,
        you can instead use :meth:`~astropy_healpix.HEALPix.cone_search_lonlat`.

        Parameters
        ----------
        skycoord : :class:`~astropy.coordinates.SkyCoord`
            The celestial coordinates to use for the cone search
        radius : :class:`~astropy.units.Quantity`
            The search radius

        Returns
        -------
        healpix_index : `~numpy.ndarray`
            1-D array with all the matching HEALPix pixel indices.
        """
        if self.frame is None:
            raise NoFrameError("cone_search_skycoord")
        skycoord = skycoord.transform_to(self.frame)
        representation = skycoord.represent_as(UnitSphericalRepresentation)
        lon, lat = representation.lon, representation.lat
        return self.cone_search_lonlat(lon, lat, radius) 

def skycoord_to_healpix(self, skycoord, return_offsets=False):
        """
        Convert celestial coordinates to HEALPix indices (optionally with offsets).

        Note that this method requires that a celestial frame was specified when
        initializing HEALPix. If you don't know or need the celestial frame, you
        can instead use :meth:`~astropy_healpix.HEALPix.lonlat_to_healpix`.

        Parameters
        ----------
        skycoord : :class:`~astropy.coordinates.SkyCoord`
            The celestial coordinates to convert
        return_offsets : bool
            If `True`, the returned values are the HEALPix pixel as well as
            ``dx`` and ``dy``, the fractional positions inside the pixel. If
            `False` (the default), only the HEALPix pixel is returned.

        Returns
        -------
        healpix_index : `~numpy.ndarray`
            1-D array of HEALPix indices
        dx, dy : `~numpy.ndarray`
            1-D arrays of offsets inside the HEALPix pixel in the range [0:1] (0.5
            is the center of the HEALPix pixels). This is returned if
            ``return_offsets`` is `True`.
        """
        if self.frame is None:
            raise NoFrameError("skycoord_to_healpix")
        skycoord = skycoord.transform_to(self.frame)
        representation = skycoord.represent_as(UnitSphericalRepresentation)
        lon, lat = representation.lon, representation.lat
        return self.lonlat_to_healpix(lon, lat, return_offsets=return_offsets) 

def healpix_to_skycoord(self, healpix_index, dx=None, dy=None):
        """
        Convert HEALPix indices (optionally with offsets) to celestial coordinates.

        Note that this method requires that a celestial frame was specified when
        initializing HEALPix. If you don't know or need the celestial frame, you
        can instead use :meth:`~astropy_healpix.HEALPix.healpix_to_lonlat`.

        Parameters
        ----------
        healpix_index : `~numpy.ndarray`
            1-D array of HEALPix indices
        dx, dy : `~numpy.ndarray`, optional
            1-D arrays of offsets inside the HEALPix pixel, which must be in
            the range [0:1] (0.5 is the center of the HEALPix pixels). If not
            specified, the position at the center of the pixel is used.

        Returns
        -------
        coord : :class:`~astropy.coordinates.SkyCoord`
            The resulting celestial coordinates
        """
        if self.frame is None:
            raise NoFrameError("healpix_to_skycoord")
        lon, lat = self.healpix_to_lonlat(healpix_index, dx=dx, dy=dy)
        representation = UnitSphericalRepresentation(lon, lat, copy=False)
        return SkyCoord(self.frame.realize_frame(representation)) 

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_healpix_to_skycoord(self):
        coord = self.pix.healpix_to_skycoord([1, 2, 3])

        assert isinstance(coord, SkyCoord)
        assert isinstance(coord.frame, Galactic)

        # Make sure that the skycoord_to_healpix method converts coordinates
        # to the frame of the HEALPix
        coord = coord.transform_to('fk5')

        index = self.pix.skycoord_to_healpix(coord)

        assert_equal(index, [1, 2, 3])

        coord = self.pix.healpix_to_skycoord([1, 2, 3],
                                             dx=[0.1, 0.2, 0.3],
                                             dy=[0.5, 0.4, 0.7])

        assert isinstance(coord, SkyCoord)
        assert isinstance(coord.frame, Galactic)

        # Make sure that the skycoord_to_healpix method converts coordinates
        # to the frame of the HEALPix
        coord = coord.transform_to('fk5')

        index, dx, dy = self.pix.skycoord_to_healpix(coord, return_offsets=True)

        assert_equal(index, [1, 2, 3])
        assert_allclose(dx, [0.1, 0.2, 0.3])
        assert_allclose(dy, [0.5, 0.4, 0.7]) 

def revise_lists(self, sInds):
        """Replaces Target List catalog attributes with filtered values, 
        and updates the number of target stars.
        
        Args:
            sInds (integer ndarray):
                Integer indices of the stars of interest
        
        """
       
        # cast sInds to array
        sInds = np.array(sInds, ndmin=1, copy=False)
        
        if len(sInds) == 0:
            raise IndexError("Target list filtered to empty.")
        
        for att in self.catalog_atts:
            if att == 'coords':
                ra = self.coords.ra[sInds].to('deg')
                dec = self.coords.dec[sInds].to('deg')
                self.coords = SkyCoord(ra, dec, self.dist.to('pc'))
            else:
                if getattr(self, att).size != 0:
                    setattr(self, att, getattr(self, att)[sInds])
        try:
            self.Completeness.revise_updates(sInds)
        except AttributeError:
            pass
        self.nStars = len(sInds)
        assert self.nStars, "Target list is empty: nStars = %r"%self.nStars 

def gal_to_shape(gal, shape):
    l = np.reshape(gal.l.deg, shape)*units.deg
    b = np.reshape(gal.b.deg, shape)*units.deg

    has_dist = hasattr(gal.distance, 'kpc')
    d = np.reshape(gal.distance.kpc, shape)*units.kpc if has_dist else None

    return coordinates.SkyCoord(l, b, distance=d, frame='galactic') 

def main():
    w,h = (2056,1024)
    l_0 = 0.

    # Create a grid of coordinates
    print('Creating grid of coordinates...')
    l = np.linspace(-180.+l_0, 180.+l_0, 2*w)
    b = np.linspace(-90., 90., 2*h+2)
    b = b[1:-1]
    l,b = np.meshgrid(l, b)

    l += (np.random.random(l.shape) - 0.5) * 360./(2.*w)
    b += (np.random.random(l.shape) - 0.5) * 180./(2.*h)

    coords = SkyCoord(l*u.deg, b*u.deg, frame='galactic')

    # Set up BH query object
    print('Loading BH map...')
    bh = BHQuery()

    print('Querying map...')
    ebv = bh.query(coords)

    # Convert the output array to a PIL image and save
    print('Saving image...')
    img = numpy2pil(ebv[::-1,::-1], 0., 1.5)
    img = img.resize((w,h), resample=PIL.Image.LANCZOS)
    fname = 'bh.png'
    img.save(fname)

    return 0 

def object_hook(self, d):
        if isinstance(d, dict):
            if ('_type' in d):
                if d['_type'] == 'astropy.coordinates.SkyCoord':
                    return deserialize_skycoord(d)
                elif d['_type'] == 'astropy.units.Quantity':
                    return deserialize_quantity(d)
                elif d['_type'] == 'np.ndarray':
                    return deserialize_ndarray(d)
                elif d['_type'] == 'np.dtype':
                    return deserialize_dtype(d)
                elif d['_type'] == 'tuple':
                    return deserialize_tuple(d)
        return d 

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)