How to use the healpy.pixelfunc function in healpy

H.mollview(sig)

print 'ini: mono=%s, dipole=%s'%(c,d)
mono,dipole = H.pixelfunc.fit_dipole(sig)
print 'fit: mono=%s, dipole=%s'%(mono,dipole)

resfit = dot(vec.T,dipole)+mono

diff = sig.copy()
diff[sig != H.UNSEEN] -= resfit[sig != H.UNSEEN]

H.mollview(diff)

# use remove_dipole
m2 = H.pixelfunc.remove_dipole(sig)
m3 = H.pixelfunc.remove_monopole(sig)

H.mollview(m2)
H.mollview(m3)

w = wcs.WCS(hdr)
x, y = np.meshgrid(np.arange(img.shape[1]),
                   np.arange(img.shape[0]),
                   indexing='xy')

# c = wcs.utils.pixel_to_skycoord(x, y, w, 0)
# Bleeding-edge astropy isn't playing nice. Here's a work-around for now.
l, b = w.all_pix2world(x, y, 0)
c = SkyCoord(l, b, frame='fk5', unit=(u.deg, u.deg))
l = c.galactic.l.radian
b = c.galactic.b.radian

planckfile = 'HFI_SkyMap_857_2048_R1.10_nominal.fits'
planckdata = hp.read_map(planckfile)
#idx = hp.ang2pix(2048,l,b)
planckmap = (hp.pixelfunc.get_interp_val(
    planckdata, np.pi / 2 - b, l, nest=False))
#planckidx = hp.pixelfunc.ang2pix(2048,b,l)


good = np.where(np.isfinite(planckmap) * np.isfinite(img))
m, b = np.polyfit(planckmap[good], img[good], 1)
plt.hexbin(planckmap.ravel(), img.ravel()-b , cmap='gray_r',
           bins='log', gridsize=[100, 100])
print('Offset is {0} MJy/sr'.format(b))
plt.plot(np.linspace(0, 1000, 20), np.linspace(0, 1000, 20))
plt.xlabel(r'Planck 350 $\mu$m')
plt.ylabel(r'Herschel 350 $\mu$m')
plt.xlim([0, np.max(planckmap[good])])
plt.ylim([0, np.max(img[good])])
plt.show()
# plt.savefig('planck_vs_aur.png')

1.-(r+1)*1./nrows+margins[1],
			  1./ncols-margins[2]-margins[0],
			  1./nrows-margins[3]-margins[1])
		extent = (extent[0]+margins[0],
			  extent[1]+margins[1],
			  extent[2]-margins[2]-margins[0],
			  extent[3]-margins[3]-margins[1])

	# Starting to draw : turn interactive off
	wasinteractive = plt.isinteractive()
	plt.ioff()
	try:
		if map is None:
			map = np.zeros(12)+np.inf
			cbar=False
		map = pixelfunc.ma_to_array(map)
		ax=PA.HpxMollweideAxes(f,extent,coord=coord,rot=rot,
						   format=format2,flipconv=flip)
		f.add_axes(ax)
		if remove_dip:
			map=pixelfunc.remove_dipole(map,gal_cut=gal_cut,
									nest=nest,copy=True,
									verbose=True)
		elif remove_mono:
			map=pixelfunc.remove_monopole(map,gal_cut=gal_cut,nest=nest,
									  copy=True,verbose=True)
		img = ax.projmap(map,nest=nest,xsize=xsize,coord=coord,vmin=min,vmax=max,
			   cmap=cmap,norm=norm)
		if cbar:
			im = ax.get_images()[0]
			b = im.norm.inverse(np.linspace(0,1,im.cmap.N+1))
			v = np.linspace(im.norm.vmin,im.norm.vmax,im.cmap.N)

def swap_scheme(self):
        """
        """
        hpx_out = self.hpx.make_swapped_hpx()
        if self.hpx.nest:
            if self.data.ndim == 2:
                data_out = np.vstack([hp.pixelfunc.reorder(
                    self.data[i], n2r=True) for i in range(self.data.shape[0])])
            else:
                data_out = hp.pixelfunc.reorder(self.data, n2r=True)
        else:
            if self.data.ndim == 2:
                data_out = np.vstack([hp.pixelfunc.reorder(
                    self.data[i], r2n=True) for i in range(self.data.shape[0])])
            else:
                data_out = hp.pixelfunc.reorder(self.data, r2n=True)
        return HpxMap(data_out, hpx_out)

Defines the convention of projection : 'astro' (default, east towards left, west towards right)
      or 'geo' (east towards roght, west towards left)
    remove_dip : bool, optional
      If :const:`True`, remove the dipole+monopole
    remove_mono : bool, optional
      If :const:`True`, remove the monopole
    gal_cut : float, scalar, optional
      Symmetric galactic cut for the dipole/monopole fit.
      Removes points in latitude range [-gal_cut, +gal_cut]
    format : str, optional
      The format of the scale label. Default: '%g'
    """
    import pylab

    # Ensure that the nside is valid
    nside = pixelfunc.get_nside(map)
    pixelfunc.check_nside(nside, nest=nest)

    # create the figure (if interactive, it will open the window now)
    f = pylab.figure(fig, figsize=(10.5, 5.4))
    extent = (0.02, 0.25, 0.56, 0.72)
    # Starting to draw : turn interactive off
    wasinteractive = pylab.isinteractive()
    pylab.ioff()
    try:
        if map is None:
            map = np.zeros(12) + np.inf
        map = pixelfunc.ma_to_array(map)
        ax = PA.HpxMollweideAxes(
            f, extent, coord=coord, rot=rot, format=format, flipconv=flip
        )
        f.add_axes(ax)

def _interpolate_cube(self, lon, lat, egy=None, interp_log=True):
        """Perform interpolation on a healpix cube.  If egy is None
        then interpolation will be performed on the existing energy
        planes.

        """

        shape = np.broadcast(lon, lat, egy).shape
        lon = lon * np.ones(shape)
        lat = lat * np.ones(shape)
        theta = np.pi / 2. - np.radians(lat)
        phi = np.radians(lon)
        vals = []
        for i, _ in enumerate(self.hpx.evals):
            v = hp.pixelfunc.get_interp_val(self.counts[i], theta,
                                            phi, nest=self.hpx.nest)
            vals += [np.expand_dims(np.array(v, ndmin=1), -1)]

        vals = np.concatenate(vals, axis=-1)

        if egy is None:
            return vals.T

        egy = egy * np.ones(shape)

        if interp_log:
            xvals = utils.val_to_pix(np.log(self.hpx.evals), np.log(egy))
        else:
            xvals = utils.val_to_pix(self.hpx.evals, egy)

        vals = vals.reshape((-1, vals.shape[-1]))

Returns:
        An array of pixel indices (integers), with the same shape as the input
        SkyCoord coordinates (:obj:`coords.shape`).

    Raises:
        :obj:`dustexceptions.CoordFrameError`: If the specified frame is not supported.
    """
    if coords.frame.name != frame:
        c = coords.transform_to(frame)
    else:
        c = coords

    if hasattr(c, 'ra'):
        phi = c.ra.rad
        theta = 0.5*np.pi - c.dec.rad
        return hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest)
    elif hasattr(c, 'l'):
        phi = c.l.rad
        theta = 0.5*np.pi - c.b.rad
        return hp.pixelfunc.ang2pix(nside, theta, phi, nest=nest)
    elif hasattr(c, 'x'):
        return hp.pixelfunc.vec2pix(nside, c.x.kpc, c.y.kpc, c.z.kpc, nest=nest)
    elif hasattr(c, 'w'):
        return hp.pixelfunc.vec2pix(nside, c.w.kpc, c.u.kpc, c.v.kpc, nest=nest)
    else:
        raise dustexceptions.CoordFrameError(
            'No method to transform from coordinate frame "{}" to HEALPix.'.format(
                frame))

Returns
    -------
    alms : array or tuple of array
      alm or a tuple of 3 alm (almT, almE, almB) if polarized input.

    Notes
    -----
    The pixels which have the special `UNSEEN` value are replaced by zeros
    before spherical harmonic transform. They are converted back to `UNSEEN`
    value, so that the input maps are not modified. Each map have its own,
    independent mask.
    """
    maps = ma_to_array(maps)
    info = maptype(maps)
    nside = pixelfunc.get_nside(maps)
    check_max_nside(nside)

    pixel_weights_filename = None
    if use_pixel_weights:
        if use_weights:
            raise RuntimeError("Either use pixel or ring weights")
        filename = "full_weights/healpix_full_weights_nside_%04d.fits" % nside
        if datapath is not None:
            pixel_weights_filename = os.path.join(datapath, filename)
            if os.path.exists(pixel_weights_filename):
                if verbose:
                    warnings.warn(
                        "Accessing pixel weights from {}".format(pixel_weights_filename)
                    )
            else:
                raise RuntimeError(

wasinteractive = plt.isinteractive()
	plt.ioff()
	try:
		if map is None:
			map = np.zeros(12)+np.inf
			cbar=False
		map = pixelfunc.ma_to_array(map)
		ax=PA.HpxMollweideAxes(f,extent,coord=coord,rot=rot,
						   format=format2,flipconv=flip)
		f.add_axes(ax)
		if remove_dip:
			map=pixelfunc.remove_dipole(map,gal_cut=gal_cut,
									nest=nest,copy=True,
									verbose=True)
		elif remove_mono:
			map=pixelfunc.remove_monopole(map,gal_cut=gal_cut,nest=nest,
									  copy=True,verbose=True)
		img = ax.projmap(map,nest=nest,xsize=xsize,coord=coord,vmin=min,vmax=max,
			   cmap=cmap,norm=norm)
		if cbar:
			im = ax.get_images()[0]
			b = im.norm.inverse(np.linspace(0,1,im.cmap.N+1))
			v = np.linspace(im.norm.vmin,im.norm.vmax,im.cmap.N)
			if matplotlib.__version__ >= '0.91.0':
				cb=f.colorbar(im,ax=ax,
						  orientation='horizontal',
						  shrink=0.5,aspect=25,ticks=PA.BoundaryLocator(),
						  pad=0.05,fraction=0.1,boundaries=b,values=v,
						  format=format)
			else:
				# for older matplotlib versions, no ax kwarg
				cb=f.colorbar(im,orientation='horizontal',