@@ -48,25 +48,25 @@ def freqInfo(fits):
4848
4949
5050def spifit (cube , prefix = None , mask = None , thresh = None ,
51- sigma = 20 , spi_image = None ,
51+ sigma = 20 , spi_image = None ,
5252 spi_err_image = None ):
53-
5453
5554 if isinstance (cube , (list , tuple )):
5655 ims = []
5756 nchan = len (cube )
58- freqs = []
57+ freqs = []
5958 prefix = prefix or cube [0 ][- 3 :]
6059
6160 for im in cube :
6261 with pyfits .open (im ) as hdu :
6362 hdr = hdu [0 ].header
63+ hdu_data = hdu [0 ].data
6464 ndim = hdr ["NAXIS" ]
65- freq_axis = fitsFreqInd (hdr )
66- freqs .append ( hdr ["CRVAL%d" % freq_axis ])
67- ims .append ( im [get_imslice (ndim )] )
65+ ind = fitsFreqInd (hdr )
66+ freqs .append ( hdr ["CRVAL%d" % ind ])
67+ ims .append ( hdu_data [get_imslice (ndim )] )
6868
69- data = numpy .vstack (ims )
69+ data = numpy .dstack (ims ). T
7070 ndim = data .ndim
7171 cnt_freq = freqs [nchan / 2 ]
7272
@@ -78,8 +78,8 @@ def spifit(cube, prefix=None, mask=None, thresh=None,
7878 ndim = data .ndim
7979
8080 freqs , cnt_freq , bw , nchan , ind = freqInfo (cube )
81-
82- freq_ind = ndim - ind
81+
82+ freq_ind = ndim - ind
8383 imslice = get_imslice (ndim )
8484 imslice [freq_ind ] = slice (None )
8585 data = data [imslice ]
@@ -90,7 +90,7 @@ def spifit(cube, prefix=None, mask=None, thresh=None,
9090 aa = []
9191 bb = []
9292 I0 = []
93-
93+
9494 if mask :
9595 with pyfits .open (mask ) as hdu :
9696 mdata = hdu [0 ].data
@@ -104,8 +104,11 @@ def spifit(cube, prefix=None, mask=None, thresh=None,
104104 thresh = sigma * noise
105105 ind = numpy .where (mfs > thresh )
106106
107+ if len (ind ) < 1 :
108+ raise RunTimeError ("No pixels above set threshold, or outside masked region" )
109+
107110 for i ,j in zip (ind [0 ], ind [1 ]):
108- x = numpy .log (freqs / cnt_freq )
111+ x = numpy .log (numpy . array ( freqs ) / cnt_freq )
109112 val = data [:,i ,j ]
110113 if val .any () <= 0 :
111114 continue
@@ -120,7 +123,6 @@ def spifit(cube, prefix=None, mask=None, thresh=None,
120123 # aa.append(intercept)
121124 # bb.append(std_err)
122125 # I0.append(numpy.log(data[nchan/2,i,j]))
123-
124126
125127 nans = numpy .isnan (alpha )
126128 alpha [nans ] = 0.0
0 commit comments