1

我正在应用卷积技术来卷积 2 个数据集，一个 nside = 256 的 healpix 图和一个形状的主光束 (256, 256)，以便测量卷积后的 healpix 图的总强度。我的问题是，在将我的地图与主光束卷积后，我在卷积后的地图中得到了环。我尝试使用 lanczos 或 Gaussian 内核对其进行归一化以处理环，但所有这些方法都失败了。

在下面的代码中，我使用 scipy 中的查询函数在给定半径内的 healpix 地图中搜索最近的像素，并使用地图坐标获取主光束中相应像素的乘积之和。我得到的最终图像中有环。请问谁能帮我解决这个问题？提前致谢。

def query_npix(nside, npix, radius):
    print 'searching for nearest pixels:......'
    t1, t2 = hp.pix2ang(nside, np.arange(npix)) 
    tree = spatial.cKDTree(zip(t1, t2))

    dist, ipix_indx = tree.query(zip(t1, t2), k = 150, distance_upper_bound = radius)
    r1, r2 = hp.pix2ang(nside, ipix_indx)
    ra = r1.T - t1
    dec = r2.T - t2
    print 'Done searching'
    return np.array(dist), np.array(ipix_indx), np.array(ra.T), np.array(dec.T)

def fullSky_convolve(healpix_map, primary_beam_fits, ipix_indx, dist, radius, r1, r2):

    measured_map = []

    hdulist = openFitsFile(primary_beam_fits)
    beam_data = hdulist[0].data
    header = hdulist[0].header  
    nside = hp.get_nside(healpix_map[0, ...])
    npix = hp.get_map_size(healpix_map[0, ...])         # total number of pixels in the map must be  12 * nside^2 

    crpix1, crval1, cdelt1 = [ header.get(x) for x in "CRPIX1", "CRVAL1", "CDELT1" ]
    crpix2, crval2, cdelt2 = [ header.get(x) for x in "CRPIX2", "CRVAL2", "CDELT2" ]

    # beam centres in pixel coordinates
    xc = crpix1-1 + (np.rad2deg(r1.ravel()) - crval1)/(256*cdelt1)
    yc = crpix2-1 + (np.rad2deg(r2.ravel()) - crval2)/(256*cdelt2)
    #xc =  (np.rad2deg(r1.ravel()) )/cdelt1

    for j in xrange(4):
        print 'started Stokes: %d' %j

        for iter in xrange(0 + j, 16, 4):
            outpt = np.zeros(shape = npix, dtype=np.float64)
            #by = outpt.copy()
            # mask beam
            bm_data = beam_data[iter]
            #masked_beam= beam_data[iter]
            shape = bm_data.shape
            rad = np.linspace(-shape[0]/2,shape[-1]/2,shape[0])
            rad2d =  np.sqrt(rad[np.newaxis,:]**2+rad[:,np.newaxis]**2)
            mask = rad2d <= radius/abs(cdelt2)
            masked_beam = bm_data*mask

            s1 = ndimage.map_coordinates(masked_beam, [xc, yc], mode = 'constant')
            bm_map = s1.reshape(dist.shape[0], dist.shape[-1])

        for itr in xrange(npix):
            g_xy = (1.0/(np.sqrt(2*np.pi)*np.std(dist[itr])))*np.exp(-(dist[itr])**2/(2*np.var(dist[itr])))
            #weighted_healpix_map = np.convolve(healpix_map[j, ...][ipix_indx[itr]],  g_xy/g_xy.sum(), mode='same')
            weighted_healpix_map = ndimage.filters.convolve(healpix_map[j, ...][ipix_indx[itr]],  g_xy/g_xy.sum(), mode='reflect')
            #outpt[itr] = np.sum(weighted_healpix_map*(bm_map[itr]/bm_map[itr].sum()))
            outpt[itr] = np.sum(weighted_healpix_map*(bm_map[itr]))
            #print 'itr', itr           
            alpha = file('pap%d.save'%iter, 'wb')
            #h_map = ndimage.filters.gaussian_filter(outpt, sigma = 3.)
            cPickle.dump(outpt, alpha, protocol = cPickle.HIGHEST_PROTOCOL)
            alpha.close()
            print 'Just dumped stripp%d.save:-------'%iter
    print 'Loading dumped files:-------'
    loaded_objects = []
    for itr4 in xrange(16):
        alpha = file('stripp%d.save'%itr4, 'rb')
        loaded_objects.append(cPickle.load(alpha))
        alpha.close()
        measured_map.append(copy.deepcopy(loaded_objects))
    return measured_map
4

1 回答 1

0

请记住，HEALPix 地图可以是“环形”或“嵌套”格式。听起来您可能需要将关键字添加nest=True到诸如hp.pix2ang. 如果您的输入映射是嵌套格式，则需要此关键字。

例如：我最近尝试使用该healpy.smoothing()功能，并在使用healpix.mollview(). mollview使用nested=True关键字运行后，环消失了，图像按我预期的方式呈现。检查您的输入文件使用的排序方案

参考： http ://healpy.readthedocs.org/en/latest/tutorial.html#creating-and-manipulating-maps

Healpix 支持两种不同的排序方案，RING 或 NESTED。默认情况下，healpy 贴图按 RING 顺序排列。为了使用 NESTED 排序，所有与地图相关的函数都支持 nest 关键字，例如： hp.mollview(m, nest=True, title="Mollview image NESTED")

于 2015-11-02T08:55:53.283 回答