How to use the healpy.nside2pixarea function in healpy

def flat_bitmap(self):
        """Return flattened HEALPix representation."""
        m = np.empty(hp.nside2npix(hp.order2nside(self.order)))
        for nside, full_nside, ipix, ipix0, ipix1, samples in self.visit():
            m[ipix0:ipix1] = len(samples) / hp.nside2pixarea(nside)
        return m

def stellarDensity(infile, nside=256, lon_field='RA', lat_field='DEC'): 
    area = hp.nside2pixarea(nside,degrees=True)
    logger.debug("Reading %s"%infile)
    data = fitsio.read(infile,columns=[lon_field,lat_field])

    lon,lat = data[lon_field],data[lat_field]
    pix = ang2pix(nside,lon,lat)
    counts = collections.Counter(pix)
    pixels, number = np.array(sorted(counts.items())).T
    density = number/area

    return pixels, density

# After overlap, get about 8 sq deg per pointing.

        if extra_features is None:
            self.extra_features = []
            self.extra_features.append(features.Coadded_depth(filtername=filtername,
                                                              nside=nside))
            self.extra_features[0].feature += 1e-5

        super(Smooth_area_survey, self).__init__(basis_functions=basis_functions,
                                                 basis_weights=basis_weights,
                                                 extra_features=self.extra_features,
                                                 smoothing_kernel=smoothing_kernel,
                                                 nside=nside)
        self.filtername = filtername
        pix_area = hp.nside2pixarea(nside, degrees=True)
        block_size = int(np.round(max_area/pix_area))
        self.block_size = block_size
        # Make the dithering solving object
        self.hpc = dithering.hpmap_cross(nside=nside)
        self.max_region_size = np.radians(max_region_size)
        self.nside = nside
        self.percentile_clip = percentile_clip

def getHealpixArea(nside):
    return hp.nside2pixarea(nside, degrees=True)

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)

        inside = False
        while not inside:
            # Apply random offset
            arcminToDegree = 1 / 60.
            resolution = arcminToDegree * healpy.nside2resol(nside_pix, arcmin=True)
            x = 2. * (numpy.random.rand() - 0.5) * resolution # Using factor 2 to be conservative
            y = 2. * (numpy.random.rand() - 0.5) * resolution

def __init__(self, infile,
                 config=None, stellar_mass=None,
                 distance_modulus_extension='DISTANCE_MODULUS', distance_modulus_field='DISTANCE_MODULUS',
                 pix_data_extension='PIX_DATA'):
        """
        Infile is the output of likelihood grid search. 
        """

        self.infile = infile

        print 'Distance Modulus'
        reader = pyfits.open(self.infile)
        self.nside = reader[pix_data_extension].header['NSIDE']
        self.area_pixel = healpy.nside2pixarea(self.nside, degrees=True)
        self.distance_modulus_array = reader[distance_modulus_extension].data.field(distance_modulus_field)
        reader.close()

        print 'Pixels'
        self.pixels = ugali.utils.skymap.readSparseHealpixMap(self.infile, 'LOG_LIKELIHOOD', construct_map=False)[0]
        print 'Log-likelihood'
        self.log_likelihood_sparse = ugali.utils.skymap.readSparseHealpixMap(self.infile, 'LOG_LIKELIHOOD',
                                                                             construct_map=False)[1]
        print 'Richness'
        self.richness_sparse = ugali.utils.skymap.readSparseHealpixMap(self.infile, 'RICHNESS',
                                                                       construct_map=False)[1]
        print 'Richness Limit'
        self.richness_lim_sparse = ugali.utils.skymap.readSparseHealpixMap(self.infile, 'RICHNESS_LIMIT',
                                                                           construct_map=False)[1]
        
        # In case there is only a single distance modulus

self.pixels_target = PixelRegion(self.config['coords']['nside_pixel'],pix)

        # Boolean arrays for selecting given pixels 
        # (Careful, this works because pixels are pre-sorted by query_disc before in1d)
        self.pixel_interior_cut = numpy.in1d(self.pixels, self.pixels_interior)

        # ADW: Updated for more general ROI shapes
        #self.pixel_annulus_cut  = ~self.pixel_interior_cut
        self.pixel_annulus_cut  = numpy.in1d(self.pixels, self.pixels_annulus)

        # # These should be unnecessary now
        # self.centers_lon, self.centers_lat = self.pixels.lon, self.pixels.lat
        # self.centers_lon_interior,self.centers_lat_interior = self.pixels_interior.lon,self.pixels_interior.lat
        # self.centers_lon_target, self.centers_lat_target = self.pixels_target.lon, self.pixels_target.lat

        self.area_pixel = healpy.nside2pixarea(self.config.params['coords']['nside_pixel'],degrees=True) # deg^2
                                     
        """
        self.centers_x = self._centers(self.bins_x)
        self.centers_y = self._centers(self.bins_y)

        self.delta_x = self.config.params['coords']['pixel_size']
        self.delta_y = self.config.params['coords']['pixel_size']
        
        # Should actually try to take projection effects into account for better accuracy
        # MC integration perhaps?
        # Throw points in a cone around full ROI and see what fraction fall in
        self.area_pixel = self.config.params['coords']['pixel_size']**2
        
        self.centers_lon, self.centers_lat = self.projector.imageToSphere(self.centers_x, self.centers_y)
        """

# After overlap, get about 8 sq deg per pointing.

        if extra_features is None:
            self.extra_features = []
            self.extra_features.append(features.Coadded_depth(filtername=filtername,
                                                              nside=nside))
            self.extra_features[0].feature += 1e-5

        super(Smooth_area_survey, self).__init__(basis_functions=basis_functions,
                                                 basis_weights=basis_weights,
                                                 extra_features=self.extra_features,
                                                 smoothing_kernel=smoothing_kernel,
                                                 nside=nside)
        self.filtername = filtername
        pix_area = hp.nside2pixarea(nside, degrees=True)
        block_size = int(np.round(max_area/pix_area))
        self.block_size = block_size
        # Make the dithering solving object
        self.hpc = dithering.hpmap_cross(nside=nside)
        self.max_region_size = np.radians(max_region_size)
        self.nside = nside
        self.percentile_clip = percentile_clip

surface_brightness_population = surface_brightness_population[cut_difficulty_population]
    ellipticity_population = ellipticity_population[cut_difficulty_population]
    position_angle_population = position_angle_population[cut_difficulty_population]
    age_population = age_population[cut_difficulty_population]
    metal_z_population = metal_z_population[cut_difficulty_population]
    mc_source_id_population = mc_source_id_population[cut_difficulty_population]
    """
    
    # Create bonus columns
    
    print("Creating bonus columns...")
    distance_modulus_population = ugali.utils.projector.distanceToDistanceModulus(distance_population)
    hpix_32_population = ugali.utils.healpix.angToPix(32, lon_population, lat_population) # Make sure this matches the dataset

    # Local stellar density
    pixarea = healpy.nside2pixarea(nside_density, degrees=True) * 60.**2 # arcmin^2
    density_population = m_density[ugali.utils.healpix.angToPix(nside_density, lon_population, lat_population)] / pixarea # arcmin^-2

    # Average fracdet within the azimuthally averaged half-light radius
    #m_fracdet_zero = np.where(m_fracdet >= 0., m_fracdet, 0.)
    #m_fracdet_zero = m_fracdet
    r_half = np.degrees(np.arctan2(r_physical_population, distance_population)) # Azimuthally averaged half-light radius in degrees
    fracdet_half_population = meanFracdet(m_fracdet, lon_population, lat_population, r_half)
    fracdet_core_population = meanFracdet(m_fracdet, lon_population, lat_population, 0.1)
    fracdet_wide_population = meanFracdet(m_fracdet, lon_population, lat_population, 0.5)

    # Magnitude limits
    nside_maglim = healpy.npix2nside(len(m_maglim_g))
    pix_population = ugali.utils.healpix.angToPix(nside_maglim, lon_population, lat_population)
    maglim_g_population = m_maglim_g[pix_population]
    maglim_r_population = m_maglim_r[pix_population]