python - 创建二进制完整性图

Question

我正在实现项目的目标，但遇到了一个问题:如何创建完整性图？我有大量数据，一个可能有 500.000 个对象的字段，这些对象在我的图中用不同缩放的点表示:

我想创建一个蒙版，我的意思是，将我的图切成小像素，并说如果我在这个像素中有一个对象，我得到的值:1(例如黑色)elif，我没有对象我的像素，我得到的值:0(例如白色)。

我将创建一个掩码，并且可以将每个字段除以该掩码。问题是我不知道如何处理才能做到这一点:/

我创建了第一个脚本，以便对我的数据进行选择。这一个:

#!/usr/bin/python
# coding: utf-8

from astropy.io import fits
from astropy.table import Table
import numpy as np
import matplotlib.pyplot as plt

                ###################################
                # Fichier contenant le champ brut #
                ###################################

filename = '/home/valentin/Desktop/Field52_combined_final_roughcal.fits'

# Ouverture du fichier à l'aide d'astropy  
field = fits.open(filename)   
print "Ouverture du fichier : " + str(filename)       

# Lecture des données fits
tbdata = field[1].data   
print "Lecture des données du fits"            

                ###############################
                # Application du tri sur PROB #
                ###############################

mask = np.bitwise_and(tbdata['PROB'] < 1.1, tbdata['PROB'] > -0.1)  
new_tbdata = tbdata[mask]   
print "Création du Masque"       

            #################################################
            # Détermination des valeurs extremales du champ #
            #################################################

# Détermination de RA_max et RA_min 
RA_max = np.max(new_tbdata['RA'])
RA_min = np.min(new_tbdata['RA'])
print "RA_max vaut :     " + str(RA_max)
print "RA_min vaut :     " + str(RA_min)

# Détermination de DEC_max et DEC_min   
DEC_max = np.max(new_tbdata['DEC'])
DEC_min = np.min(new_tbdata['DEC'])
print "DEC_max vaut :   " + str(DEC_max)
print "DEC_min vaut :   " + str(DEC_min)

            #########################################
            # Calcul de la valeur centrale du champ #
            #########################################

# Détermination de RA_moyen et DEC_moyen
RA_central = (RA_max + RA_min)/2.
DEC_central = (DEC_max + DEC_min)/2.

print "RA_central vaut : " + str(RA_central)
print "DEC_central vaut : " + str(DEC_central)

print " "
print " ------------------------------- "
print " "

        ##############################
        # Détermination de X et de Y #
        ##############################


# Creation du tableau
new_col_data_X = array = (new_tbdata['RA'] - RA_central) * np.cos(DEC_central)
new_col_data_Y = array = new_tbdata['DEC'] - DEC_central
print 'Création du tableau'


# Creation des nouvelles colonnes
col_X = fits.Column(name='X', format='D', array=new_col_data_X)
col_Y = fits.Column(name='Y', format='D', array=new_col_data_Y)
print 'Création des nouvelles colonnes X et Y'


# Creation de la nouvelle table
tbdata_final = fits.BinTableHDU.from_columns(new_tbdata.columns + col_X + col_Y)

# Ecriture du fichier de sortie .fits
tbdata_final.writeto('{}_{}'.format(filename,'mask'))
print 'Ecriture du nouveau fichier mask'

field.close()

好的，工作正常了!但现在，第二部分是这样的:

###################################################
###################################################
###################################################

filename = '/home/valentin/Desktop/Field52_combined_final_roughcal.fits_mask'

print 'Fichier en cours de traitement' + str(filename) + '\n'

# Ouverture du fichier à l'aide d'astropy
field = fits.open(filename)         

# Lecture des données fits  
tbdata = field[1].data    

figure = plt.figure(1)
plt.plot (tbdata['X'], tbdata['Y'], '.')
plt.show()

你知道如何处理吗？我如何在小垃圾箱中切割我的地 block ？

谢谢!更新:

在armatita回答后，我更新了我的脚本:

###################################################
###################################################
###################################################

filename = '/home/valentin/Desktop/Field52_combined_final_roughcal.fits_mask'

print 'Fichier en cours de traitement' + str(filename) + '\n'

# Opening file with astropy
field = fits.open(filename)         

# fits data reading 
tbdata = field[1].data    

##### BUILDING A GRID FOR THE DATA ########
nodesx,nodesy = 360,360   # PIXELS IN X, PIXELS IN Y
firstx,firsty = np.min(tbdata['X']),np.min(tbdata['Y'])
sizex = (np.max(tbdata['X'])-np.min(tbdata['X']))/nodesx
sizey = (np.max(tbdata['Y'])-np.min(tbdata['Y']))/nodesy
grid = np.zeros((nodesx+1,nodesy+1),dtype='bool') # PLUS 1 TO ENSURE ALL DATA IS INSIDE GRID

# CALCULATING GRID COORDINATES OF DATA
indx = np.int_((tbdata['X']-firstx)/sizex)
indy = np.int_((tbdata['Y']-firsty)/sizey)
grid[indx,indy] = True  # WHERE DATA EXISTS SET TRUE

# PLOT MY FINAL IMAGE
plt.imshow(grid.T,origin='lower',cmap='binary',interpolation='nearest')
plt.show()

我发现这个情节:

因此，当我调整 bin 大小时，我可以看到或多或少的空白，这表明我的像素中是否存在对象:)

Answer 1

这通常是将数据插入网格(像素方式或节点方式)的过程。以下示例构建网格(二维数组)并计算样本数据的“网格坐标”。一旦它有了这些网格坐标(实际上只是数组索引)，您就可以将这些元素设置为 True。检查以下示例:

import numpy as np
import matplotlib.pyplot as plt

x = np.random.normal(0,1,1000)
y = np.random.normal(0,1,1000)

##### BUILDING A GRID FOR THE DATA ########
nodesx,nodesy = 100,100   # PIXELS IN X, PIXELS IN Y
firstx,firsty = x.min(),y.min()
sizex = (x.max()-x.min())/nodesx
sizey = (y.max()-y.min())/nodesy
grid = np.zeros((nodesx+1,nodesy+1),dtype='bool') # PLUS 1 TO ENSURE ALL DATA IS INSIDE GRID

# CALCULATING GRID COORDINATES OF DATA
indx = np.int_((x-firstx)/sizex)
indy = np.int_((y-firsty)/sizey)
grid[indx,indy] = True  # WHERE DATA EXISTS SET TRUE

# PLOT MY FINAL IMAGE
plt.imshow(grid.T,origin='lower',cmap='binary',interpolation='nearest')
plt.show()

，结果是:

请注意，我正在使用 imshow 显示图像。如果我减少像素数(20,20 = 节点x，节点y)，我会得到:

另外，对于 matplotlib 中更自动的绘图，您可以考虑 hexbin .