1

I'd like to smooth a map that doesn't cover the full sky. This map is not Gaussian nor has it mean zero, so the default behavior of healpy in which it fills missing values with 0 leads to a bias towards lower values at the edges of this mask:

import healpy as hp

nside = 128
npix = hp.nside2npix(nside)

arr = np.ones(npix)
mask = np.zeros(npix, dtype=bool)

mask[:mask.size//2] = True

arr[~mask] = hp.UNSEEN
arr_sm = hp.smoothing(arr, fwhm=np.radians(5.))

hp.mollview(arr, title='Input array')
hp.mollview(arr_sm, title='Smoothed array')

enter image description here enter image description here

I would like to preserve the sharp edge by setting the weight of the masked values to zero, instead of setting the values to zero. This appears to be difficult because healpy performs the smoothing in harmonic space.

To be more specific, I'd like to mimic the mode keyword in scipy.gaussian_filter(). healpy.smoothing() implicitly uses mode=constant with cval=0, but I'd require something like mode=reflect.

Is there any reasonable way to overcome this issue?

Daniel Lenz
  • 3,334
  • 17
  • 36

2 Answers2

2

The easiest way to handle this is to remove the mean of the map, perform the smoothing with hp.smoothing, then add the offset back. This works around the issue because now the map is zero-mean so zero-filling does not create a border effect.

def masked_smoothing(m, fwhm_deg=5.0): #make sure m is a masked healpy array m = hp.ma(m) offset = m.mean() smoothed=hp.smoothing(m - offset, fwhm=np.radians(fwhm_deg))
return smoothed + offset

The other option I can think of is some iterative algorithm to fill the map in "reflect" mode before smoothing, possibly to be implemented in cythonor numba , the main issue is how complex is your boundary. If it is easy like a latitude cut then all of this is easy, for a general case is very complex and could have a lot of corner cases you need to handle:

Identify "border layers"

  • get all the missing pixels
  • find the neighbors and find which one has a valid neighbor and mark it as the "first border"
  • repeat this algorithm and find pixels that have a "first border" pixel neighbor and mark it as "second border"
  • repeat until you have all the layers you need

Fill reflected values

  • loop on border layers
  • loop on each layer pixel
  • find the valid neighbors, compute their barycenter, now assume that the line between the border pixel center and the barycenter is going perpendicular through the mask boundary and the mask boundary is halfway
  • now extend this line by doubling it in the direction inside the mask, take the interpolated value of the map at that location and assign it to the current missing pixel
  • repeat this for the other layers by playing with the length of the line.
Andrea Zonca
  • 8,378
  • 9
  • 42
  • 70
  • Thanks @Andrea Zonca, these are really good suggestions. I'll play around with that and likely accept the answer later! The mask is indeed complex (think intensity cut + point sources) – Daniel Lenz Apr 28 '18 at 22:34
1

This problem is related to the following question and answer (disclaimer: from me):

https://stackoverflow.com/a/36307291/5350621

It can be transferred to your case as follows:

import numpy as np
import healpy as hp

nside = 128
npix = hp.nside2npix(nside)

# using random numbers here to see the actual smoothing
arr = np.random.rand(npix) 
mask = np.zeros(npix, dtype=bool)
mask[:mask.size//2] = True

def masked_smoothing(U, rad=5.0):     
    V=U.copy()
    V[U!=U]=0
    VV=hp.smoothing(V, fwhm=np.radians(rad))    
    W=0*U.copy()+1
    W[U!=U]=0
    WW=hp.smoothing(W, fwhm=np.radians(rad))    
    return VV/WW

# setting array to np.nan to exclude the pixels in the smoothing
arr[~mask] = np.nan
arr_sm = masked_smoothing(arr)
arr_sm[~mask] = hp.UNSEEN

hp.mollview(arr, title='Input array')
hp.mollview(arr_sm, title='Smoothed array')

David
  • 1,909
  • 1
  • 20
  • 32