How to use the healpy.pix2ang function in healpy

def setUp(self) :
        self.nside=64
        self.lmax=3*self.nside-1
        self.ntemp=5
        self.npix=int(hp.nside2npix(self.nside))
        self.msk=np.ones(self.npix)
        self.mps=np.zeros([3,self.npix])
        self.tmp=np.zeros([self.ntemp,3,self.npix])
        self.beam=np.ones(self.lmax+1)

        th,ph=hp.pix2ang(self.nside,np.arange(self.npix))
        sth=np.sin(th); cth=np.cos(th)
        self.mps[0]=np.sqrt(15./2./np.pi)*sth**2*np.cos(2*ph) #Re(Y_22)
        self.mps[1]=-np.sqrt(15./2./np.pi)*sth**2/4.          #_2Y^E_20 + _2Y^B_30
        self.mps[2]=-np.sqrt(105./2./np.pi)*cth*sth**2/2.
        for i in range(self.ntemp) :
            self.tmp[i][0]=np.sqrt(15./2./np.pi)*sth**2*np.cos(2*ph) #Re(Y_22)
            self.tmp[i][1]=-np.sqrt(15./2./np.pi)*sth**2/4.          #_2Y^E_20 + _2Y^B_30
            self.tmp[i][2]=-np.sqrt(105./2./np.pi)*cth*sth**2/2.

def setUp(self) :
        self.nside=256
        self.th0=np.pi/4
        self.msk=np.zeros(hp.nside2npix(self.nside),dtype=float)
        self.th,ph=hp.pix2ang(self.nside,np.arange(hp.nside2npix(self.nside),dtype=int))
        self.msk[self.th

def AzimuthalRotation(hmap):
    """
    Azimuthal rotation of a healpix map by pi/2 about the z-axis
    """
    npix = len(hmap)
    nside= hp.npix2nside(npix)
    hpxidx = np.arange(npix)
    t2,p2 = hp.pix2ang(nside, hpxidx)

    p = p2 - np.pi/2
    p[p < 0] += 2. * np.pi
    t = t2

    idx = hp.ang2pix(nside, t, p)

    hout = hmap[idx]
    return hout

def __init__(self):
        """ Initialize the Observer object.

        Calls ephem.Observer.__init__ function and adds on gsm
        """
        super(GSMObserver2016, self).__init__()
        self.gsm = GlobalSkyModel2016(freq_unit='MHz')
        self.observed_sky = None

        # Generate mapping from pix <-> angles
        self.gsm.generate(1000)
        self._n_pix  = hp.get_map_size(self.gsm.generated_map_data)
        self._n_side = hp.npix2nside(self._n_pix)
        self._theta, self._phi = hp.pix2ang(self._n_side, np.arange(self._n_pix))

JK= (object-wide default) J-Ks color or an array of samples of the J-Ks color 
           relative= (False) if True, compute the effective volume completeness = effective volume / true volume; computed using the same integration grid, so will be more robust against integration errors (especially due to the finite HEALPix grid for the angular integration). For simple volumes, a more precise effective volume can be computed by using relative=True and multiplying in the correct true volume
           ndists= (101) number of distances to use in the distance integration
           linearDist= (False) if True, integrate in distance rather than distance modulus
           ncpu= (None) if set to an integer, use this many CPUs to compute the effective selection function (only for non-zero extinction)
        OUTPUT
           effective volume
        HISTORY:
           2017-01-18 - Written - Bovy (UofT/CCA)
        """           
        # Pre-compute coordinates for integrand evaluation            
        if not hasattr(self,'_ra_cen_4vol') or \
                (hasattr(self,'_ndists_4vol') and 
                 (ndists != self._ndists_4vol or 
                  linearDist != self._linearDist_4vol)):
            theta,phi= healpy.pix2ang(\
                _BASE_NSIDE,numpy.arange(_BASE_NPIX)\
                    [True^self._tgasSel._exclude_mask_skyonly],nest=True)
            self._ra_cen_4vol= 180./numpy.pi*phi
            self._dec_cen_4vol= 90.-180./numpy.pi*theta
            if linearDist:
                dists= numpy.linspace(0.001,10.,ndists)
                dms= 5.*numpy.log10(dists)+10.
                self._deltadm_4vol= dists[1]-dists[0]
            else:
                dms= numpy.linspace(0.,18.,ndists)
                self._deltadm_4vol= (dms[1]-dms[0])*numpy.log(10.)/5.
            self._dists_4vol= 10.**(0.2*dms-2.)
            self._tiled_dists3_4vol= numpy.tile(\
                self._dists_4vol**(3.-linearDist),(len(self._ra_cen_4vol),1))
            self._tiled_ra_cen_4vol= numpy.tile(self._ra_cen_4vol,
                                                 (len(self._dists_4vol),1)).T

# read probability skymap
     hpx = hp.read_map( sky_map, verbose = False )

     # number of pixels
     npix = len( hpx )

     # nside: resolution for the HEALPix map
     nside = hp.npix2nside( npix )

     # pixel position
     ipix_max = np.argmax( hpx )
     hpx[ ipix_max ]

     # sky coordinate
     theta, phi = hp.pix2ang( nside, ipix_max )

     ra_max = np.rad2deg( phi )

     dec_max = np.rad2deg( 0.5 * np.pi - theta )

     print ' or insert the highest probability pixel located \n at RA =', str(( '% .5f' % ra_max))+'°', 'and Dec =',str(( '% .5f' % dec_max))+'°.'

vals = []
            for i, ax in enumerate(self.axes):
                bins += [pix[1 + i]]
                vals += [ax.pix_to_coord(pix[1 + i])]

            idxs = pix_tuple_to_idx(bins)

            if self.nside.size > 1:
                nside = self.nside[idxs]
            else:
                nside = self.nside

            ipix = np.round(pix[0]).astype(int)
            m = ipix == INVALID_INDEX.int
            ipix[m] = 0
            theta, phi = hp.pix2ang(nside, ipix, nest=self.nest)
            coords = [np.degrees(phi), np.degrees(np.pi / 2.0 - theta)]
            coords = tuple(coords + vals)
            if np.any(m):
                for c in coords:
                    c[m] = INVALID_INDEX.float
        else:
            ipix = np.round(pix[0]).astype(int)
            theta, phi = hp.pix2ang(self.nside, ipix, nest=self.nest)
            coords = (np.degrees(phi), np.degrees(np.pi / 2.0 - theta))

        return coords

max_order = np.max(order)
    max_nside = hp.order2nside(max_order)
    max_ipix = ipix << np.uint64(2 * (max_order - order))
    ipix = ipix.astype(np.int64)
    max_ipix = max_ipix.astype(np.int64)
    true_theta = 0.5 * np.pi - true_dec
    true_phi = true_ra
    true_pix = hp.ang2pix(max_nside, true_theta, true_phi, nest=True)
    # At this point, we could sort the dataset by max_ipix and then do a binary
    # (e.g., np.searchsorted) to find true_pix in max_ipix. However, would be
    # slower than the linear search below because the sort would be N log N.
    i = np.flatnonzero(max_ipix <= true_pix)
    true_idx = i[np.argmax(max_ipix[i])]

    # Find the angular offset between the mode and true locations.
    mode_theta, mode_phi = hp.pix2ang(
        hp.order2nside(order[0]), ipix[0].astype(np.int64), nest=True)
    offset = np.rad2deg(
        angle_distance(true_theta, true_phi, mode_theta, mode_phi))

    # Calculate the cumulative area in deg2 and the cumulative probability.
    area = moc.uniq2pixarea(sky_map['UNIQ'])
    prob = np.cumsum(sky_map['PROBDENSITY'] * area)
    area = np.cumsum(area) * np.square(180 / np.pi)

    # Construct linear interpolants to map between probability and area.
    # This allows us to compute more accurate contour areas and probabilities
    # under the approximation that the pixels have constant probability
    # density.
    prob_padded = np.concatenate(([0], prob))
    area_padded = np.concatenate(([0], area))
    prob_for_area = interp1d(area_padded, prob_padded, assume_sorted=True)

def __init__(self):
        """ Initialize the Observer object.

        Calls ephem.Observer.__init__ function and adds on gsm
        """
        super(GSMObserver, self).__init__()
        self.gsm = GlobalSkyModel()
        self.observed_sky = None

        # Generate mapping from pix <-> angles
        self.gsm.generate(100)
        self._n_pix  = hp.get_map_size(self.gsm.generated_map_data)
        self._n_side = hp.npix2nside(self._n_pix)
        self._theta, self._phi = hp.pix2ang(self._n_side, np.arange(self._n_pix))

def __init__(self,comp=1.,ramin=None,ramax=None,keepexclude=False,
                 **kwargs):
        self._comp= comp
        if ramin is None: self._ramin= -1.
        else: self._ramin= ramin
        if ramax is None: self._ramax= 361.
        else: self._ramax= ramax
        gaia_tools.select.tgasSelect.__init__(self,**kwargs)
        if not keepexclude:
            self._exclude_mask_skyonly[:]= False
        if not ramin is None:
            theta,phi= healpy.pix2ang(2**5,
                                      numpy.arange(healpy.nside2npix(2**5)),
                                      nest=True)
            ra= 180./numpy.pi*phi
            self._exclude_mask_skyonly[(ra < ramin)+(ra > ramax)]= True
        return None