How to use the healpy.npix2nside function in healpy

def match_hp_resolution(in_map, nside_out, UNSEEN2nan=True):
    """Utility to convert healpix map resolution if needed and change hp.UNSEEN values to
    np.nan.

    Parameters
    ----------
    in_map : np.array
        A valie healpix map
    nside_out : int
        The desired resolution to convert in_map to
    UNSEEN2nan : bool (True)
        If True, convert any hp.UNSEEN values to np.nan
    """
    current_nside = hp.npix2nside(np.size(in_map))
    if current_nside != nside_out:
        out_map = hp.ud_grade(in_map, nside_out=nside_out)
    else:
        out_map = in_map
    if UNSEEN2nan:
        out_map[np.where(out_map == hp.UNSEEN)] = np.nan
    return out_map

"""

     global  percentage_poly
     
     import healpy as hp
     import numpy as np

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

     # number of pixels
     npix = len(hpx)

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

     # Spherical polar coordinates of vertices in radians
     theta = 0.5 * np.pi - np.deg2rad(dec_vertices)

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

     # Array of indices of pixels inside polygon
     ipix_poly = hp.query_polygon(nside, xyz)

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

def __expand_valid(self, min_p=1e-7):
        #
        # Determine what the 'quanta' of probabilty is
        #
        if self._massp == 1.0:
            # This is to ensure we don't blow away everything because the map
            # is very spread out
            min_p = min(min_p, max(self.skymap))
        else:
            # NOTE: Only valid if CDF descending order is kept
            min_p = self.pseudo_pdf(*self.valid_points_decra[-1])

        self.valid_points_hist = []
        ns = healpy.npix2nside(len(self.skymap))

        # Renormalize first so that the vector histogram is properly normalized
        self._renorm = 0
        # Account for probability lost due to cut off
        for i, v in enumerate(self.skymap >= min_p):
            self._renorm += self.skymap[i] if v else 0

        for pt in self.valid_points_decra:
            th, ph = HealPixSampler.decra2thph(pt[0], pt[1])
            pix = healpy.ang2pix(ns, th, ph)
            if self.skymap[pix] < min_p:
                continue
            self.valid_points_hist.extend([pt]*int(round(self.pseudo_pdf(*pt)/min_p)))
        self.valid_points_hist = numpy.array(self.valid_points_hist).T

def orfFromMap_fast(psr_locs, usermap, response=None):
    """
    Calculate an ORF from a user-defined sky map.

    @param psr_locs:    Location of the pulsars [phi, theta]
    @param usermap:     Provide a healpix map for GW power

    Note: GW directions are in direction of GW propagation
    """
    if response is None:
        pphi = psr_locs[:, 0]
        ptheta = psr_locs[:, 1]

        # Create the pixels
        nside = hp.npix2nside(len(usermap))
        npixels = hp.nside2npix(nside)
        pixels = hp.pix2ang(nside, np.arange(npixels), nest=False)
        gwtheta = pixels[0]
        gwphi = pixels[1]

        # Create the signal response matrix
        F_e = signalResponse_fast(ptheta, pphi, gwtheta, gwphi)
    elif response is not None:
        F_e = response

    # Double the power (one for each polarization)
    sh = np.array([usermap, usermap]).T.flatten()

    # Create the cross-pulsar covariance
    hdcov_F = np.dot(F_e * sh, F_e.T)

import json
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import healpy as hp
import lal
from lalinference.io import fits
from lalinference import plot
from lalinference.bayestar import postprocess

fig = plt.figure(frameon=False)
ax = plt.axes(projection='mollweide' if opts.geo else 'astro hours mollweide')
ax.grid()

skymap, metadata = fits.read_sky_map(opts.input.name, nest=None)
nside = hp.npix2nside(len(skymap))

if opts.geo:
    dlon = -lal.GreenwichMeanSiderealTime(lal.LIGOTimeGPS(metadata['gps_time'])) % (2*np.pi)
else:
    dlon = 0

# Convert sky map from probability to probability per square degree.
probperdeg2 = skymap / hp.nside2pixarea(nside, degrees=True)

# Plot sky map.
vmax = probperdeg2.max()
plot.healpix_heatmap(
    probperdeg2, dlon=dlon, nest=metadata['nest'], vmin=0., vmax=vmax)

# Add colorbar.
if opts.colorbar:

Parameters
    -----------
    goal_dict : dict of healpy maps
        The target goal map(s) being used
    radius : float (1.75)
        Radius of the FoV (degrees)

    Returns
    -------
    Value to use as Target_map_basis_function norm_factor kwarg
    """
    all_maps_sum = 0
    for key in goal_dict:
        good = np.where(goal_dict[key] > 0)
        all_maps_sum += goal_dict[key][good].sum()
    nside = hp.npix2nside(goal_dict[key].size)
    hp_area = hp.nside2pixarea(nside, degrees=True)
    norm_val = radius**2*np.pi/hp_area/all_maps_sum
    return norm_val

def randomPositions(input, nside_pix, n=1):
    """
    Generate n random positions within a full HEALPix mask of booleans, or a set of (lon, lat) coordinates.

    nside_pix is meant to be at coarser resolution than the input mask or catalog object positions
    so that gaps from star holes, bleed trails, cosmic rays, etc. are filled in. 
    Return the longitude and latitude of the random positions and the total area (deg^2).

    Probably there is a faster algorithm, but limited much more by the simulation and fitting time
    than by the time it takes to generate random positions within the mask.
    """
    input = numpy.array(input)
    if len(input.shape) == 1:
        subpix = numpy.nonzero(input)[0] # All the valid pixels in the mask at the NSIDE for the input mask
        lon, lat = pix2ang(healpy.npix2nside(len(input)), subpix)
    elif len(input.shape) == 2:
        lon, lat = input[0], input[1] # All catalog object positions
    else:
        logger.warning('Unexpected input dimensions for skymap.randomPositions')
    pix = surveyPixel(lon, lat, nside_pix)

    # Area with which the random points are thrown
    area = len(pix) * healpy.nside2pixarea(nside_pix, degrees=True)
    
    lon = []
    lat = []
    for ii in range(0, n):
        # Choose an unmasked pixel at random, which is OK because HEALPix is an equal area scheme
        pix_ii = pix[numpy.random.randint(0, len(pix))]
        lon_ii, lat_ii = ugali.utils.projector.pixToAng(nside_pix, pix_ii)
        projector = ugali.utils.projector.Projector(lon_ii, lat_ii)

coordsys = self.config['coords']['coordsys']
        if coordsys.lower() == 'gal':
            print("GAL coordintes")
            objs['GLON'],objs['GLAT'] = lon,lat
            objs['RA'],objs['DEC'] = gal2cel(lon,lat)
        else:
            print("CEL coordintes")
            objs['RA'],objs['DEC'] = lon,lat
            objs['GLON'],objs['GLAT'] = cel2gal(lon,lat)

        modulus = objects['Z_MAX']
        objs['MODULUS'] = modulus
        objs['DISTANCE'] = mod2dist(modulus)

        nside = healpy.npix2nside(len(self.nannulus))
        pix = ang2pix(nside,lon,lat)

        richness = self.richness[objects['IDX_MAX'],objects['ZIDX_MAX']]
        objs['RICHNESS'] = richness
        objs['MASS'] = richness * self.stellar[pix]

        objs['NANNULUS']  = self.nannulus[pix].astype(int)
        objs['NINTERIOR'] = self.ninterior[pix].astype(int)

        # Default name formatting
        # http://cdsarc.u-strasbg.fr/ftp/pub/iau/
        # http://cds.u-strasbg.fr/vizier/Dic/iau-spec.htx
        fmt = "J%(hour)02i%(hmin)04.1f%(deg)+03i%(dmin)02i"
        for obj,_ra,_dec in zip(objs,objs['RA'],objs['DEC']):
            hms = dec2hms(_ra); dms = dec2dms(_dec)
            params = dict(hour=hms[0],hmin=hms[1]+hms[2]/60.,

def contour(m, levels, nest=False, degrees=False, simplify=True):
    try:
        import networkx as nx
    except:
        raise RuntimeError('This function requires the networkx package.')

    # Determine HEALPix resolution
    npix = len(m)
    nside = hp.npix2nside(npix)
    min_area = 0.4 * hp.nside2pixarea(nside)

    # Compute faces, vertices, and neighbors.
    # vertices is an N X 3 array of the distinct vertices of the HEALPix faces.
    # faces is an npix X 4 array mapping HEALPix faces to their vertices.
    # neighbors is an npix X 4 array mapping faces to their nearest neighbors.
    faces = np.ascontiguousarray(
        np.rollaxis(hp.boundaries(nside, np.arange(npix), nest=nest), 2, 1))
    dtype = faces.dtype
    faces = faces.view(np.dtype((np.void, dtype.itemsize * 3)))
    vertices, faces = np.unique(faces.ravel(), return_inverse=True)
    faces = faces.reshape(-1, 4)
    vertices = vertices.view(dtype).reshape(-1, 3)
    neighbors = hp.get_all_neighbours(nside, np.arange(npix), nest=nest)[::2].T

    # Loop over the requested contours.

def count_modes(m, nest=False):
    """Count the number of modes in a binary HEALPix image by repeatedly
    applying the flood-fill algorithm.

    WARNING: The input array is clobbered in the process."""
    npix = len(m)
    nside = hp.npix2nside(npix)
    for nmodes in range(npix):
        nonzeroipix = np.flatnonzero(m)
        if len(nonzeroipix):
            flood_fill(nside, nonzeroipix[0], m, nest=nest)
        else:
            break
    return nmodes