astropy.io适合大表的高效元素访问

2 人关注

我正试图用Python和astropy.io从FITS文件的二进制表中提取数据。该表包含一个有超过200万个事件的事件数组。我想做的是将某些事件的时间值存储在一个数组中，这样我就可以对该数组进行分析。我遇到的问题是，在fortran中（使用FITSIO），同样的操作在一个慢得多的处理器上可能只需要几秒钟，而在Python中使用astropy.io进行完全相同的操作却需要几分钟。我想知道瓶颈到底在哪里，是否有更有效的方法来访问各个元素，以确定是否将每个时间值存储在新数组中。以下是我到目前为止的代码。

from astropy.io import fits
minenergy=0.3
maxenergy=0.4
xcen=20000
ycen=20000
radius=50
datafile=fits.open('datafile.fits')
events=datafile['EVENTS'].data
datafile.close()
times=[]
for i in range(len(events)):
    energy=events['PI'][i]
    if energy<maxenergy*1000:
        if energy>minenergy*1000:
            x=events['X'][i]
            y=events['Y'][i]
            radius2=(x-xcen)*(x-xcen)+(y-ycen)*(y-ycen)
            if radius2<=radius*radius:
                times.append(events['TIME'][i])
print times
希望得到任何帮助。我在其他语言中是个不错的程序员，但我以前还没有真正担心过Python的效率问题。我现在选择用Python做这件事的原因是，我在FITSIO和PGPLOT中都使用了fortran，还有一些来自Numerical Recipes的程序，但是我在这台机器上使用的新的fortran编译器不能被说服来产生一个正常工作的程序（有一些32位与64位的问题，等等）。Python似乎有我需要的所有功能（FITS I/O、绘图等），但如果访问一个列表中的单个元素要花很长时间，我将不得不找到另一个解决方案。
非常感谢。


         2
         
         个评论


           
            这很慢的原因是，对于表中的每一行，你都要从表中提取一列（
            
             events['PI']
            
            ），返回一个新的
            
             ndarray
            
            对象，代表从整个表中切出的那一列--这是一个非常慢的操作。  对'X'、'Y'和'TIME'也一样。  然后你用一个Python整数做索引，它必须进入数组，寻找到那个索引，然后把数字值复制出来，转换成一个Python的
            
             float
            
            对象。  另外，如果你使用的是Python 2，
            
             range()
            
            函数正在创建一个200万Python
            
             int
            
            的列表。


           
            user4322779
           
           ：


           
            在Python for Data Analysis第12章中，Wes McKinney建议确保数组在Fortran或C语言中的排序是连续的，以便在对列或行的操作中获得最佳性能。  默认情况下，numpy数组是行连续的，但一些操作如转置可以改变这一点。这可以通过ndarray的flags属性来检查。如果一个数组不在最佳的内存顺序中，可以用copy来改变，如arr.copy('C').flags将Arr改为C（行）内存顺序，或者arr.copy('F').flags将其改为F（列）顺序。


         python


         arrays


         fits


         astropy


        
         
         
          dragoncat16
         
        
        
         发布于
         
         2015-07-09


        1
        
        个回答


          
           
           
            Tom Aldcroft
           
          
          
           发布于
           
           2015-07-09


          已采纳


         0
         
         人赞同


          
           你需要使用 numpy 矢量操作来做这个。  如果没有像numba这样的特殊工具，在Python中做像你所做的那样的大循环总是很慢，因为它是一种解释型语言。  你的程序应该看起来更像。
          
          energy = events['PI'] / 1000.
e_ok = (energy > min_energy) & (energy < max_energy)
rad2 = (events['X'][e_ok] - xcen)**2 + (events['Y'][e_ok] - ycen)**2
r_ok = rad2 < radius**2
times = events['TIMES'][e_ok][r_ok]