How to use the healpy.read_map function in healpy

import healpy as hp
import matplotlib.pyplot as plt
import pymaster as nmt

# This script showcases the apodization routine included with pymaster
# and the three apodization modes supported.


# Read input binary mask
mask_raw = hp.read_map("mask.fits", verbose=False)

# The following function calls create apodized versions of the raw mask
# with an apodization scale of 2.5 degrees using three different methods

# Apodization scale in degrees
aposcale = 2.5

# C1 and C2: in these cases, pixels are multiplied by a factor f
#            (with 0<=f<=1) based on their distance to the nearest fully
#            masked pixel. The choices of f in each case are documented in
#            Section 3.4 of the C API documentation. All pixels separated
#            from any masked pixel by more than the apodization scale are
#            left untouched.
mask_C1 = nmt.mask_apodization(mask_raw, aposcale, apotype="C1")
mask_C2 = nmt.mask_apodization(mask_raw, aposcale, apotype="C2")

import healpy as hp
import matplotlib.pyplot as plt
from scipy.special import jv
import os

def opt_callback(option, opt, value, parser):
    setattr(parser.values, option.dest, value.split(','))
parser = OptionParser()
parser.add_option('--plot', dest='plot_stuff', default=False, action='store_true',
                  help='Set if you want to produce plots')
parser.add_option('--nside', dest='nside_out', default=512, type=int,
                  help='Resolution parameter')
(o, args) = parser.parse_args()

l,cltt,clee,clbb,clte,nltt,nlee,nlbb,nlte=np.loadtxt("cls_lss.txt",unpack=True)
msk=hp.ud_grade(hp.read_map("mask_lss.fits",verbose=False),nside_out=o.nside_out); fsky=np.mean(msk)

larr=np.arange(3*o.nside_out)

if not os.path.isfile("cont_lss_nvar_ns%d.fits"%o.nside_out) :
    #First, variations in \bar{n} due to stars
    starmap=hp.read_map("star_template.fits",verbose=False)
    starmap[starmap<=0]=1E-15
    starmap_low=hp.smoothing(hp.ud_grade(hp.ud_grade(starmap,nside_out=64),nside_out=o.nside_out),fwhm=3.5*np.pi/180,verbose=False) #Smooth out the star map to 3.5deg FWHM
    r_starmap=starmap_low/np.amax(starmap_low[msk>0]) #Normalize by maximum value
    shotnoise_factor=(r_starmap-4)/(2*r_starmap-4) #Compute shot noise factor. This interpolates between 1 (no stars) and 1.5 (lots of stars)

    #Now, impose additional depth variations on scales of 3.5 degrees
    cl_extravar=np.exp(-larr*(larr+1)*(3.5*np.pi/180/2.355)**2)
    norm=0.67E-3/np.sum(cl_extravar*(2*larr+1)/4/np.pi); cl_extravar*=norm; #This yields <10% variations
    depth_factor=1+hp.synfast(cl_extravar,o.nside_out,new=True,verbose=False);

def show_healpix_skymap(self):
        """Inserting a valid LIGO/Virgo skymap for sending to Aladin plane."""
        
        try:
            prob = hp.read_map(self.entry_skymap.get(), verbose = False)
            aladin.send_file(self.entry_skymap.get()) 
            aladin.rename(self.entry_skymap.get())
        except ValueError as value_error:
            messagebox.showerror ('Error: 1. Load local file...',
                                    value_error)
        except IOError as io_error:
            messagebox.showerror ('Error: 1. Load local file...',
                                    io_error)
        except SAMPHubError as samphub_error:
            messagebox.showerror ('Error: 1. Load local file...',
                                    samphub_error)

import aladinSAMP
     import aladin_console
     from print_area_prob import print_area_prob
     from table_ipix_percentage import table_ipix_percentage
     from find_highest_probability_pixel import find_highest_probability_pixel
     from instrument_FOV import instrument_FOV
     from FOV_sequence import add_FOV
     from progress_bars import progress_bar
     
     # load a probability healpix sky map
     while True: 
         try:
              print ' Load a probability skymap '
              input_skymap = raw_input( ' : ' ).strip( )
              hpx_test = hp.read_map( input_skymap, verbose = False )
         except IOError as io_error:
               print ''
               print '    -------------------------------------------------------------'
               print '               ***Oops!  No sky map has been found.***'
               print '    -------------------------------------------------------------'
               print ''
               print '', io_error
               print ''
         except ValueError as value_error:  
               print ''
               print '    -------------------------------------------------------------'
               print '               ***Oops!  No sky map has been found.***'
               print '    -------------------------------------------------------------'
               print ''
               print '', value_error
         else:

ax.set_xlabel("$\\ell$",fontsize=17)
    ax.set_ylabel("$C_\\ell$",fontsize=17)
    ax.set_xlim([2,3E4])
    ax.set_ylim([5E-13,5E-6])
    tickfs(ax)
    plt.loglog()
    plt.savefig("../plots_paper/cls_lss.pdf",bbox_inches='tight')
    plt.show()

if (o.whichfig==4) or (o.whichfig==-1) :
    #Plotting LSS contaminant power spectra
    l,cltt,clee,clbb,clte,nltt,nlee,nlbb,nlte=np.loadtxt("cls_lss.txt",unpack=True)
    ll=l[:3*nside_lss]; cl_x=cltt[:3*nside_lss]; cw_x=clee[:3*nside_lss]
    
    msk_s=hp.read_map("mask_lss_ns%d.fits"%nside_lss,verbose=False); fsky=np.mean(msk_s)
    l_s=hp.read_map("cont_lss_star_ns%d.fits"%nside_lss,verbose=False)
    l_d=hp.read_map("cont_lss_dust_ns%d.fits"%nside_lss,verbose=False)
    w_pq,w_pu=hp.read_map("cont_wl_psf_ns%d.fits"%nside_lss,verbose=False,field=[0,1])
    w_sq,w_su=hp.read_map("cont_wl_ss_ns%d.fits"%nside_lss ,verbose=False,field=[0,1])
    cl_d,cw_pe,cw_pb,_,_,_=hp.anafast([l_d*msk_s,w_pq*msk_s,w_pu*msk_s],pol=True,iter=0)/fsky
    cl_s,cw_se,cw_sb,_,_,_=hp.anafast([l_s*msk_s,w_sq*msk_s,w_su*msk_s],pol=True,iter=0)/fsky

    def rebin_and_plot(x,y,b,lt,c,ax,lb=None) :
        xp=np.mean(x.reshape([len(x)//b,b]),axis=1)
        yp=np.mean(y.reshape([len(x)//b,b]),axis=1)
        if lb is None :
            ax.plot(xp,yp,lt,c=c)
        else :
            ax.plot(xp,yp,lt,c=c,label=lb)
    plt.figure(figsize=(7,8))
    plt.subplots_adjust(hspace=0)
    ax1=plt.subplot(311)

def _find_highest_pixel(self, infile):
        """Finding the sky position of the highest probability pixel."""

        hpx = hp.read_map(infile, verbose = False)
        npix = len(hpx)
        nside = hp.npix2nside(npix) 

        ipix_max = np.argmax(hpx)
        hpx[ipix_max] 
        theta, phi = hp.pix2ang(nside, ipix_max)
        
        ra_max = np.rad2deg(phi) 
        dec_max = np.rad2deg(0.5 * np.pi - theta)

        return round(ra_max,5), round(dec_max,5)

def healpix_skymap(self):
        """Inserting a valid LIGO/Virgo probability skymap."""

        connect_with_samp = False
        while not connect_with_samp:
            try:
                print ('--> Load a probability skymap')
                self._skymap = raw_input( ' from a local directory : ' ).strip( )
                print('')
                prob = hp.read_map(self._skymap, verbose = False)
            except IOError as io_error:
                print (io_error)
            except ValueError as value_error:
                print (value_error)
            else:
                connect_with_samp = True

        npix = len(prob) # get healpix resolution
        nside = hp.npix2nside(npix)

        return self.GWsky_config.update({"skymap" : self._skymap,
                                       "nside" : nside})

def probability_inside_FOV( input_skymap, ra_vertices, dec_vertices ):

     """

        Give the probability contained in a specific FOV
     
     """

     
     import healpy as hp
     import numpy as np

     hpx = hp.read_map( input_skymap, verbose = False )

     npix = len( hpx )

     nside = hp.npix2nside( npix )

     theta = 0.5 * np.pi - np.deg2rad( dec_vertices )

     phi = np.deg2rad( ra_vertices )
     
     xyz = hp.ang2vec( theta, phi )

     ipix_poly = hp.query_polygon( nside, xyz )

     # probability contains in a specific FOV
     probability_inside_polygon = hpx[ipix_poly].sum()

Path to the optionally gzip-compressed FITS file.

    nest: bool, optional
        If omitted or False, then detect the pixel ordering in the FITS file
        and rearrange if necessary to RING indexing before returning.

        If True, then detect the pixel ordering and rearrange if necessary to
        NESTED indexing before returning.

        If None, then preserve the ordering from the FITS file.

        Regardless of the value of this option, the ordering used in the FITS
        file is indicated as the value of the 'nest' key in the metadata
        dictionary.
    """
    prob, header = hp.read_map(filename, h=True, verbose=False, nest=nest)
    header = dict(header)

    metadata = {}

    ordering = header['ORDERING']
    if ordering == 'RING':
        metadata['nest'] = False
    elif ordering == 'NESTED':
        metadata['nest'] = True
    else:
        raise ValueError(
            'ORDERING card in header has unknown value: {0}'.format(ordering))

    try:
        value = header['OBJECT']
    except KeyError:

def read_sky_map(filename):
    prob, header = hp.read_map(filename, h=True)
    header = dict(header)

    metadata = {}

    try:
        value = header['OBJECT']
    except KeyError:
        pass
    else:
        metadata['objid'] = value

    try:
        value = header['REFERENC']
    except KeyError:
        pass
    else: