1
2
3
4
5
6
7
8 import numpy as np
9 import pyfits as pf
10 import white_2dmpfit as fitter
11 import matplotlib.pyplot
12 import datetime as datetime
13 import pickle
14 import shelve
15
16 import matplotlib.pyplot as pyplot
17 import scipy.special as sp
18 import scipy.integrate as intg
19 import scipy.interpolate as interpol
20 from operator import itemgetter
21 import collections
22 import math
23
25 """
26 A class to hold all the information from a NIFS-produced data cube. The
27 class is capable of handling both single- and multie-extension FITS files.
28
29 Attributes:
30 scidata - A three-dimensional masked array containing the science data.
31 This has been loaded from either the PRIMARY header (for a
32 single-extension FITS file) or the science header (for a multi-
33 extenstion file.) scidata follows the Python convetion for
34 reference indices, that is, if x and y are the spatial
35 dimensions, and z is the spectral dimension, then data values
36 can be found using scidata[z][y][x].
37 vars - A three-dimensional masked array containing the variances of the
38 data points.
39 dq - A three-dimensional array holding the pixel quality flags for
40 the data. Follows the NIFS convention of 0 = good pixel, 1 = bad
41 (flagged) pixel)
42 xs, dx - The first pixel value, and pixel increment, in the first
43 spatial direction (x). These are simply the values of CRVAL1
44 and CDELT1 from the FITS header, respectively.
45 ys, dy, - As above, but for the second spatial dimension (y) and spectral
46 zs, dz - dimension (z) respectively. They are the values of
47 CRVAL2/CDELT2 and CRVAL3/CDELT3, respectively.
48 xpts - A single-dimensional numpy array containing the list of x-
49 values for this data. These values are computed from the CRVAL1
50 and CDELT1 in the FITS file header, and are stored in the same
51 values as used in the header.
52 ypts, zpts- As above, but for the second spatial dimension (y) and the
53 spectral dimension (z). Computed from CRVAL2/CDELT2 and
54 CRVAL3/CDELT3 respectively.
55 nz,ny,nz - Number of pixels in each dimension.
56 xu, yu - Strings denoting the units of the spatial values. Almost always
57 'arcsec'.
58 zu - String denoting the units of the spectral values. Generally
59 'A' (angstroms) or 'um' (microns).
60 pxoff - Pixel coordinate offset (for use in plotting)
61 pyoff - Pixel coordinate offset (for use in plotting)
62 xptsoff - List of offset pixel coordinates (for use in plotting)
63 yptsoff - List of offset pixel coordinates (for use in plotting)
64
65 du - Units of the data
66
67 Methods:
68 getdata(self) - Returns the data portion of the fitscube as a masked
69 array.
70 save(self) - Saves a copy of the current nifscube object to disk,
71 using Python's pickle module. This is not the fastest
72 or most memory efficient method of doing this (writing
73 a long function to write a text file would be faster),
74 but is the simplest and less prone to read/write errors.
75 """
76
78
79 self.scidata = np.ma.masked_array(data=[], mask=np.ma.nomask)
80 self.vars = np.ma.masked_array(data=[], mask=np.ma.nomask)
81 self.dq = np.array([], dtype=int)
82
83 self.xs = 0.0
84 self.dx = 0.0
85 self.ys = 0.0
86 self.dy = 0.0
87 self.zs = 0.0
88 self.dz = 0.0
89 self.xpts = np.array([])
90 self.ypts = np.array([])
91 self.zpts = np.array([])
92 self.xu = "arcsec"
93 self.yu = "arcsec"
94 self.zu = "A"
95
96 self.pxoff = 0.0
97 self.pyoff = 0.0
98 self.xptsoff = self.xpts - self.pxoff
99 self.yptsoff = self.ypts - self.pyoff
100
101 self.du = 'erg/cm2/s/A'
102
104 """
105 Returns the data content of scidata as a numpy array.
106
107 Inputs:
108 filename - File to save to. Defaults to None, in which case a filename
109 based on a date/time stamp will be used.
110
111 Returns:
112 nscidata - The scidata from the cube, as an (unasked) numpy array.
113 """
114 nscidata = np.ma.getdata(self.scidata)
115 return nscidata
116
117 - def save(self, filename=None):
118 """
119 Save the current fits cube to file. Will be saved using a datetime
120 stamp as the file identifier.
121
122 Inputs:
123 filename - File to save to. Defaults to None, in which case a filename
124 based on a date/time stamp will be used.
125
126 Returns:
127 Nil. A file is saved to the current working directory containing this
128 fitscube instance.
129 """
130 if filename == None:
131 filename = 'nifscube '+str(datetime.datetime.now())+'.nfc'
132
133 outfile = open(filename, 'w')
134 pickle.dump(self, outfile)
135 print 'nifscube.save: File written to '+filename
136 outfile.close()
137
139 """
140 Mask any pixels in the 3D data structure that are assigned negative
141 variance values (this will occur occasionally, if the DQ pixel array
142 hasn't caught all the bad pixels.)
143
144 Inputs:
145 Nil.
146
147 Returns:
148 Nil. The fitscube instance is updated in-situ.
149 """
150
151 self.scidata = np.ma.masked_where(self.vars <= 0.0, self.scidata)
152 self.vars = np.ma.masked_where(self.vars <= 0.0, self.vars)
153
154 - def returnSpaxSpec(self, i, j, wrange='all', vrange=None, restwavl=None):
155 """
156 Returns a spectrum from the specified spaxel, over the designated
157 wavelength or velocity range.
158
159 Inputs:
160 i, j - Spaxel x- and y-indices, respectively.
161 wrange - A two-component array, giving the start and end wavelengths
162 to be plotted. Wavelengths must be passed in the same units
163 as used in the fitscube instance. If set to None,
164 plotSpaxSpec will then look for a designated velocity range
165 and rest wavelength to plot. Can also be passed a special
166 argument, 'all', in which case the entire recorded spectrum
167 of the spaxel will be shown. Defaults to 'all'.
168 vrange - A two-component array, giving the start and end radial
169 velocities to plot. Must be passed in conjunction with a
170 restwavl. If both a wrange and vrange are passed, the wrange
171 will take precendence. Defaults to None.
172 restwavl - The rest wavelength to be used in determination of radial
173 velocities, in m. Will not be considered if a wavelength range
174 is provided instead. Must be passed in conjunction with a
175 velocity range - if no rest wavelength is specified,
176 plotSpaxSpec will return an error. Defaults to None.
177
178 Returns:
179 spec - A one-dimensional masked array, containing the data values
180 for this spaxel spectrum.
181 specvars - The corresponding one-dimensional variance array for this
182 spaxel spectrum.
183 specwavl - A one-dimensional array, holding the wavelength markers
184 (i.e. the x-axis values) for this spaxel spectrum. specwavl
185 will be returned in the same units as the fitscube instance.
186 """
187
188
189 if wrange == None:
190 if vrange == None:
191 raise ValueError('returnSpaxSpec: You have not specified\
192 velocity/wavelength range to retrieve this\
193 spectrum from!')
194 else:
195 if restwavl == None:
196 raise ValueError('returnSpaxSpec: You have passed a \
197 velocity range, but failed to specify \
198 a rest wavelength!')
199 else:
200 if vrange[1] < vrange[0]:
201 vtemp = vrange[1]
202 vrange[1] = vrange[0]
203 vrange[0] = vtemp
204 if self.zu == 'A':
205 wavlfactor = 1e10
206 elif self.zu == 'um':
207 wavlfactor = 1e6
208 startw = velToWavl(vrange[0], restwavl) * wavlfactor
209 stopw = velToWavl(vrange[1], restwavl) * wavlfactor
210 print 'returnSpaxSpec: Converted velocity range:'
211 print str(vrange)
212 print 'to wavelength range:'
213 print str(startw)+' ,'+str(stopw)
214 elif wrange == 'all':
215 if vrange != None:
216 raise ValueError("returnSpaxSpec: You've tried to specify a \
217 velocity range without setting the wavelength \
218 range to None!")
219 startw = self.zpts[0]
220 stopw = self.zpts[-1]
221 else:
222 if wrange[1] < wrange[0]:
223 wtemp = wrange[1]
224 wrange[1] = wrange[0]
225 wrange[0] = wtemp
226 startw = wrange[0]
227 stopw = wrange[1]
228
229 if (startw > self.zpts[-1]) or (stopw < self.zpts[0]):
230 raise ValueError('returnSpaxSpec: The specified \
231 wavelength/velocity range does not \
232 exist in this data cube! Perhaps you \
233 should make sure the units you used \
234 are correct.')
235
236
237 list1 = np.where(self.zpts >= startw)
238 list2 = np.where(self.zpts <= stopw)
239 print list1
240 print list2
241 p = max([np.amin(list1), np.amin(list2)])
242 q = min([np.amax(list1), np.amax(list2)])
243
244 spec = self.scidata[p:q+1, j, i]
245 specvars = self.vars[p:q+1, j, i]
246 specwavl = self.zpts[p:q+1] / wavlfactor
247
248 return spec, specvars, specwavl
249
250 - def plotSpaxSpec(self, i, j, wrange='all', vrange=None, restwavl=None,
251 linewidth=2):
252 """
253 Plots the spectrum from the selected spaxel over the designated
254 wavelength or velocity range.
255
256 Inputs:
257 i, j - Spaxel x- and y-indices, respectively.
258 wrange - A two-component array, giving the start and end wavelengths
259 to be plotted. Wavelengths must be passed in the same units
260 as used in the fitscube instance. If set to None,
261 plotSpaxSpec will then look for a designated velocity range
262 and rest wavelength to plot. Can also be passed a special
263 argument, 'all', in which case the entire recorded spectrum
264 of the spaxel will be shown. Defaults to 'all'.
265 vrange - A two-component array, giving the start and end radial
266 velocities to plot. Must be passed in conjunction with a
267 restwavl. If both a wrange and vrange are passed, the wrange
268 will take precendence. Defaults to None.
269 restwavl - The rest wavelength to be used in determination of radial
270 velocities, in m. Will not be considered if a wavelength range
271 is provided instead. Must be passed in conjunction with a
272 velocity range - if no rest wavelength is specified,
273 plotSpaxSpec will return an error. Defaults to None.
274 linewidth - The line width to be used in the resulting plots. Defaults
275 to 2.
276
277 Returns:
278 Nil. A plot is made to the current figure.
279 """
280
281 spec, specvars, specwavl = self.returnSpaxSpec(i, j, wrange=wrange,
282 vrange=vrange,
283 restwavl = restwavl)
284
285 fig = pyplot.gcf()
286 ax = pyplot.gca()
287 pyplot.errorbar(specwavl, spec,
288 yerr=np.sqrt(specvars),
289 ls='-',marker='o', mec='k', mfc='k', color='k',
290 linewidth=linewidth)
291
293 """
294 Extract an average spectrum from a circular aperture of radius r
295 centered on coordinates (x,y). Spaxels are included inside the aperture
296 based on their centre position.
297
298 No checking is done on the coordinate inputs. This means that you can
299 specify an aperture outside the extent of the cube. If the aperture
300 does not include any spaxel centres, an array of zeros will be returned.
301
302 Inputs:
303 x, y - The x- and y-coordinates of the aperture centre, in
304 units of fitscube.xu and fitscube.yu.
305 r - The aperture radius, in units fitscube.zu.
306 useoffsets - Optional. Boolean value specifying whether to use the
307 offset coordinates (True) or not (False). Defaults
308 to False.
309
310 Returns:
311 spec - The average spectrum in the aperture. This is determined by
312 summing together the spectra from the selected spaxels, and
313 then dividing by the number of selected spaxels. Will be
314 returned in units of fitscube.zu. This will be a
315 one-dimensional array with the same dimensions as fitscube.zpts.
316 """
317
318
319 x = float(x)
320 y = float(y)
321 r = float(r)
322
323
324 if useoffsets:
325 xpts = self.xptsoff
326 ypts = self.yptsoff
327 else:
328 xpts = self.xptsoff
329 ypts = self.yptsoff
330
331
332 i0 = findXIndex(self, x, useoffsets=useoffsets)
333 j0 = findYIndex(self, y, useoffsets=useoffsets)
334
335
336
337
338
339
340
341
342
343
344
345
346
347 workdata = np.ma.copy(self.scidata)
348 X, Y = np.meshgrid(xpts, ypts)
349 for k in range(len(workdata[:,0,0])):
350 workdata[k,:,:] = np.ma.masked_where((X-x)**2 + (Y-y)**2 > r**2,
351 workdata[k,:,:])
352 spec = np.ma.average(np.ma.average(workdata, axis=2), axis=1)
353 return spec
354
357 """
358 Subtract a continuum spectrum from every spaxel in the data cube.
359 Continuum spectrum is formed by averaging over a number of circular
360 apertures, scaling the continuum spectrum to match the science
361 spectrum over the region of interest, then subtracting.
362
363 Inputs:
364 restwavl - The rest wavelength for determining the continuum windows.
365 Must be given in units of m.
366 v1, v2 - The start velocity of the velocity slices. Given in
367 units of km/s.
368 vdelt - The width of the velocity slices. Given in km/s.
369 apers - List of tuples describing each aperture to be used.
370 Should be given as [(x1, y1, r1), (x2, y2, r2), ...].
371 useoffsets -Optional. Boolean value specifying whether to use the
372 offset coordinates (True) or not (False). Defaults to True.
373
374 Returns:
375 Nil. scidata array of fitscube instance is modified in place. Nothing
376 is done to the vars or dq arrays, as it is assumed that the continuum
377 is perfectly known from this procedure.
378 """
379
380
381
382
383
384 w1a = velToWavl(v1, restwavl)
385 w1b = velToWavl(v1+vdelt, restwavl)
386 w2a = velToWavl(v2, restwavl)
387 w2b = velToWavl(v2+vdelt, restwavl)
388 w1a = findZIndex(self, w1a)
389 w1b = findZIndex(self, w1b)
390 w2a = findZIndex(self, w2a)
391 w2b = findZIndex(self, w2b)
392 print 'subtractContSpectrum: Forming continuum over indices ['+str(w1a)+','+str(w1b)+'] and ['+str(w2a)+','+str(w2b)+']'
393
394
395 print 'subtractContSpectrum: Forming continuum....'
396 contspec = np.ma.zeros(self.zpts.shape)
397 for a in range(len(apers)):
398 contspec = contspec + self.extractApertureSpec(apers[a][0],
399 apers[a][1],
400 apers[a][2],
401 useoffsets=useoffsets)
402 contspec = contspec / float(len(apers))
403
404 contlvl = (np.ma.average(contspec[w1a:w1b+1])
405 +np.ma.average(contspec[w2a:w2b+1])) / 2.0
406 print 'subtractContSpectrum: Continuum formed!'
407 print 'Average value in entire continuum spectrum: '+str(np.ma.average(contspec))
408 print 'Average continuum level: '+str(contlvl)
409
410 print 'subtractContSpectrum: Subtracting continuum....'
411
412
413
414
415
416
417
418
419
420
421
422
423
424 scilvl = (np.ma.average(self.scidata[w1a:w1b+1,:,:],axis=0) +
425 np.ma.average(self.scidata[w2a:w2b+1,:,:],axis=0)) / 2.0
426 specfact = scilvl / contlvl
427 contsubtr = np.ma.zeros(self.scidata.shape)
428 for j in range(len(self.ypts)):
429 for i in range(len(self.xpts)):
430 contsubtr[:,j,i] = contspec * specfact[j,i]
431 self.scidata = self.scidata - contsubtr
432
433
435 """
436 A class to hold information by fitting multi-component Gaussians to a
437 NIFS data cube.
438
439 Attributes:
440 restwavl - The rest wavelength used to determine line velocities.
441 vmin - The minimum velocity value used for generating the spectrum to
442 fit to.
443 vmax - The maximum velocity value used for generating the spectrum to
444 fit to.
445 scidata - The scidata from the fitscube, trimmed down so that it only
446 includes the data between vmin and vmax
447 vars - The variance array from the fitscube, trimmed down to the
448 velocity range of interest
449 dq - The DQ pixel quality array from the fitscube, trimmed down to
450 the velocity range of interest
451 uncertest - A 2D masked array, of the same shape as the last two axes of
452 scidata, which contains the uncertainty estimate for each
453 spaxel. This uncertainty estimate is derived by determining the
454 1-sigma standard deviation of the science data in the original
455 data cube, just outside the region of interest.
456 xpts - The x-coordinates from the nifscube
457 ypts - The y-coordinates from the nifscube
458 vpts - The z-coordinates from the nifscube, re-expressed as km/s
459 velocities.
460 fitno - A 2D masked array, of the same shape as the last two axes of
461 scidata, which holds the optimal number of components fit to
462 each spaxel. This has been determined by a statistical F-test.
463 fitorig - A 2D masked array, of the same shape as fitno. This is a copy
464 of fitno, recorded BEFORE any masking if performed (i.e. by
465 maskFits.)
466 dof - A 2D masked array, of the same shape as the spatial axes of
467 scidata, holding lists of length fitmax. Holds the degrees of
468 freedom (dof) for each spaxel when fit with a number of
469 Gaussian components, up to (and including) fitmax components.
470 For example, the dof for a one-Gaussian fit to a given spaxel
471 will be given by dof[j,i][0].
472 chisq - As for dof, however, it holds the chi-squared value (NOT
473 REDUCED) for each fit to the spaxel.
474 fits - A 2D masked array, of the same shape as scidata, holding the
475 fit information for each spaxel. The elements of fits take the
476 form of a list. The first element will be the fit parameters
477 corresponding to a one-Gaussian fit, the next the parameters
478 for a two-Gaussian fit, and so on.
479 perror - As for fits, but holding the formal perror covariance matrix
480 values for each fit.
481 ftestlim - The maximum value allowed for an F-test ratio to be before the
482 P-value is taken. Defaults to 1.4.
483 ftestp - The maximum P-value permitted for two fits to be considered
484 statistically different under the F-test. Defaults to 0.05 (5%).
485
486 Methods:
487 fitVelSpax - Computes a multi-dimensional Gaussian fit to one particular
488 spaxel.
489 plotSpaxSpec - Display the data, and fit components, for a particular
490 spaxel.
491 """
492
494 self.restwavl = 0.0
495 self.vmin = 0.0
496 self.vmax = 0.0
497 self.scidata = np.ma.empty(())
498 self.vars = np.ma.empty(())
499 self.summed = np.ma.empty(())
500 self.sumvars = np.ma.empty(())
501 self.uncertest = np.ma.empty(())
502 self.dq = np.asarray([])
503 self.xpts = np.asarray([])
504 self.dx = 0.0
505 self.xs = 0.0
506 self.ypts = np.asarray([])
507 self.dy = 0.0
508 self.ys = 0.0
509 self.vpts = np.asarray([])
510 self.vs = 0.0
511 self.fitno = np.ma.empty(())
512 self.fitorig = np.ma.empty(())
513 self.dof = np.ma.empty(())
514 self.chisq = np.ma.empty(())
515 self.fits = np.ma.empty(())
516 self.fitsorig = np.ma.empty(())
517 self.perror = np.ma.empty(())
518 self.fitstat = np.ma.empty(())
519 self.fluxes = np.ma.empty(())
520 self.fluxerr = np.ma.empty(())
521 self.ftestlim = 1.4
522 self.ftestp = 0.05
523 self.fitmax = 1
524
525 - def fitVelSpax(self, i, j, comps=1, fixno=False, ignoreamp=False,
526 quiet=False, params=None, vlim=None,
527 wmin=10.0, wpercent=0.10, wperspace=0.25, wmaxcut=50.0,
528 alim=None, fixbg=False, snthresh=5):
529 """
530 Compute a multi-component Gaussian fit to spaxel x=i, y=j in the velfits
531 instance. Modifies the velfits instance in place.
532
533 If params are not provided, the function uses the following algorithm
534 to produce initial guesses:
535 - Make a copy of the spaxel data
536 - Find the highest data value within the range defined by vlim, use to
537 give guess velocity shift and amplitude
538 - Compute estimate for line width
539 - Find nearest point to peak <= half the guess amplitude, compute
540 where the line between the that point, and the next point
541 towards the peak, is = half the guess amplitude
542 - If the next point towards the peak is masked, use the line
543 between the peak and the half-amplitude point instead
544 - Compute the maximum allowed width for this fitted line, based on the
545 guess velocity centroid
546 - If guess width is above this value, adjust guess to 5.0 km/s
547 below max allowed value
548 - If guess width is below wmin, adjust guess to 5.0 km/s above max
549 allowed value
550 - Subtract the guessed Gaussian profile from copy of spaxel data
551 - Repeat the procedure until there are enough guess parameters to fit
552 comps components
553
554 Inputs:
555 velfits - The velfits instance to be acted upon.
556 i - The x-index of the spaxel to be fitted.
557 j - The y-index of the spaxel to be fitted.
558 comps - The number of Gaussian components to fit to the spaxel
559 spectrum. Defaults to 1.
560 fixno - Boolean value denoting whether to force the spectrum to be
561 fit with comps components (True), or to incrementally try a
562 number of components from 1 to comps, stopping when the
563 limiting reduced chi-squared is acheived (False). Defaults
564 to False.
565 ignoreamp - Boolean value specifying whether to consider the nifscube's
566 stated amplitude limit when fitting. Defaults to False (any
567 fits not conforming to the limit will be considered bad).
568 quiet - Boolean value denoting whether or not the MPFIT fitting
569 routine should be quiet or not. Defaults to False.
570 params - Fitting parameters to be passed to the fitting routine.
571 Intended for use when the automatic paramter finder produces
572 unacceptable results. Requires fixno = True.
573 Takes the form of a list composed of
574 [background, amplitude, shift, width, amplitude....].
575 Defaults to None (at which point the
576 automatic first-guess finding will be used.)
577 vlim - A two-element list giving the minimum and maximum values
578 allowed for ALL component centroids in the fits, in units
579 of km/s. Defaults to None (where the minimum and maximum
580 velocities present in the fitvels instance will be used as
581 the limits.
582 wmin - The minimum allowed line sigma to be fit. Note this does NOT
583 represent a FWHM value. Defaults to 10.0
584 wpercent - The percentage of total line intensity allowed to be on the
585 opposite side of v=0 to the line peak. Used to calculate an
586 upper limit on the fitted line width. Must be a float
587 between 0.0 and 1.0. Defaults to 0.10.
588 wperspace - Expected percentage error on the initial guess for the line
589 centroid velocity. Used in the calculation of maximum
590 allowed fitted line width. Must be a float between 0.0 and
591 1.0. Defaults to 0.25.
592 wmaxcut - Lower limit on maximum allowed fitted line width. Avoids the
593 initial guess routine from trying to fit peaks of noise.
594 Defaults to 50.0 km/s.
595 alim - A two-element list giving the minimum and maximum values
596 allowed for the amplitude of ALL fit components, in units
597 km/s. Defaults to None (at which point, the amplitudes
598 will only be limited to be positive).
599 fixbg - Boolean value, specifying if the background height of the
600 fits should be forced to 0.0 (True) or not. Defaults to
601 False.
602 snthresh - Signal-to-noise threshold - if spaxel is below this
603 threshold, fits will still be made, but fitno[j,i] will then
604 be set to 0 afterwards. Masking calculated by the
605 maskFits function.
606
607 Returns:
608 Nil. The velfits instance is updated in place. The fit information for
609 the i-th, j-th pixel will be updated as follows:
610 If fixno = False, the fit information for the first multi-component
611 fit that satisfies the rchisqlim for this velfits instance will be
612 recorded. If all fits, up to and including the comps-component fit,
613 do not satisfy this limit, no fit info (and empty list) will be
614 recorded.
615 If fixno = True, then the fit will be made blindly and recorded,
616 with no regard to the reduced chi-squared limit.
617 """
618
619 if (snthresh < 0.0):
620 raise ValueError('fitVelSpax: negative S/N threshold value passed!')
621
622
623
624
625 if (params != None) and (len(params) != (3*comps + 1)):
626 raise ValueError('fitVelSpax: length of supplied params list does\
627 not correspond to number of components specified!')
628
629
630 if vlim == None:
631 vlim = [self.vpts[0],self.vpts[-1]]
632 elif vlim[1] < vlim[0]:
633 vlim.reverse()
634 if wmin == None:
635 wmin = 10.0
636 if alim == None:
637 alim = [0.0,0.0]
638 alimbool = [True,False]
639 else:
640 alimbool = [True,True]
641 if alim[1] < alim[0]:
642 alim.reverse()
643
644
645
646
647 peakinds = []
648 guessv = []
649 guessa = []
650 guesss = []
651 wlims = []
652 spaxdata = np.ma.copy(self.scidata[:,j,i])
653 for k in range(0, comps):
654
655
656 ind = np.ma.argmax(np.ma.masked_where(self.vpts < vlim[0],
657 np.ma.masked_where(self.vpts>vlim[1], spaxdata)))
658 v = self.vpts[ind]
659
660 a = spaxdata[ind]
661
662
663
664
665 lesshalf = np.ma.where(spaxdata <= (a/2.0))[0]
666 if len(lesshalf) == 0:
667 fwhmguess = wmaxcut
668 else:
669 halfpt = np.ma.argmin(np.ma.abs(lesshalf - ind))
670 halfind = lesshalf[halfpt]
671
672
673
674 if halfind < ind:
675 indshift = 1
676 elif halfind > ind:
677 indshift = -1
678 else:
679 indshift = 1
680
681
682 if spaxdata[halfind+indshift] is np.ma.masked:
683 rise = indshift*(spaxdata[ind] - spaxdata[halfind])
684 run = indshift*(self.vpts[ind] - self.vpts[halfind])
685 else:
686
687
688 rise = indshift*(spaxdata[halfind+indshift] - spaxdata[halfind])
689 run = indshift*(self.vpts[halfind+indshift] - self.vpts[halfind])
690 m = rise / run
691 c = spaxdata[halfind] - m*self.vpts[halfind]
692 intercept = (1.0/m) * ((a/2.0) - c)
693 fwhmguess = 2*np.abs(intercept - v)
694 sigguess = fwhmguess / fwhm
695 wlim = maxWidthAllowed(v, wpercent, wperspace)
696 if sigguess > wlim:
697 sigguess = wlim - 5.0
698 elif sigguess < wmin:
699 sigguess = wmin + 5.0
700 guessv.append(v)
701 guessa.append(a)
702 guesss.append(sigguess)
703 wlims.append(wlim)
704
705
706 spaxdata = spaxdata - fitter.n_gaussian(pars=[0.0,
707 a,v,
708 sigguess])(self.vpts)
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782 if not(quiet):
783 print 'fitVelSpax: Paramater initial guesses:'
784 print 'Amplitude(s): '+str(guessa)
785 print 'Line centroid(s): '+str(guessv)
786 print 'Line sigma(s): '+str(guesss)
787 print 'Line sigma upper limit(s): '+str(wlims)
788
789
790 self.fits[j,i,:,:] = 0.0
791
792 for k in range(2, 3*self.fitmax, 3):
793 self.fits[j,i,:,k] = 1.0
794 self.dof[j,i,:] = 0.0
795 self.chisq[j,i,:] = 0.0
796 self.perror[j,i,:,:] = 0.0
797 self.fitstat[j,i,:] = 0.0
798 self.fluxes[j,i,:] = 0.0
799 self.fluxerr[j,i,:] = 0.0
800
801
802
803 for n in range(1, comps+1):
804
805
806
807 if params == None:
808 guesspars = [0.0]
809 for l in range(0, n):
810 guesspars.append(guessa[l])
811 guesspars.append(guessv[l])
812 guesspars.append(guesss[l])
813 if not(quiet):
814 print 'Generated params: '+str(guesspars)
815 pars = guesspars
816 else:
817 if len(params) > ((3*n)+1):
818 pars = params[0:4*n]
819 else:
820 pars = params
821 if not(quiet):
822 print 'n being used: '+str(n)
823 print 'Generated pars: '+str(pars)
824
825
826 minpars = [0.0]
827 maxpars = [0.0]
828 limitedmin = [False]
829 limitedmax = [False]
830 for k in range(0,n):
831 limitedmin.append(alimbool[0])
832 limitedmin.append(True)
833 limitedmin.append(True)
834 limitedmax.append(alimbool[1])
835 limitedmax.append(True)
836 limitedmax.append(True)
837 minpars.append(alim[0])
838 minpars.append(vlim[0])
839 minpars.append(wmin)
840 maxpars.append(alim[1])
841 maxpars.append(vlim[1])
842 maxpars.append(max([maxWidthAllowed(pars[3*k + 2],wpercent,
843 wperspace),wmaxcut]))
844
845
846
847
848
849
850
851
852
853
854
855
856
857 if not(quiet):
858 print 'minpars, maxpars and pars:'
859 print str(minpars)
860 print str(maxpars)
861 print str(pars)
862 fit = fitter.multigaussfit(self.vpts, self.scidata[:,j,i],
863 err=np.sqrt(self.vars[:,j,i]),
864 quiet=quiet,
865 ngauss=n, params=pars,
866 limitedmin=limitedmin,
867 minpars=minpars,
868 limitedmax=limitedmax,
869 maxpars=maxpars,
870 fixed=[fixbg,False,False,False])
871 if fit.status > 0:
872
873 fits = fit.params.tolist()
874 errs = fit.perror.tolist()
875 bg = fits.pop(0)
876 bgerr = errs.pop(0)
877 fitpars = []
878 fitvels = []
879 errpars = {}
880 errlist = []
881 for q in range(0,n):
882 fitpars.append(fits[3*q:3*(q+1)])
883 fitvels.append(fits[(3*q)+1])
884 errpars[fitvels[q]] = errs[3*q:3*(q+1)]
885 fitpars.sort(key=itemgetter(1))
886 fitvels.sort()
887 for q in range(0,n):
888 errlist.append(errpars[fitvels[q]])
889 fitpars.insert(0, bg)
890 errlist.insert(0, bgerr)
891 fitpars = flatten(fitpars)
892 errlist = flatten(errlist)
893 else:
894 fitpars = fit.params.tolist()
895 errlist = None
896
897 if errlist != None:
898
899
900
901
902
903
904 for k in range(0,3*n+1):
905 self.fits[j,i,n-1,k] = fitpars[k]
906 self.perror[j,i,n-1,k] = errlist[k]
907 if fixbg:
908 self.dof[j,i,n-1] = len(self.vpts) - (3*n)
909 else:
910 self.dof[j,i,n-1] = len(self.vpts) - (3*n) - 1
911 self.chisq[j,i,n-1] = fit.fnorm
912 self.fitstat[j,i,n-1] = fit.status
913
914
915
916
917
918
919 if 0 in self.fitstat[j,i]:
920 self.fitno[j,i] = 0
921 elif fixno:
922 self.fitno[j,i] = comps
923 else:
924 if not(quiet):
925 print 'fitVelSpax: Performing F-test...'
926 n = 1
927 success = False
928 while (n < comps) and (success == False):
929
930 fnum = (self.chisq[j,i][n-1]-self.chisq[j,i][n])
931 fnum = fnum / (self.dof[j,i][n-1]-self.dof[j,i][n])
932 fden = self.chisq[j,i][n] / self.dof[j,i][n]
933 frat = fnum / fden
934 if not(quiet):
935 print 'fitVelSpax: for n '+str(n)+'->'+str(n+1)+', F-test ratio is '+str(frat)
936 if frat < self.ftestlim:
937 success = True
938 self.fitno[j,i] = n
939 else:
940
941
942 fp = sp.fdtrc(self.dof[j,i][n-1]-self.dof[j,i][n],
943 self.dof[j,i][n], frat)
944 if not(quiet):
945 print 'fitVelSpax: for n '+str(n)+'->'+str(n+1)+', F-test prob is '+str(fp)
946 if fp > self.ftestp:
947 success = True
948 self.fitno[j,i] = n
949 else:
950 n += 1
951
952
953 if success == False:
954 self.fitno[j,i] = comps
955
956
957
958
959 for k in range(0, self.fitno[j,i]):
960 lineflux, lineferr = lineFluxKms(self.fits[j,i][self.fitno[j,i]-1][1+3*k],
961 self.fits[j,i][self.fitno[j,i]-1][3+3*k],
962 self.restwavl*1e10,
963 amperr=self.perror[j,i][self.fitno[j,i]-1][1+3*k],
964 widerr=self.perror[j,i][self.fitno[j,i]-1][3+3*k])
965
966
967 self.fluxes[j,i,k] = lineflux
968 self.fluxerr[j,i,k] = lineferr
969
970
971 self.fitorig[j,i] = self.fitno[j,i]
972
973
974 self.maskFits(i, j, snthresh=snthresh)
975
977 """
978 'Mask' spaxel if signal-to-noise is under the specified value. By
979 'mask', the no. of fits recorded for that spaxel will be set to 0.
980 Signal-to-noise is determined as the ratio of the peak flux value
981 recorded in the spaxel to the corresponding error of that data point.
982
983 Inputs:
984 i, j - Spaxel x- and y- indices, respectively.
985 snthresh - The S/N threshold value. Any spaxels with an S/N less than this
986 value will be 'masked'. Defaults to 7.5.
987
988 Returns:
989 Nil. fitvels instance is updated in-situ.
990 """
991
992 self.fitno[j,i] = self.fitorig[j,i]
993
994
995 maxval = self.scidata[:,j,i].max()
996
997
998
999 maxerr = np.sqrt(self.vars[self.scidata[:,j,i].argmax(),j,i])
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010 if (maxval / maxerr) < snthresh:
1011 self.fitno[j,i] = 0
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1025 """
1026 Mask a particular pixel in the velfit instance.
1027
1028 Inputs:
1029 i, j, k - The pixel's x-, y- and spectral indices respectively.
1030
1031 Returns:
1032 Nil. The velfit instance is updated in-situ.
1033 """
1034 maskval = self.scidata[k,j,i]
1035 self.scidata = np.ma.masked_where(self.scidata == maskval, self.scidata)
1036 self.vars = np.ma.masked_where(self.scidata == maskval, self.vars)
1037
1038 - def remask(self, snthresh=5):
1039 """
1040 Apply maskFits to every spaxel in the velfit instances.
1041
1042 Inputs:
1043 snthresh - The S/N threshold for masking. Defaults to 5.
1044
1045 Returns:
1046 Nil. fitvels instance updated in-situ.
1047 """
1048 for i in range(len(self.xpts)):
1049 for j in range(len(self.ypts)):
1050 self.maskFits(i,j,snthresh)
1051
1052
1053
1054 - def plotSpaxSpec(self, i, j, linewidth=2, plotpeak=True, highres=False,
1055 plotresid=True, rsf=0.2, ms=5, arcsecs=False):
1056 """
1057 Plot the spectrum of a given spaxel, along with any fits made.
1058
1059 Inputs:
1060 i - x-index of the spaxel in the fitvel instance
1061 j - y-index of the spaxel in the fitvel instance
1062 linewidth - Line widths of the plot components. Defaults to 2.
1063 plotpeak - Boolean value, denoting whether to plot the peak position
1064 (with errorbars) for each fit component. Defaults to True.
1065 highres - Boolean value, denoting whether to plot the fits at high
1066 velocity resolution (True), or a with the same sampling as
1067 the data (False). Defaults to False.
1068 plotresid - Boolean value, denoting whether to plot the fitting
1069 residual or not. Defaults to True. If plotted, residuals
1070 will be scaled such that they sit between -rsf times the
1071 maximum data value, and 0.0.
1072 rsf - Float value, indicating the scaling that should be applied
1073 to the residuals data. Residuals will then be scaled to
1074 be plotted between -rsf*data.max() and 0.0. Defaults to 0.2.
1075 ms - Float value, denoting the size of the markers used to plot
1076 the data. Defaults to 5.
1077 arcsecs - Boolean value, denoting whether to convert the data values
1078 into flux density per arcsec^2 (True), or leave them as
1079 flux density per pixel (False). Defaults to False.
1080
1081 Outputs:
1082 Nil. The plot is made to the plotting region specified before the
1083 function is called (i.e. a whole figure, a subplot, etc.)
1084 """
1085
1086 if arcsecs:
1087 arcsecfactor = self.dx * self.dy
1088 else:
1089 arcsecfactor = 1.0
1090
1091
1092
1093 fig = pyplot.gcf()
1094 ax = pyplot.gca()
1095 pyplot.errorbar(self.vpts, self.scidata[:,j,i] / arcsecfactor,
1096 yerr=np.sqrt(self.vars[:,j,i]) / arcsecfactor,
1097 ls='None',marker='o', mec='k', mfc='w', color='k',
1098 linewidth=linewidth, ms=ms)
1099
1100 if highres:
1101 xvels = np.asarray(range(int(self.vpts[0]), int(self.vpts[-1]), 1))
1102 else:
1103 xvels = self.vpts
1104
1105
1106 stylelist = ['b-.', 'g:', 'y--']
1107 colorlist = ['b', 'g', 'y']
1108 if self.fitno[j,i] != 0:
1109 pyplot.plot(xvels,
1110 fitter.n_gaussian(pars=self.fits[j,i][self.fitno[j,i]-1])(xvels) / arcsecfactor,
1111 'r--', linewidth=linewidth)
1112 if self.fitno[j,i] > 1:
1113 for k in range(0,self.fitno[j,i]):
1114 parmin = 1 + 3*k
1115 parmax = 4 + 3*k
1116 pyplot.plot(xvels,
1117 fitter.n_gaussian(pars=[self.fits[j,i][self.fitno[j,i]-1][0]] +
1118 self.fits[j,i][self.fitno[j,i]-1][parmin:parmax].tolist())(xvels) / arcsecfactor,
1119 stylelist[k], linewidth=linewidth)
1120 if plotpeak:
1121 peaka = []
1122 peakhalf = []
1123 peakv = []
1124 peakaerr = []
1125 peakverr = []
1126 widths = []
1127 widerr = []
1128 for k in range(0, self.fitno[j,i]):
1129 peaka.append(self.fits[j,i][self.fitno[j,i]-1][1 + 3*k] +
1130 self.fits[j,i][self.fitno[j,i]-1][0])
1131 peakv.append(self.fits[j,i][self.fitno[j,i]-1][2 + 3*k])
1132 widths.append(self.fits[j,i][self.fitno[j,i]-1][3 + 3*k])
1133 peakaerr.append(self.perror[j,i][self.fitno[j,i]-1][1+(3*k)])
1134 peakverr.append(self.perror[j,i][self.fitno[j,i]-1][2+(3*k)])
1135 widerr.append(self.perror[j,i][self.fitno[j,i]-1][3+(3*k)])
1136 pyplot.errorbar(peakv, np.asarray(peaka)/arcsecfactor,
1137 yerr=np.asarray(peakaerr)/arcsecfactor,
1138 xerr=peakverr,
1139 marker='^', mec='red', mfc='w',
1140 color='red', ecolor='red', linewidth=0,
1141 elinewidth=linewidth, ms=ms)
1142 pyplot.errorbar(peakv, ((np.asarray(peaka)+self.fits[j,i][self.fitno[j,i]-1][0])/2.0) / arcsecfactor, yerr=None,
1143 xerr=(fwhm/2.0)*(np.asarray(widths)+np.asarray(widerr)),
1144 fmt=None, mec='red', mfc='w', color='red',
1145 ecolor='red', linewidth=0, elinewidth=linewidth,
1146 capsize=5, ms=0)
1147 pyplot.errorbar(peakv, ((np.asarray(peaka)+self.fits[j,i][self.fitno[j,i]-1][0])/2.0) / arcsecfactor,
1148 yerr=None,
1149 xerr=(fwhm/2.0)*(np.asarray(widths)-np.asarray(widerr)),
1150 fmt=None, mec='red', mfc='w', color='red',
1151 ecolor='red', linewidth=0, elinewidth=linewidth,
1152 capsize=5)
1153
1154
1155 if (plotresid and self.fitno[j,i] > 0):
1156
1157 maxval = self.scidata[:,j,i].max()
1158
1159 resids = self.scidata[:,j,i] - fitter.n_gaussian(pars=self.fits[j,i][self.fitno[j,i]-1])(self.vpts)
1160
1161 sf = rsf * (maxval / resids.max())
1162 resids = sf * resids
1163
1164 resids = resids - resids.max()
1165
1166 pyplot.plot(self.vpts, resids/arcsecfactor, color='purple', ls='-', linewidth=1)
1167
1168
1169
1170
1171 ax.set_title(r'Spectrum for Spaxel ('+str(i)+r','+str(j)+r')')
1172
1173
1174
1175 - def plotImg(self, z=100, colorbar=True, colormap=matplotlib.cm.jet,
1176 log=False, edged=False):
1177 """
1178 Plots a pseudocolor image of the specified parameter in the velfit
1179 instance.
1180
1181 Inputs:
1182 z - An integer code, corresponding to the parameter from the
1183 velfit that we wish to plot. Defaults to 0 (corresponds to
1184 plotting the total flux per spaxel, i.e. self.summed. The
1185 codes are:
1186 0 - Number of fit components
1187 1 - Spaxel chi-squared
1188 2 - Spaxel degrees of freedom
1189 3 - Spaxel reduced chi-squared
1190 100 - Spaxel total flux
1191 101 - Fit background height
1192 102 - Spaxel peak flux
1193 103 - Spaxel summed error
1194 104 - Ratio of spaxel peak flux/spaxel average error
1195 105 - Ratio of spaxel summed flux/spaxel summed error
1196 106 - Ratio of spaxel peak flux/corresponding pixel error
1197 111 - First component total flux
1198 112 - Second component total flux
1199 ...and so on
1200 121 - Spaxel fitted total flux
1201 122 - Spaxel fitted total flux error
1202 123 - Spaxel fitted total flux / total flux error
1203 11 - First component peak flux
1204 12 - First component line velocity
1205 13 - First component line width
1206 21 - Second component peak flux
1207 22 - Second component line velocity
1208 23 - Second component line width
1209 - ....and so on. Codes are currently included for up to
1210 three components.
1211 colorbar - A boolean value, denoting whether or not to include a
1212 colorbar. Defaults to True.
1213 colormap - A colormap instance, to be used as the colormap for the
1214 pseudocolor diagram. Must be passed as matplotlib.cm.*.
1215 Defaults to matplotlib.cm.jet.
1216 log - A boolean value, specifying whether to plot the data as is
1217 (False), or the base-10 log of the data (True). If True,
1218 the data parameter to be plotted will be masked wherever
1219 it is less than or equal to zero. Defaults to False.
1220 edged - Boolean vale, denoting whether to draw edges around each
1221 spaxel in the returned image. Defaults to False.
1222
1223 Returns:
1224 plotted - The colorbar instance which has been plotted. This is
1225 required in order to add to the diagram later on, e.g.
1226 colorbar(plotted).
1227 Will plot a pseudocolor diagram of the requested parameter to the
1228 current figure/axes when called. If the requested parameter does not
1229 exist (i.e. asked to plot the third component amplitude on a two-
1230 component velfit), a message will be returned to the terminal, and
1231 nothing else will occur.
1232 """
1233
1234
1235
1236
1237 if (z < 100) and (z/10 > self.fitmax):
1238 print 'plotImg: The requested fit parameter '+str(z)+' does not\
1239 exist, as the maximum number of components fit = '+str(self.fitmax)\
1240 +'. Aborting plotting.'
1241 return
1242
1243 if edged == True:
1244 edgec = 'k'
1245 else:
1246 edgec = None
1247
1248 fig = pyplot.gcf()
1249 ax = pyplot.gca()
1250
1251
1252 X, Y = gridPtsIndep(self.xpts, self.ypts, self.dx, self.dy, mesh=True)
1253
1254
1255 toplot = np.ma.empty_like(self.fitno)
1256 if z == 100:
1257 toplot = self.summed
1258 plttitle = r'Summed flux'
1259
1260 elif z == 101:
1261 toplot = self.genFitInfoArray(0)
1262 plttitle = r'Fit background height'
1263
1264 elif z == 102:
1265 toplot = self.scidata.max(axis = 0)
1266 plttitle = r'Spaxel peak flux'
1267
1268 elif z == 103:
1269 toplot = np.sqrt(self.vars).sum(axis=0)
1270 plttitle = r'Spaxel summed error'
1271
1272 elif z == 104:
1273 toplot = self.scidata.max(axis=0)/ \
1274 (np.sqrt(self.vars).sum(axis=0)/np.sqrt(self.vars).count(axis=0))
1275 plttitle = r'Spaxel peak flux / spaxel average error'
1276 elif z == 105:
1277 toplot = self.summed / np.sqrt(self.vars).sum(axis=0)
1278 plttitle = r'Spaxel total flux / spaxel total error'
1279 elif z == 106:
1280 numer = self.scidata.max(axis = 0)
1281 denom = np.ma.empty_like(numer)
1282 for i in range(0, len(self.xpts)):
1283 for j in range(0, len(self.ypts)):
1284 denom[j,i] = self.vars[np.where(self.scidata[:,j,i] ==
1285 numer[j,i]),j,i]
1286 toplot = numer / np.sqrt(denom)
1287 plttitle = r'Spaxel peak flux / corresponding pixel error'
1288 elif (z > 110) and (z < 120):
1289 toplot = self.genFitInfoArray(-(z-110))
1290 component = str(z%110)
1291 plttitle = r'Component '+component+r' fitted line intensity'
1292
1293 elif (z > -1) and (z < 4):
1294 toplot = self.genInfoArray(z)
1295 if z == 0:
1296 plttitle = r'Number of fit components'
1297 elif z == 1:
1298 plttitle = r'Spaxel $\chi^2$'
1299 elif z == 2:
1300 plttitle = r'Spaxel DOF'
1301 elif z == 3:
1302 plttitle = r'Spaxel reduced $\chi^2$'
1303 elif (z > 10) and (z < 100) and (z%10 > 0) and (z%10 < 4):
1304 infoindex = 3*(z/10 - 1) + z%10
1305 toplot = self.genFitInfoArray(infoindex)
1306 component = str(z/10)
1307 if z%10 == 1:
1308 plttitle = r'Component '+component+r' fitted amplitude'
1309
1310 elif z%10 == 2:
1311 plttitle = r'Component '+component+r' fitted velocity cent.'
1312
1313 elif z%10 == 3:
1314 plttitle = r'Component '+component+r' fitted line width'
1315
1316 elif z == 121:
1317 toplot = self.sumFitFlux()[0]
1318 plttitle = r'Total fitted line intensity'
1319
1320 elif z == 122:
1321 toplot = self.sumFitFlux()[1]
1322 plttitle = r'Total fitted line intensity error'
1323
1324 elif z == 123:
1325 toplot = self.sumFitFlux()[0] / self.sumFitFlux()[1]
1326 plttitle = r'Total fitted line intensity S/N'
1327 else:
1328 raise ValueError('Unknown parameter code '+str(z)+'passed.')
1329
1330
1331
1332
1333 plotted = ax.pcolor(X, Y, toplot, cmap=colormap, edgecolor=edgec, antialiased=True)
1334
1335
1336
1337
1338
1339
1340 ax.set_xlim(X[0,0], X[-1,-1])
1341 ax.set_ylim(Y[0,0], Y[-1,-1])
1342 ax.set_aspect(1)
1343
1344
1345 ax.set_xlabel(r'x (")')
1346 ax.set_ylabel(r'y (")')
1347 ax.set_title(plttitle)
1348
1349
1350 return plotted
1351
1353 """
1354 Sums over the velocity range of the cube, and save the result to the
1355 summed array.
1356
1357 Inputs: Nil.
1358
1359 Returns: Nil. The summed attribute of the velfit instance is updated.
1360 """
1361 nx = len(self.xpts)
1362 ny = len(self.ypts)
1363
1364 self.summed = np.ma.zeros((ny, nx))
1365 self.sumvars = np.ma.zeros((ny,nx))
1366
1367 for i in range(0,nx):
1368 for j in range(0,ny):
1369 self.summed[j,i] = self.scidata[:,j,i].sum()
1370 self.sumvars[j,i] = np.sqrt(self.vars[:,j,i]).sum()
1371
1372
1373 if (np.ma.getmask(self.scidata)[:,j,i]).all():
1374 self.summed[j,i] = np.ma.masked
1375 self.sumvars[j,i] = np.ma.masked
1376 print 'sumCube: completed row '+str(i+1)+' of '+str(nx)
1377
1378 self.sumvars = self.sumvars**2
1379
1381 """
1382 Collapse the fluxes and fluxerr arrays to represent total flux, summed
1383 over all components that have been fit.
1384
1385 Inputs:
1386 Nil.
1387
1388 Returns:
1389 fluxsum - The array holding the summed flux values.
1390 fluxsumerr - The array holding the summex flux values.
1391 """
1392
1393 fluxsum = np.ma.empty_like(self.fitno, dtype='float')
1394 fluxsumerr = np.ma.empty_like(self.fitno, dtype='float')
1395
1396 for i in range(0, len(self.xpts)):
1397 for j in range(0, len(self.ypts)):
1398 fluxsum[j,i] = sum(self.fluxes[j,i])
1399 fluxsumerr[j,i] = sum(self.fluxerr[j,i])
1400
1401
1402 fluxsum = np.ma.masked_where(self.fitno == 0, fluxsum)
1403 fluxsumerr = np.ma.masked_where(self.fitno == 0, fluxsumerr)
1404
1405 return fluxsum, fluxsumerr
1406
1407
1409 """
1410 Returns the F-test information for all fits made to spaxel (i,j).
1411
1412 Inputs:
1413 i, j - The x- and y-indices of the spaxel to be inspected.
1414
1415 Returns:
1416 ftest - A list of list containing the F test information. Each element
1417 of the master list will be a three-element list, containing:
1418 - The number of components in the more complex fit
1419 - The F-ratio of this complex fit to the simpler fit before
1420 - The F-test P-ratio for the same comparison.
1421 """
1422 if self.fitmax < 2:
1423 print 'returnFTestInfo: Only one component has been fit to this cube!'
1424 return
1425
1426 ftest = []
1427 for k in range(1, self.fitmax):
1428 felem = []
1429 felem.append(k+1)
1430 fnum = (self.chisq[j,i][k-1]-self.chisq[j,i][k])
1431 fnum = fnum / (self.dof[j,i][k-1]-self.dof[j,i][k])
1432 fden = self.chisq[j,i][k] / self.dof[j,i][k]
1433 frat = fnum / fden
1434 felem.append(frat)
1435 fp = sp.fdtrc(self.dof[j,i][k-1]-self.dof[j,i][k],
1436 self.dof[j,i][k], frat)
1437 felem.append(fp)
1438 ftest.append(felem)
1439
1440 return ftest
1441
1442
1444 """
1445 Returns fit information (e.g. chi-sq, degrees of freedom) for
1446 the comps-component fit of the given spaxel.
1447
1448 Inputs:
1449 i - The spaxel x-index
1450 j - The spaxel y-index
1451 comps - The parameter will be extracted from the fit information with
1452 this number of components.
1453 info - An integer code corresponding to the information wanted. Options
1454 are:
1455 0 - number of fit components
1456 1 - spaxel chi-squared
1457 2 - spaxel degrees of freedom
1458 3 - spaxel reduced chi-squared
1459 4 - spaxel F-test information
1460
1461 Returns:
1462 value - The requested information for this spaxel
1463 """
1464 comps = int(comps)
1465 info = int(info)
1466 if (comps < 0) or (comps > self.fitmax):
1467 raise ValueError('returnInfo: This velfit instance does not have fit info\
1468 for '+str(comps)+' components (max. no of\
1469 components = '+str(self.fitmax)+')')
1470 if (info < 0) or (info > 3):
1471 raise ValueError('returnInfo: Invalid info code: '+str(info))
1472
1473
1474
1475 if comps == 0:
1476 return -999
1477
1478 if info == 0:
1479 value = self.fitno[j,i]
1480 elif info == 1:
1481 value = self.chisq[j,i][self.fitno[j,i]-1]
1482 elif info == 2:
1483 value = self.dof[j,i][self.fitno[j,i]-1]
1484 elif info == 3:
1485 value = (self.chisq[j,i][self.fitno[j,i]-1] /
1486 self.dof[j,i][self.fitno[j,i]-1])
1487
1488 return value
1489
1490
1491
1493 """
1494 Returns a particular fit parameter for a given spaxel.
1495
1496 Inputs:
1497 i - The spaxel x-index
1498 j - The spaxel y-index
1499 comps - The parameter will be extracted from the fit information with
1500 this number of components.
1501 info - An integer code corresponding to the information wanted. Options
1502 are:
1503 0 - background amplitude
1504 1 - fit amplitude, component A
1505 2 - fit displacement (i.e. line velocity), component A
1506 3 - fit FWHM (i.e. line width), component A
1507 4 - fit amplitude, component B
1508 5 - fit displacement (i.e. line velocity), component B
1509 ...and so on.
1510
1511 Returns:
1512 parameter - The requested parameter from the relevant fit for this
1513 spaxel.
1514 perror - The formal covariance matrix error associated with this
1515 parameter for this spaxel. Note that this is NOT the
1516 formal error - the relevant chi-squared value and
1517 degrees of freedom will still need to be taken into account.
1518 """
1519
1520
1521 comps = int(comps)
1522 info = int(info)
1523
1524
1525
1526 if comps == 0:
1527 return -999, -999
1528
1529 if (comps < 0) or (comps > self.fitmax):
1530 raise ValueError('returnFitInfo: This velfit instance does not have fit info\
1531 for '+str(comps)+' components (max. no of\
1532 components = '+str(self.fitmax)+')')
1533 if (info < -102) or (info > 3*comps) or (info < -comps and info > -100):
1534 raise ValueError('returnFitInfo: Invalid info code: '+str(info))
1535
1536 if info >= 0:
1537 parameter = self.fits[j,i][comps-1][info]
1538 perror = self.perror[j,i][comps-1][info]
1539 elif info > -100:
1540 parameter = self.fluxes[j,i][-info-1]
1541 perror = self.fluxerr[j,i][-info-1]
1542 else:
1543 parameter = sum(self.fluxes[j,i])
1544 perror = sum(self.fluxerr[j,i])
1545 print parameter, perror
1546 return parameter, perror
1547
1549 """
1550 Generates a masked array of fit information, as returned by returnInfo.
1551 The information returned will correspond to the best fit for each
1552 spaxel (as denoted in self.fitno).
1553
1554 Inputs:
1555 info - An integer code denoting what information we want returned.
1556 These codes correspond to those found in returnInfo().
1557
1558 Returns:
1559 infoarr - A two-dimensional masked array, containing the requested
1560 information for each spaxel. Will take the same shape as the
1561 spatial dimension of the fitvels instance (i.e. the same
1562 shape as, say, fitno).
1563 """
1564
1565 info = int(info)
1566 print 'genInfoArray: Info code '+str(info)
1567 if (info < 0) or (info > 3):
1568 raise ValueError('genInfoArray: Invalid info code: '+str(info))
1569
1570 infoarr = np.ma.empty_like(self.fitno, dtype=float)
1571 for i in range(len(self.xpts)):
1572 for j in range(len(self.ypts)):
1573 infoarr[j,i] = self.returnInfo(i, j, self.fitno[j,i], info)
1574
1575
1576 infoarr = np.ma.masked_where(infoarr == -999, infoarr)
1577
1578 print infoarr
1579 return infoarr
1580
1582 """
1583 Generates a masked array of fit parameters, as returned by
1584 returnFitInfo. The fit parameter returned will correspond to the best
1585 fit for each spaxel (as denoted in self.fitno).
1586
1587 Inputs:
1588 info - An integer code denoting what information we want returned.
1589 These codes correspond to those found in returnFitInfo().
1590
1591 Returns:
1592 paramarr - A two-dimensional masked array, containing the requested
1593 information for each spaxel. Will take the same shape as the
1594 spatial dimension of the fitvels instance (i.e. the same
1595 shape as, say, fitno).
1596 paramperr - A two-dimensional masked array, containing the formal
1597 1-sigma errors for the requested information.
1598 """
1599 info = int(info)
1600 print 'genFitInfoArray: Info code '+str(info)
1601
1602
1603
1604 paramarr = np.ma.empty_like(self.fitno, dtype=float)
1605
1606 print str(self.returnFitInfo(0, 0, self.fitno[0,0], info))
1607
1608 for i in range(len(self.xpts)):
1609 for j in range(len(self.ypts)):
1610
1611
1612 if (info > 3*self.fitno[j,i]) or (info < -self.fitno[j,i]):
1613 paramarr[j,i] = -1.0e10
1614 else:
1615 paramarr[j,i] = (self.returnFitInfo(i, j, self.fitno[j,i], info))[0]
1616 paramarr = np.ma.masked_where(paramarr <= -1.0e9, paramarr)
1617 paramarr = np.ma.masked_where(paramarr == -999, paramarr)
1618
1619 return paramarr
1620
1621 - def save(self, filename=None):
1622 """
1623 Save the current velfit to file.
1624
1625 Inputs:
1626 filename - Filename to be saved to. Defaults to None, in which case
1627 a datetime stamp will be used to generate a file name.
1628
1629 Returns:
1630 Nil. A file is saved to the current working directory containing this
1631 velfit instance.
1632 """
1633 if filename == None:
1634 filename = 'velfit_'+str(datetime.datetime.now())+'.vft'
1635 outfile = open(filename, 'wb')
1636 pickle.dump(self, outfile)
1637 print 'velfit.save: File written to '+filename
1638 outfile.close()
1639
1641 """
1642 Find the spaxel index corresponding to the given velocity value.
1643
1644 Inputs:
1645 z - The x-coordinate that we wish to find the spaxel index for
1646
1647 Returns:
1648 k - The spaxel index corresponding to the passed velocity.
1649 """
1650 vpts = self.vpts
1651 dv = vpts[1]-vpts[0]
1652
1653
1654 if (v < (vpts[0] - dv/2.0)) or (v > (vpts[-1] + dv/2.0)):
1655 raise ValueError('The velocity '+str(v)+' is outside the range of this velfit instance!')
1656
1657
1658
1659 vbases = vpts - (dv/2.0)
1660
1661
1662 coords = np.where(vbases < v)[0]
1663 k = coords[-1]
1664 return k
1665
1666
1667
1668
1669
1670 -def loadFITS(filename, pxoff=0.0, pyoff=0.0, scihead='SCI', varhead='VAR',
1671 dqhead='DQ', xu='arcsec', yu='arcsec', zu='A', scalefact = 1.0,
1672 flipx=False,
1673 quiet=False):
1674 """
1675 Loads a NIFS FITS file into computer memory.
1676
1677 Inputs:
1678 filename - The FITS file to be loaded into memory. May be a single- or
1679 multi-extenstion FITS file.
1680 pxoff - X-pixel offset to be used for plotting purposes. Defaults to 0.0.
1681 pyoff - Y-pixel offset to be used for plotting purposes. Defaults to 0.0.
1682 scihead - The header name of the science data extentsion. Defaults to
1683 'SCI' for multi-extenstion FITS files, and will be automatically
1684 be set to 'PRIMARY' for a single-extension file.
1685 varhead - The header name of the variance data extenstion. Defaults to
1686 'VAR'.
1687 dqhead - The header name of the pixel quality flag array. Defaults to
1688 'DQ'.
1689 xu, yu - Strings defining the units of the x and y axes. Defaults to
1690 'arcsec' for both.
1691 zu - String defining the units of the x and y axes. Defaults to 'A'
1692 (angstrom). This module currently only recognises two options:
1693 'A' - angstrom
1694 'um'- micron
1695 scalefact - Float giving the value to scale the data and variances by. This
1696 is useful for pulling data values into something close to
1697 reasonable human numbers. Defaults to 1.0 (no scaling).
1698 flipx - Boolean value, specifying whether to invert the x-axis (i.e.
1699 flip the sign on the x pixel increment.) Designed for use with
1700 files coming out of IRAF that have had their CDELT1 value sign
1701 flipped for some reason. Defaults to False.
1702
1703
1704 Outputs:
1705 data - A fitscube instance containing the data cube.
1706 """
1707 if not(quiet):
1708 print 'loadFITS: Beginning loading...'
1709
1710 hdulist = pf.open(filename)
1711
1712 if len(hdulist) == 1:
1713 single = True
1714 scihead = 'PRIMARY'
1715 else:
1716 single = False
1717
1718 scalefact = float(scalefact)
1719
1720
1721 data = nifscube()
1722
1723 if not(quiet):
1724 print 'loadFITS: Populating data arrays...'
1725
1726 data.scidata = np.ma.masked_array(data = hdulist[scihead].data,
1727 mask=np.ma.nomask)
1728 if not(single):
1729 data.dq = np.array(hdulist[dqhead].data, dtype=bool)
1730 mask = data.dq
1731 data.vars = np.ma.masked_array(data = hdulist[varhead].data,
1732 mask = mask)
1733 data.scidata = np.ma.masked_array(np.ma.getdata(data.scidata),
1734 mask = mask)
1735 else:
1736 data.vars = np.ma.zeros(data.scidata.shape)
1737 data.dq = np.ma.zeros(data.scidata.shape, dtype=bool)
1738
1739 data.scidata = data.scidata * scalefact
1740 data.vars = data.vars * scalefact * scalefact
1741 if not(quiet):
1742 print 'loadFITS: Computing data coordinates...'
1743
1744 data.nz = data.scidata.shape[0]
1745 data.ny = data.scidata.shape[1]
1746 data.nx = data.scidata.shape[2]
1747 if ((data.nx != hdulist[scihead].header['NAXIS1']) or
1748 (data.ny != hdulist[scihead].header['NAXIS2']) or
1749 (data.nz != hdulist[scihead].header['NAXIS3'])):
1750 raise ValueError("loadFITS: The size of your science data array does\
1751 not match the axis length(s) in the header! I suggest\
1752 you manually inspect your data cube to determine the\
1753 issue.")
1754 try:
1755 data.dx = hdulist[scihead].header['CDELT1']
1756 except KeyError:
1757 data.dx = hdulist[scihead].header['CD1_1']
1758 try:
1759 data.dy = hdulist[scihead].header['CDELT2']
1760 except KeyError:
1761 data.dy = hdulist[scihead].header['CD2_2']
1762 try:
1763 data.dz = hdulist[scihead].header['CDELT3']
1764 except KeyError:
1765 data.dz = hdulist[scihead].header['CD3_3']
1766 if flipx:
1767 dx = -data.dx
1768 data.dx = dx
1769
1770
1771
1772 try:
1773 ix = hdulist[scihead].header['CRVAL1']
1774 except KeyError:
1775 ix = 0.0
1776 ixpix = hdulist[scihead].header['CRPIX1']
1777 try:
1778 iy = hdulist[scihead].header['CRVAL2']
1779 except KeyError:
1780 iy = 0.0
1781 iypix = hdulist[scihead].header['CRPIX2']
1782 try:
1783 iz = hdulist[scihead].header['CRVAL3']
1784 except KeyError:
1785 iz = 0.0
1786 izpix = hdulist[scihead].header['CRPIX3']
1787
1788
1789 data.xs = ix - ((ixpix - 1) * data.dx)
1790 data.ys = iy - ((iypix - 1) * data.dy)
1791 data.zs = iz - ((izpix - 1) * data.dz)
1792
1793
1794 xpts = np.arange(data.xs, data.xs + (data.dx*data.nx), data.dx)
1795 data.xpts = xpts[0:data.nx]
1796 ypts = np.arange(data.ys, data.ys + (data.dy*data.ny), data.dy)
1797 data.ypts = ypts[0:data.ny]
1798 zpts = np.arange(data.zs, data.zs + (data.dz*data.nz), data.dz)
1799 data.zpts = zpts[0:data.nz]
1800
1801
1802 data.xu = xu
1803 data.yu = yu
1804 data.zu = zu
1805
1806 data.pxoff = pxoff
1807 data.pyoff = pyoff
1808
1809
1810 data.xptsoff = data.xpts - data.pxoff
1811 data.yptsoff = data.ypts - data.pyoff
1812
1813
1814
1815 hdulist.close()
1816
1817 if not(quiet):
1818 print 'loadFITS: load completed!'
1819 return data
1820
1822 """
1823 Adds a variance array to a pre-existing fitscube instance. The mask from
1824 fitscube.scidata will be applied to the variance array after loading.
1825
1826 Inputs:
1827 fitscube - The fitscube instance the variance array is to be added to
1828 vars - The variance array to be added to the fitscube instance. It must
1829 have the same shape as fitscube.scidata. Can be a nested list,
1830 a numpy array, etc etc.
1831
1832 Returns:
1833 data - A fitscube instance identical to the original fitscube, but
1834 with the new variance array
1835 """
1836
1837 if vars.shape != fitscube.scidata.shape:
1838 raise ValueError("applyVars: The variance array you're trying to add\
1839 does not have the same dimensions as the scidata!")
1840
1841 data = fitscube
1842 data.vars = np.ma.masked_array(data=vars,
1843 mask=np.ma.getmask(fitscube.scidata))
1844 return data
1845
1847 """
1848 Applies mask to the scidata and variance arrays in the fitscube.
1849
1850 Inputs:
1851 fitscube - The fitscube instance the mask is to be added to
1852 mask - The mask to be applied. The mask should take the same form as
1853 the initial DQ pixel quality array - that is, a True/1 value
1854 denotes a pixel that is bad, False/0 denotes good.
1855
1856 Returns:
1857 data - A fitscube instance identical to the original fitscube, but with
1858 the new mask applied.
1859 """
1860 data = fitscube
1861
1862
1863 mask = np.array(mask, dtype=bool)
1864
1865 data.scidata = np.ma.masked_array(data=np.ma.getdata(fitscube.scidata),
1866 mask=mask)
1867 if len(fitscube.vars != 0):
1868 data.vars = np.ma.masked_array(data=np.ma.getdata(fitscube.vars),
1869 mask=mask)
1870
1871 data.dq = mask
1872
1873 return data
1874
1875
1877 """
1878 Function to load a DQ pixel quality mask onto data that has already been
1879 loaded into a fitscube() class instance.
1880
1881 Inputs:
1882 fitscube - The pre-existing fitscube instance. If the fitscube already has
1883 a mask applied to the scidata and/or variance arrays, that mask
1884 will be overwritten.
1885 dqfile - The FITS file (possibly multi-extension) which containts the DQ
1886 pixel quality array.
1887 dqhead - The header of dqfile containing the DQ pixel quality array.
1888 Defaults to 'DQ', however, if dqfile is detected to be single-
1889 extension, it will automatically be set to 'PRIMARY'.
1890
1891 Returns:
1892 data - A fitscube instance which is identical to the input fitscube,
1893 but with the new masking applied to the scidata and vars arrays.
1894 """
1895
1896 data = fitscube
1897
1898
1899 dqhdulist = pf.open(dqfile)
1900
1901
1902 if len(dqhdulist) == 1:
1903 dqhead = 'PRIMARY'
1904
1905 data = applyMask(data, dqhdulist[dqhead].data)
1906 dqhdulist.close()
1907 return data
1908
1910 """
1911 Function to load a variance array into a fitscube() instance that has
1912 already been loaded into memory. The mask from fitscube.scidata will then
1913 be applied to the new variance array.
1914
1915 Inputs:
1916 fitscube - The pre-existing fitscube instance. If this fitscube already has
1917 a variance array in it, this array will be overwritten.
1918 varfile - The FITS file (possibly multi-extension) which contains the
1919 variance array.
1920 varhead - The header of the FITS file extension containing the variance
1921 array. Defaults to 'VAR'; however, if varfile is detected to be
1922 single-extension, it will automatically be set to 'PRIMARY'.
1923 zu - Units the sigmas are represented in (i.e. the variances will be
1924 in units of zu**2.) Options are angstrom 'A' or micron 'um'.
1925 Defaults to 'A'.
1926
1927 Returns:
1928 data - A fitscube instance identical to the input fitscube, but with
1929 the new variance data.
1930 """
1931
1932 varhdulist = pf.open(varfile)
1933
1934 if len(varhdulist) == 1:
1935 varhead = 'PRIMARY'
1936
1937
1938
1939 if fitscube.scidata.shape() != varhdulist[varhead].data.shape():
1940 raise ValueError("loadFITSvars: The shape of the variance array you're\
1941 trying to load does not match the shape of the pre-\
1942 existing data!")
1943
1944
1945
1946
1947 data = applyVars(fitscube, varhdulist[varhead].data)
1948
1949
1950 varhdulist.close()
1951 return data
1952
1953 -def loadFITSsupp(fitscube, filename, varhead='VAR', dqhead='DQ'):
1954 """
1955 Function to simultaneously load a variance array and a pixel quality mask
1956 into a pre-existing fitscube instance. This function simply calls
1957 loadFITSvars, and then applyMask, and is provided for convenience.
1958
1959 Inputs:
1960 fitscube - The pre-existing fitscube instance. If this instance already
1961 contains a variance array and/or a pixel quality mask, these
1962 will both be replaced.
1963 filename - the multi-extenstion FITS file containing the variance array and
1964 pixel quality array.
1965 varhead - The header of the extension containing the variance array.
1966 Defaults to 'VAR'.
1967 dqhead - The header of the FITS extension containing the pixel quality
1968 array. Defaults to 'DQ'.
1969
1970 Returns:
1971 data - A fitscube instance identical to the initial instance, but with
1972 the new variance array and pixel quality mask.
1973 """
1974
1975
1976 hdulist = pf.open(filename)
1977
1978
1979 data = applyVars(fitscube, hdulist[varhead].data)
1980
1981
1982 data = applyMask(data, hdulist[dqhead].data)
1983
1984
1985 hdulist.close()
1986 return data
1987
1989 """
1990 Open a saved nifscube or velfit instance, and load it into Python memory.
1991
1992 Inputs:
1993 filename - The name of the saved nifscube/velfit file to be opened.
1994
1995 Returns:
1996 velfit - The nifscube/velfit instance read out of the file.
1997 """
1998 infile = open(filename, 'rb')
1999 pninst = pickle.load(infile)
2000 print 'openPyNifs: file '+filename+' opened!'
2001 infile.close()
2002 return pninst
2003
2004
2005
2006
2007
2009 """
2010 Integreates the raw data in the nifcscube instance over the spectral
2011 dimension between the defined wavelength limits. Only unmasked pixels
2012 will be used. Integral is calculated histogram-style - that is, each data
2013 points represents a vertical column of height (data valae) and width
2014 (delta lambda).
2015
2016 Inputs:
2017 nifscube - The nifscube instance to be summed over.
2018 startw, stopw - The beginning and end wavelengths of the wavelength range to
2019 be integrated over. Given in units of m.
2020
2021 Returns:
2022 intdata - A two-dimensional masked array containing the integrated
2023 values. The integral will only be of those pixels which are
2024 unmasked in the full three-dimensional data cube.
2025 The mask will only suppress spaxels which have no unmasked
2026 value in the full three-dimensional data cube.
2027 interr - A two-dimensional masked array, containing the error on
2028 intdata related to the summation/integration process. Note
2029 that this error won't account for systematics introduced by,
2030 e.g., masked pixels within the integration range.
2031 nopoints - A two-dimensional array, containing the number of points
2032 used to compute the integrated flux density in that spaxel.
2033 """
2034
2035 intdata = np.ma.zeros(nifscube.scidata[0,:,:].shape)
2036 intvar = np.ma.zeros(nifscube.scidata[0,:,:].shape)
2037 nopoints = np.ma.zeros(nifscube.scidata[0,:,:].shape, dtype=int)
2038
2039 if nifscube.zu == 'A':
2040 startw = startw * 1.0e10
2041 stopw = stopw * 1.0e10
2042 elif nifscube.zu == 'um':
2043 startw = startw * 1.0e6
2044 stopw = stopw * 1.0e6
2045 else:
2046 raise ValueError("sumVelRange: I don't know what the z unit in this\
2047 dataset is! It's "+nifscube.zu)
2048
2049
2050
2051 list1 = np.where(nifscube.zpts > startw)
2052 list2 = np.where(nifscube.zpts < stopw)
2053
2054
2055
2056 p = max([np.amin(list1), np.amin(list2)])
2057 q = min([np.amax(list1), np.amax(list2)])
2058
2059
2060 print 'intRange: integrating over ['+str(startw)+','+str(stopw)+'] '+nifscube.zu
2061 print 'intRange: integrating over z indices '+str(p)+' to '+str(q)
2062
2063 for i in range(len(nifscube.xptsoff)):
2064 for j in range(len(nifscube.yptsoff)):
2065
2066
2067 intvals = nifscube.scidata[p:q+1, j, i]
2068 intvars = np.ma.masked_where(intvals is np.ma.masked,
2069 nifscube.vars[p:q+1, j, i])
2070 nopoints[j,i] = intvals.count()
2071
2072
2073 intvals = np.ma.filled(intvals, fill_value=0.0)
2074 intvars = np.ma.filled(intvars, fill_value=0.0)
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090 for k in range(0,len(intvals)):
2091 intdata[j,i] = intdata[j,i] + nifscube.dz * intvals[k]
2092 intvar[j,i] = intvar[j,i] + nifscube.dz**2 * intvars[k]
2093
2094
2095
2096 print 'intVelRange: Completed row '+str(i+1)+' of '+str(len(nifscube.xptsoff))
2097
2098
2099
2100
2101 if nifscube.zu == 'um':
2102 intdata = 1e4 * intdata
2103 intvar = 1e4 * 1e4 * intvar
2104
2105 interr = np.ma.sqrt(intvar)
2106 return intdata, interr, nopoints
2107
2109 """
2110 Integrates the raw data in the nifscube instance over the spectral dimension
2111 between the defined velocity limits. Only unmasked pixels will be used.
2112 Integral is calculated histogram-style - that is, each data point represents
2113 a vertical column of height (data value) and width (delta lambda).
2114
2115 Inputs:
2116 nifscube - The nifscube instance to be summed over.
2117 restwavl - The rest wavelength, given in units of m.
2118 startv, stopv - The beginning and end velocities of the velocity range to
2119 be integrated over. Given in units of km/s.
2120
2121 Returns:
2122 intdata - A two-dimensional masked array containing the integrated
2123 values. The integral will only be of those pixels which are
2124 unmasked in the full three-dimensional data cube.
2125 The mask will only suppress spaxels which have no unmasked
2126 value in the full three-dimensional data cube.
2127 interr - A two-dimensional masked array, containing the error on
2128 intdata related to the summation/integration process. Note
2129 that this error won't account for systematics introduced by,
2130 e.g., masked pixels within the integration range.
2131 nopoints - A two-dimensional array, containing the number of points
2132 used to compute the integrated flux density in that spaxel.
2133 """
2134 intdata = np.ma.zeros(nifscube.scidata[0,:,:].shape)
2135 intvar = np.ma.zeros(nifscube.scidata[0,:,:].shape)
2136 nopoints = np.ma.zeros(nifscube.scidata[0,:,:].shape, dtype=int)
2137
2138
2139
2140 startw = velToWavl(startv, restwavl)
2141 stopw = velToWavl(stopv, restwavl)
2142 if nifscube.zu == 'A':
2143 startw = startw * 1.0e10
2144 stopw = stopw * 1.0e10
2145 elif nifscube.zu == 'um':
2146 startw = startw * 1.0e6
2147 stopw = stopw * 1.0e6
2148 else:
2149 raise ValueError("sumVelRange: I don't know what the z unit in this\
2150 dataset is! It's "+nifscube.zu)
2151
2152
2153
2154 list1 = np.where(nifscube.zpts > startw)
2155 list2 = np.where(nifscube.zpts < stopw)
2156
2157 p = max([np.amin(list1), np.amin(list2)])
2158 q = min([np.amax(list1), np.amax(list2)])
2159
2160
2161 print 'intVelRange: integrating over ['+str(startv)+','+str(stopv)+'] km/s'
2162 print 'intVelRange: integrating over z indices '+str(p)+' to '+str(q)
2163
2164 for i in range(len(nifscube.xptsoff)):
2165 for j in range(len(nifscube.yptsoff)):
2166
2167
2168 intvals = nifscube.scidata[p:q+1, j, i]
2169 intvars = np.ma.masked_where(intvals is np.ma.masked,
2170 nifscube.vars[p:q+1, j, i])
2171 nopoints[j,i] = intvals.count()
2172
2173
2174 intvals = np.ma.filled(intvals, fill_value=0.0)
2175 intvars = np.ma.filled(intvars, fill_value=0.0)
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191 for k in range(0,len(intvals)):
2192 intdata[j,i] = intdata[j,i] + nifscube.dz * intvals[k]
2193 intvar[j,i] = intvar[j,i] + nifscube.dz**2 * intvars[k]
2194
2195
2196
2197 print 'intVelRange: Completed row '+str(i+1)+' of '+str(len(nifscube.xptsoff))
2198
2199
2200
2201
2202 if nifscube.zu == 'um':
2203 intdata = 1e4 * intdata
2204 intvar = 1e4 * 1e4 * intvar
2205
2206 interr = np.ma.sqrt(intvar)
2207 return intdata, interr, nopoints
2208
2210 """
2211 Sums the data in the nifscube instance over the spectral dimension between
2212 the defined velocity limits. Only unmasked pixels will be summed.
2213
2214 Inputs:
2215 nifscube - The nifscube instance to be summed over.
2216 restwavl - The rest wavelength, given in units of m.
2217 startv, stopv - The beginning and end velocities of the velocity range to
2218 be summed over. Given in units of km/s.
2219
2220 Returns:
2221 sumdata - A two-dimensional masked array containing the summed values.
2222 The sum will only be of those pixels which are unmasked in
2223 the full three-dimensional data cube.
2224 The mask will only suppress spaxels which have no unmasked
2225 value in the full three-dimensional data cube.
2226 vardata - A two-dimensional masked array containing the variances of
2227 the summed data. The variance of each spaxel is computed
2228 by adding the sqrt(variance) of each pixel, and then taking
2229 the square of the result. The same mask applied to sumdata
2230 will be applied to vardata. If there is no variance
2231 information in the input nifscube, a two-dimensional masked
2232 array of zeros in the correct shape will be returned.
2233 """
2234
2235
2236
2237 startw = velToWavl(startv, restwavl)
2238 stopw = velToWavl(stopv, restwavl)
2239 if nifscube.zu == 'A':
2240 startw = startw * 1.0e10
2241 stopw = stopw * 1.0e10
2242 elif nifscube.zu == 'um':
2243 startw = startw * 1.0e6
2244 stopw = stopw * 1.0e6
2245 else:
2246 raise ValueError("sumVelRange: I don't know what the z unit in this\
2247 dataset is! It's "+nifscube.zu)
2248
2249
2250 list1 = np.where(nifscube.zpts >= startw)
2251 list2 = np.where(nifscube.zpts <= stopw)
2252
2253 p = max([np.amin(list1), np.amin(list2)])
2254 q = min([np.amax(list1), np.amax(list2)])
2255 print 'sumVelRange: summing over ['+str(startv)+','+str(stopv)+'] km/s'
2256 print 'sumVelRange: summing over z indices '+str(p)+' to '+str(q)
2257
2258
2259 sumdata = np.ma.zeros((nifscube.ny, nifscube.nx))
2260 vardata = np.ma.zeros((nifscube.ny, nifscube.nx))
2261
2262 if np.ma.getdata(nifscube.vars) == None:
2263 for i in range(0,nifscube.nx):
2264 for j in range(0,nifscube.ny):
2265 sumdata[j,i] = nifscube.scidata[p:q+1,j,i].sum()
2266
2267
2268 if (np.ma.getmask(nifscube.scidata)[p:q+1,j,i]).all():
2269 sumdata[j,i] = np.ma.masked
2270 print 'sumVelRange: completed row '+str(i+1)+' of '+str(nifscube.nx)
2271 else:
2272 for i in range(0,nifscube.nx):
2273 for j in range(0,nifscube.ny):
2274 sumdata[j,i] = nifscube.scidata[p:q+1,j,i].sum()
2275 vardata[j,i] = nifscube.vars[p:q+1,j,i].sum()
2276
2277
2278 if (np.ma.getmask(nifscube.scidata)[p:q,j,i]).all():
2279 sumdata[j,i] = np.ma.masked
2280 vardata[j,i] = np.ma.masked
2281 print 'sumVelRange: completed row '+str(i+1)+' of '+str(nifscube.nx)
2282
2283 print 'sumVelRange: summation completed!'
2284 return sumdata, vardata
2285
2314
2315 -def fitVels(nifscube, restwavl, vmin=-1000.0, vmax=1000.0,
2316 ftestlim=1.4, ftestp=0.05, fitmax=1, xlim=None, ylim=None,
2317 vlim=None, wmin=10.0, wpercent=0.05, wperspace=0.25, wmaxcut=50.0,
2318 alim=None, fixbg=False, snthresh=5, noisepts=100):
2319 """
2320 Fits a multi-component Gaussian to every spaxel in the given NIFS data cube.
2321
2322 Inputs:
2323 nifscube - The nifscube instance to be acted uopn
2324 restwavl - The rest wavelength of interest (in m).
2325 vmin - The lowest velocity data value to be included in the
2326 fitting (in km/s). Defaults to -1000 km/s.
2327 vmax - The highest velocity data value to be included in the
2328 fitting (in km/s). Defaults to 1000 km/s.
2329 ftestlim - The upper limit allowed for the F-test, to determine whether
2330 extra Gaussian components are required or not. Defaults to
2331 1.4.
2332 ftestp - The limiting P-value for the F-test, to determine whether
2333 extra Gaussian components are required or not. Defaults to
2334 0.05.
2335 fitmax - The maximum number of Gaussian components that can be fit
2336 to any one spaxel. Defaults to 1.
2337 xlim - Either none, OR a tuple denoting the start and end x-coords
2338 of the region of interest. Defaults to None.
2339 ylim - As for xlim, but with respect to the y coordinates.
2340 vlim - A two-element list, containing min and max values for
2341 ALL fit component centroid values, in units km/s. Defaults
2342 to None (in which case, the centroid values will be forced
2343 to lie within the velocity range of interest.)
2344 wmin - Minimum allowed fitted line width. Defaults to 10.0 km/s.
2345 wpercent - Percentage of total line intensity permitted to be on the
2346 opposite side of 0 km/s to the line peak. Used to calculate
2347 maximum permitted fitted line width for each feature. Must
2348 be a float between 0.0 and 1.0. Defaults to 0.10.
2349 wperspace - Expected percentage error on initial guess line velocities.
2350 Used for calculating maximum allowed fitted line widths.
2351 Must be a float between 0.0 and 1.0. Defaults to 0.25.
2352 wmaxcut - Lower limit on maximum allowed fitted line width. Defaults
2353 to 50.0.
2354 fixbg - Boolean value, denoting whether to fix the background
2355 height of the fit to 0.0 (True) or not. Defaults to False.
2356 snthresh - Signal-to-noise threshold - if spaxel is below this
2357 threshold, fits will still be made, but fitno[j,i] will then
2358 be set to 0 afterwards. S/N is calculated as the peak flux
2359 value in the spaxel, divided by the average noise across
2360 the entire region of interest (generally speaking, the
2361 errors in any one spaxel are roughly the same in adjacent
2362 pixels.) Must be expressed as a positive number. Defaults
2363 to 5.
2364 noisepts - The number of points to be used just outside the velocity
2365 region of interest to compute an uncertainty estimate for
2366 each spaxel (to be used with signal-to-noise rejection).
2367 Points will be taken from the side of the region of interest
2368 furthest from v = 0 km/s.
2369
2370 Returns:
2371 velfit - A new velfit instance containing the fit information.
2372
2373 """
2374
2375 velfits = velfit()
2376
2377
2378
2379 print 'fitVels: Re-formatting data coordinates...'
2380 velfits.restwavl = restwavl
2381 velfits.vmin = vmin
2382 velfits.vmax = vmax
2383 velfits.ftestlim = ftestlim
2384 velfits.ftestp = ftestp
2385 velfits.fitmax = fitmax
2386
2387 velfits.xpts = nifscube.xptsoff
2388 velfits.xs = velfits.xpts[0]
2389 velfits.dx = nifscube.dx
2390 velfits.ypts = nifscube.yptsoff
2391 velfits.ys = velfits.ypts[0]
2392 velfits.dy = nifscube.dy
2393
2394 if nifscube.zu == 'A':
2395 velfits.vpts = wavlToVel(nifscube.zpts * 1.0e-10, restwavl)
2396 elif nifscube.zu == 'um':
2397 velfits.vpts = wavlToVel(nifscube.zpts * 1.0e-6, restwavl)
2398 else:
2399 raise ValueError("sumVelRange: I don't know what the z unit in this\
2400 dataset is! It's "+nifscube.zu)
2401 print velfits.vpts
2402
2403 print 'fitVels: Trimming data points to velocity region of interest'
2404
2405
2406 list1 = np.where(velfits.vpts >= vmin)
2407 list2 = np.where(velfits.vpts <= vmax)
2408 print list1
2409 print list2
2410 p = max([np.amin(list1), np.amin(list2)])
2411 q = min([np.amax(list1), np.amax(list2)])
2412
2413 velfits.vpts = velfits.vpts[p:q+1]
2414 print 'fitVels: Trimming from vmin='+str(velfits.vpts[0])+' to vmax='+\
2415 str(velfits.vpts[-1])
2416 print 'fitVels: Corresponds to indices '+str(p)+' to '+str(q)
2417
2418 velfits.vs = velfits.vpts[0]
2419 velfits.scidata = nifscube.scidata[p:q+1,:,:]
2420 velfits.vars = nifscube.vars[p:q+1,:,:]
2421 velfits.dq = nifscube.dq[p:q+1,:,:]
2422
2423
2424 velfits.uncertest = np.ma.ones(velfits.scidata[0,:,:].shape)
2425 if abs(vmin) > abs(vmax):
2426 ustart = p - noisepts
2427 ustop = p
2428 else:
2429 ustart = q + 1
2430 ustop = q + 1 + noisepts
2431 for j in range(len(nifscube.yptsoff)):
2432 for i in range(len(nifscube.xptsoff)):
2433 velfits.uncertest[j,i] = np.ma.std(nifscube.scidata[ustart:ustop
2434 ,j,i])
2435
2436
2437 if (xlim != None) or (ylim != None):
2438 if xlim == None:
2439 print 'fitVels: Trimming data to y = '+str(ylim)
2440 xlim = [velfits.xpts[0]-1.0, velfits.xpts[-1]+1]
2441 elif ylim == None:
2442 print 'fitVels: Trimming data to x = '+str(xlim)
2443 ylim = [velfits.ypts[0]-1.0, velfits.ypts[-1]+1]
2444 else:
2445 print 'fitVels: Trimming data to x = '+str(xlim)\
2446 +', y = '+str(ylim)
2447 velfits.scidata = trimPts3D(velfits.scidata, velfits.xpts, velfits.ypts,
2448 xlim[0],xlim[1],ylim[0],ylim[1],incall=True,
2449 dx=velfits.dx, dy=velfits.dy, quiet=True)[0]
2450 velfits.uncertest = trimPts(velfits.uncertest, velfits.xpts,
2451 velfits.ypts,xlim[0],xlim[1],ylim[0],
2452 ylim[1],incall=True,dx=velfits.dx,
2453 dy=velfits.dy)[0]
2454 velfits.vars = trimPts3D(velfits.vars, velfits.xpts,
2455 velfits.ypts, xlim[0], xlim[1], ylim[0],
2456 ylim[1], incall=True, dx=velfits.dx,
2457 dy=velfits.dy, quiet=True)[0]
2458 velfits.dq, velfits.xpts, velfits.ypts = trimPts3D(velfits.dq,
2459 velfits.xpts,
2460 velfits.ypts,
2461 xlim[0], xlim[1],
2462 ylim[0], ylim[1],
2463 incall=True,
2464 dx=velfits.dx,
2465 dy=velfits.dy,
2466 quiet=True)
2467 velfits.xs = velfits.xpts[0]
2468 velfits.ys = velfits.ypts[0]
2469
2470
2471
2472
2473 velfits.fitno = np.ma.empty_like(velfits.scidata[0,:,:], dtype=int)
2474 velfits.fitorig = np.ma.empty_like(velfits.scidata[0,:,:], dtype=int)
2475
2476
2477
2478
2479
2480
2481 velfits.fits = np.ma.zeros((velfits.scidata[0,:,:].shape + (velfits.fitmax,
2482 velfits.fitmax*3 + 1)), dtype=np.float64)
2483
2484 for k in range(2, 3*velfits.fitmax, 3):
2485 velfits.fits[:,:,:,k] = 1.0
2486 velfits.perror = np.ma.zeros((velfits.scidata[0,:,:].shape + (velfits.fitmax,
2487 velfits.fitmax*3 + 1)), dtype=np.float64)
2488 velfits.dof = np.ma.zeros((velfits.scidata[0,:,:].shape + (velfits.fitmax ,)), dtype=int)
2489 velfits.chisq = np.ma.zeros((velfits.scidata[0,:,:].shape + (velfits.fitmax ,)), dtype=np.float64)
2490 velfits.fitstat = np.ma.zeros((velfits.scidata[0,:,:].shape + (velfits.fitmax ,)), dtype=int)
2491 velfits.fluxes = np.ma.zeros((velfits.scidata[0,:,:].shape + (velfits.fitmax ,)), dtype=np.float64)
2492 velfits.fluxerr = np.ma.zeros((velfits.scidata[0,:,:].shape + (velfits.fitmax ,)), dtype=np.float64)
2493
2494
2495
2496
2497 print 'fitVels: Generating velocity fits...'
2498 for i in range(0,len(velfits.xpts)):
2499 for j in range(0, len(velfits.ypts)):
2500 print 'fitVels: Fitting spaxel ('+str(i)+','+str(j)+')'
2501 velfits.fitVelSpax(i, j, comps=fitmax, quiet=True, vlim=vlim,
2502 wmin=wmin, alim=alim, fixbg=fixbg,
2503 snthresh=snthresh, wpercent = wpercent,
2504 wperspace = wperspace, wmaxcut=wmaxcut)
2505 print 'fitVels: Completed fitting for row '+str(i+1)+' of '\
2506 +str(len(velfits.xpts))
2507
2508
2509
2510 velfits.fitsorig = velfits.fits.copy()
2511
2512
2513 print 'fitVels: Computing summed data for each spaxel...'
2514 velfits.sumCube()
2515 print 'fitVels: Computing fitted flux for each spaxel...'
2516 velfits.sumFitFlux()
2517
2518 return velfits
2519
2520
2521
2522
2523
2524
2525
2526
2527 c = 2.99792458e8
2528
2529 fwhm = 2.0 * np.sqrt(math.log(2.0))
2530
2532 """Converts a wavelength to a radial velocity.
2533
2534 Inputs:
2535 wavl - Wavelength to be converted (units m).
2536 restwavl - Rest wavelength (units m).
2537 wavs - Wave speed (units m s^-1). Defaults to speed of light c = 2.9979e8 m s^-1.
2538
2539 Returns:
2540 radvel - Corresponding radial velocity (units km s^-1)."""
2541
2542 radvel = ((wavl/restwavl) - 1.0) * wavs
2543 return (radvel/1000.0)
2544
2546 """Converts a radial velocity to a wavelength.
2547
2548 Inputs:
2549 radvel - Radial velocity (units km s^-1).
2550 restwavl - Rest wavelength (units m).
2551 wavs - Wave speed (units m s^-1). Defaults to speed of light c = 2.9979e8 m s^-1.
2552
2553 Returns:
2554 wavl - Corresponding wavelength (units m)."""
2555
2556 radvel = radvel * 1000.0
2557 wavl = ((wavs + radvel) / wavs) * restwavl
2558 return wavl
2559
2561 """
2562 Masks all pixels in the cube with the specified x- and y-index.
2563
2564 Inputs:
2565 cube - Cube to be acted upon
2566 xind - The x-index of pixels to be masked
2567 yind - The y-index of pixels to be masked
2568
2569 Returns:
2570 Nil.
2571 """
2572
2573 if (xind >= cube.nx) or (yind > cube.ny):
2574 raise ValueError('maskPix: Given index is out of bounds!')
2575
2576 print 'maskPix: Masking pixel...'
2577 for k in range(0,cube.nz):
2578 np.ma.masked_where(cube.scidata==cube.scidata[k,yind,xind],
2579 cube.scidata)
2580 np.ma.masked_where(cube.scidata==cube.scidata[k,yind,xind],
2581 cube.vars)
2582 print 'Completed pixel '+str(k)+' of '+str(cube.nz)
2583
2584 print 'maskPix: Pixel ('+str(xind)+','+str(yind)+') masked!'
2585
2587 """
2588 Converts the arrays of pixel centroids xpts and ypts in the passed fitscube
2589 into arrays of pixel corner values. Either returns the two arrays, or a
2590 meshgrid of them. Uses the offset position arrays
2591
2592 Inputs:
2593 fitscube - The fitscube instance to be acted upn.
2594 mesh - Denotes whether to return lists of pixel corner x- and y-
2595 coordinates (False), or a meshgrid of the two for plotting
2596 purposes (True).
2597
2598 Returns:
2599 EITHER
2600 xcors, ycors - The lists of pixel corner x- and y-coordinates,
2601 respectively.
2602 OR
2603 XX, YY - The result of meshgrid(xcors, ycors)
2604 """
2605
2606 xcors = fitscube.xptsoff - (fitscube.dx/2.0)
2607 xcors = np.append(xcors, xcors[-1]+fitscube.dx)
2608 ycors = fitscube.yptsoff - (fitscube.dy/2.0)
2609 ycors = np.append(ycors, ycors[-1]+fitscube.dy)
2610
2611 if mesh:
2612 XX, YY = np.meshgrid(xcors, ycors)
2613 return XX, YY
2614 else:
2615 return xcors, ycors
2616
2618 xcors = xpts - dx/2.0
2619 xcors = np.append(xcors, xcors[-1]+dx)
2620 ycors = ypts - dy/2.0
2621 ycors = np.append(ycors, ycors[-1]+dy)
2622
2623 if mesh:
2624 X, Y = np.meshgrid(xcors, ycors)
2625 return X, Y
2626 else:
2627 return xcors, ycors
2628
2629 -def gridRad(data, xpts, ypts, vars=None, xcen = 0.0, ycen=0.0):
2630 """
2631 Re-grids a two-dimensional data array (with coordinates given in x- and y-
2632 space) into a one-dimensional data array with coordinates given in radius.
2633 The radius is computed from the given x and y center coordinates.
2634
2635 Inputs:
2636 data - The two-dimensional (potentially masked) data array to be
2637 re-gridded. Must have shape(len(ypts), len(xpts)).
2638 xpts - The list of x-coordinates corresponding to the data.
2639 ypts - The list of y-coordinates corresponding to the data.
2640 vars - The array of data variances corresponding to the data. Defaults to
2641 None. Must be the same shape as data.
2642 xcen - The center x-coordinate to be used for computing the radius.
2643 Defaults to 0.0.
2644 ycen - The center x-coordinate to be used for computing the radius.
2645 Defaults to 0.0.
2646
2647 Returns:
2648 datarad - A one-dimensional array (potentially masked), of length
2649 len(xpts)*len(ypts), containing the input data re-gridded into
2650 radius space.
2651 varsrad - A one-dimensional array (potentially masked) containing the
2652 variance information of the input data re-gridded to radius
2653 space. If no variances were passed when the function was called,
2654 an array of zeros of the appropriate size will be returned.
2655 rads - A one-dimensional array containing the radius coordinates of the
2656 re-gridded data. Specified such that rads[i] is the radial
2657 coordinates of datarad[i]. Points in rad will almost certainly
2658 NOT be sorted in any particular order.
2659 """
2660
2661 usevars = not(vars == None)
2662
2663
2664 if (vars != None) and (data.shape != vars.shape):
2665 raise ValueError('Cannot process data and variance arrays of different\
2666 shapes: '+str(data.shape)+ 'and '+str(vars.shape))
2667 elif (data.shape[0] != len(ypts)):
2668 raise ValueError('Mismatched data array and y-coordinate array.')
2669 elif (data.shape[1] != len(xpts)):
2670 raise ValueError('Mismatched data array and y-coordinate array.')
2671
2672
2673 datarad = np.ma.zeros(np.ma.size(data))
2674 varsrad = np.ma.zeros(np.ma.size(data))
2675 rads = np.zeros(len(xpts)*len(ypts))
2676
2677
2678 print 'gridRad: Re-gridding data...'
2679 for i in range(0, len(xpts)):
2680 for j in range(0, len(ypts)):
2681
2682 rads[i*len(ypts) + j] = np.sqrt( (xpts[i] - xcen)**2
2683 +(ypts[j] - ycen)**2 )
2684
2685 datarad[i*len(ypts) + j] = data[j][i]
2686 if usevars:
2687 varsrad[i*len(ypts) + j] = vars[j][i]
2688 print 'gridRad: Completed row '+str(i+1)+' of '+str(len(xpts))
2689
2690
2691
2692
2693 return datarad, varsrad, rads
2694
2695
2696 -def trimPts(data, xpts, ypts, xmin, xmax, ymin, ymax,
2697 incall=False, dx=0.0, dy=0.0):
2698 """
2699 Function that trims a two-dimensional data array so that it only retains
2700 points that are within the ranges [xmin,xmax] and [ymin,max].
2701
2702 Inputs:
2703 data - The two-dimensional data array to be acted upon
2704 xpts, ypts - The list of x- and y-coordinates corresponding to the data
2705 array
2706 xmin, xmax - The minimum and maximum x-values to be left in the array
2707 ymin, ymax - The minimum and maximum y-values to be left in the array
2708 incall - Boolean value describing whether points are included based
2709 on their centroid coordinate (False), or whether any point
2710 in the pixel is within the specified ranges (True). Defaults
2711 to False.
2712 dx, dy - The pixel sizes of the data array.
2713
2714 Returns:
2715 trimdata - The trimmed down data array
2716 trimxpts - The pixel centroid x-coordinates corresponding to trimdata
2717 trimypts - The pixel centroid y-coordinates corresponding to trimdata
2718 """
2719
2720 if (data.shape) != (len(ypts), len(xpts)):
2721 raise ValueError('Data array shape '+str((data.shape))+\
2722 'does not match the coordinate array sizes '+\
2723 str(len(ypts))+', '+str(len(xpts)))
2724
2725 if incall:
2726 filterx = np.logical_and(xpts >= (xmin - dx/2.0), xpts <= (xmax + dx/2.0))
2727 filtery = np.logical_and(ypts >= (ymin - dy/2.0), ypts <= (ymax + dy/2.0))
2728 else:
2729 filterx = np.logical_and(xpts >= (xmin), xpts <= (xmax))
2730 filtery = np.logical_and(ypts >= (ymin), ypts <= (ymax))
2731
2732 print 'trimPts: Data trimmed using following limits:'
2733 print 'xmin:'+str(xmin)
2734 print 'xmax:'+str(xmax)
2735 print 'ymin:'+str(ymin)
2736 print 'ymax:'+str(ymax)
2737
2738
2739
2740 interdata = data[filtery, :]
2741 trimdata = interdata[:, filterx]
2742 trimxpts = xpts[filterx]
2743 trimypts = ypts[filtery]
2744
2745 return trimdata, trimxpts, trimypts
2746
2747 -def trimPts3D(data, xpts, ypts, xmin, xmax, ymin, ymax,
2748 incall=False, dx=0.0, dy=0.0, quiet=False):
2749 """
2750 Function that trims a three-dimensional data array so that it only retains
2751 points that are within the ranges [xmin,xmax] and [ymin,max].
2752
2753 Inputs:
2754 data - The two-dimensional data array to be acted upon
2755 xpts, ypts - The list of x- and y-coordinates corresponding to the data
2756 array
2757 xmin, xmax - The minimum and maximum x-values to be left in the array
2758 ymin, ymax - The minimum and maximum y-values to be left in the array
2759 incall - Boolean value describing whether points are included based
2760 on their centroid coordinate (False), or whether any point
2761 in the pixel is within the specified ranges (True). Defaults
2762 to False.
2763 dx, dy - The pixel sizes of the data array.
2764
2765 Returns:
2766 trimdata - The trimmed down data array
2767 trimxpts - The pixel centroid x-coordinates corresponding to trimdata
2768 trimypts - The pixel centroid y-coordinates corresponding to trimdata
2769 """
2770
2771 if data.shape[1:] != (len(ypts), len(xpts)):
2772 raise ValueError('Data array shape '+str(data.shape[1:])+\
2773 'does not match the coordinate array sizes'+\
2774 str(len(ypts))+', '+str(len(xpts)))
2775
2776 if incall:
2777 filterx = np.logical_and(xpts >= (xmin - dx/2.0), xpts <= (xmax + dx/2.0))
2778 filtery = np.logical_and(ypts >= (ymin - dy/2.0), ypts <= (ymax + dy/2.0))
2779 else:
2780 filterx = np.logical_and(xpts >= (xmin), xpts <= (xmax))
2781 filtery = np.logical_and(ypts >= (ymin), ypts <= (ymax))
2782
2783 if not(quiet):
2784 print 'trimPts3D: Data trimmed using following limits:'
2785 print 'xmin:'+str(xmin)
2786 print 'xmax:'+str(xmax)
2787 print 'ymin:'+str(ymin)
2788 print 'ymax:'+str(ymax)
2789
2790
2791
2792 interdata = data[:, filtery, :]
2793
2794 trimdata = interdata[:, :, filterx]
2795
2796 trimxpts = xpts[filterx]
2797 trimypts = ypts[filtery]
2798
2799 return trimdata, trimxpts, trimypts
2800
2802 """
2803 Function flattens the input list. Input list can have arbitrary, changing
2804 depth throughout.
2805
2806 Inputs:
2807 x - The list to be flattened.
2808
2809 Returns:
2810 The flattened version of x
2811 """
2812 if isinstance(x, collections.Iterable):
2813 return [a for i in x for a in flatten(i)]
2814 else:
2815 return [x]
2816
2817 -def lineFluxKms(amp, width, restwavl, amperr=0.0, widerr=0.0):
2818 '''
2819 Computes the total flux incident from a Gaussian line profile, which has
2820 base units of km/s. Assumes the line profile has a background height of zero
2821 (otherwise, the integral over +-infinity would result in infinity).
2822
2823 Note that it has not been necessary to fully transform the Gaussian function
2824 for integration over wavelength space. The distribution is actually already
2825 *in* wavelength space, being expressed in units of, say, erg/cm^2/s/lambda.
2826 However, given that its defined over velocity space (i.e. line velocity
2827 is the x-unit), we simply have to convert the width of the Gaussian into
2828 wavelength via sigma_lambda = sigma*restwavl / c.
2829
2830 Inputs:
2831 amp - The amplitude of the line profile. Can be passed in any units
2832 desired - however, the 'spectral unit' must match the unit which
2833 restwavl is expressed in. For example, if the line peak is
2834 expressedin units of erg/cm^2/s/A, the rest wavelength MUST be
2835 expressed in A to get the correct result.
2836 width - Line profile width, in km/s.
2837 restwavl - Rest wavelength of the line being investigated. Must be
2838 expressed in similar units to the line profile amplitude (see
2839 above).
2840 amperr - The uncertainty in the amplitude of the line profile, in the
2841 same units as amp. Defaults to 0.0.
2842 widerr - The uncertainty in the width of the line profile, in units km/s.
2843 Defaults to 0.0.
2844
2845 Returns:
2846 lineflux - The computed line flux. The units will mirror those used for
2847 the line profile amplitude, less the 'spectral unit' - e.g.
2848 if amplitude is passed as erg/cm^2/s/A, the flux returned will
2849 be in units of erg/cm^2/s.
2850 lineferr - The error in the computed line flux, in the same units as
2851 lineflux. Will come back as zero if amperr = widerr = 0.0.
2852 '''
2853
2854
2855
2856
2857
2858
2859
2860 lineflux = (np.sqrt(np.pi) * amp * restwavl * width) / (c / 10.0**3)
2861
2862
2863
2864 q = (amperr / amp)**2 + (widerr / width)**2
2865 lineferr = lineflux * np.sqrt(q)
2866
2867 return lineflux, lineferr
2868
2870 """
2871 Recomputes the fluxes and fluxerrs arrays of the designated velfit interest.
2872 Designed for use when the velfit instance has been modified/damaged,
2873 or on velfit instances made before the sqrt(2) discrepancy in the original
2874 calculation was identified (26 Sep 12).
2875
2876 Inputs:
2877 velfits - The velfit instance to be updated. Fluxes are re-calculated
2878 based on the velfit.fits array.
2879 wavlu - The wavelength units in the flux density values. Accepted values
2880 are 'A' and 'um'. Defaults to 'A'.
2881
2882 Returns:
2883 Nil. The flux and fluxerr arrays of the velfit instance are updated in-situ.
2884 """
2885
2886
2887 velfits.fluxes = np.ma.zeros((velfits.scidata[0,:,:].shape
2888 + (velfits.fitmax ,)), dtype=np.float64)
2889 velfits.fluxerr = np.ma.zeros((velfits.scidata[0,:,:].shape
2890 + (velfits.fitmax ,)), dtype=np.float64)
2891 if wavlu == 'A':
2892 restwavl = velfits.restwavl * 1.0e10
2893 elif wavlu == 'um':
2894 restwavl = velfits.restwavl * 1.0e6
2895 else:
2896 raise ValueError('Unknown wavelength unit '+str(wavlu))
2897
2898 for i in range(len(velfits.xpts)):
2899 for j in range(len(velfits.ypts)):
2900 for f in range(velfits.fitno[j,i]):
2901 velfits.fluxes[j,i,f], velfits.fluxerr[j,i,f] = lineFluxKms(velfits.fits[j,i,velfits.fitno[j,i]-1,1+(f*3)],
2902 velfits.fits[j,i,velfits.fitno[j,i]-1,3+(f*3)],
2903 restwavl,
2904 velfits.perror[j,i,velfits.fitno[j,i]-1,1+(f*3)],
2905 velfits.perror[j,i,velfits.fitno[j,i]-1,3+(f*3)])
2906 print 'regenFlux: Completed row '+str(j+1)+' of '+str(len(velfits.ypts))
2907
2908
2910 """
2911 Computes the maximum allowed width for a Gaussian centered on centroid, such
2912 that at most intpercent of the total Gaussian intensity is beyond x=0
2913 (relative to which side of x=0 the line centroid is on). Assumes a
2914 Gaussian of the form A.exp(-(x-x0)^2/s^2) (note no factor of 2 in the
2915 exponent denominator).
2916
2917 Inputs:
2918 centroid - The value of the line centroid.
2919 intpercent - The maximum allowable percentage of the Gaussian function
2920 intensity (i.e. the area under the Gaussian) can be on the
2921 opposite side of x=0 relative to the centroid position. Must
2922 be a float between 0 and 1.
2923 errfact - The percentage error expected in the position of the line
2924 centroid position, which will allow for some 'wiggle room'
2925 in the returned width limit. Must be a float between 0 (no
2926 error) and 1 (error is comparable to the centroid value
2927 itself).
2928
2929 Returns:
2930 widlimit - Upper limit on Gaussian sigma (NOT FWHM). Will be returned
2931 in the same units as centroid.
2932 """
2933 intpercent = float(intpercent)
2934 errfact = float(errfact)
2935
2936 if (intpercent < 0.0) or (intpercent > 1.0):
2937 raise ValueError('maxWidthAllowed: Illegal value '+str(intpercent)+' for intpercent passed!')
2938 if (errfact < 0.0) or (errfact > 1.0):
2939 raise ValueError('maxWidthAllowed: Illegal value '+str(errfact)+' for errfact passed!')
2940
2941
2942
2943 errfact = 1.0 + errfact
2944
2945 widlimit = abs(centroid*errfact) / np.sqrt(math.log(1.0/intpercent))
2946 return widlimit
2947
2949 """
2950 Find the spaxel index corresponding to the given z (wavelength) value.
2951
2952 Inputs:
2953 instance - The nifscube/velfit instance to be used.
2954 z - The wavelength value that we wish to find the spaxel index
2955 for. Specified in m.
2956
2957 Returns:
2958 j - The spaxel y-index corresponding to the passed y-coordinate.
2959 """
2960
2961 zpts = instance.zpts
2962 if instance.zu == "A":
2963 z = z * 1.0e10
2964 elif instance.zu == "um":
2965 z = z * 1.0e6
2966
2967
2968 if (z < (zpts[0] - instance.dz/2.0)) or (z > (zpts[-1] + instance.dz/2.0)):
2969 raise ValueError('The wavelength coordinate '+str(z)+' is outside the range of this fitscube/velfit instance!')
2970
2971
2972
2973 zbases = zpts - (instance.dz/2.0)
2974
2975
2976 coords = np.where(zbases < z)[0]
2977 k = coords[-1]
2978 return k
2979
2981 """
2982 Find the spaxel index corresponding to the given y value.
2983
2984 Inputs:
2985 instance - The nifscube/velfit instance to be used.
2986 y - The y-coordinate that we wish to find the spaxel index for
2987 useoffsets - Specifies whether to use the raw ypts array in the instance,
2988 or the offset values (fitscube instance only). Defaults to
2989 False. An error will be thrown if useoffsets=True is passed
2990 in concert with a velfit instance.
2991
2992 Returns:
2993 j - The spaxel y-index corresponding to the passed y-coordinate.
2994 """
2995
2996 if useoffsets:
2997 ypts = instance.yptsoff
2998 else:
2999 ypts = instance.ypts
3000
3001
3002 if (y < (ypts[0] - instance.dy/2.0)) or (y > (ypts[-1] + instance.dy/2.0)):
3003 raise ValueError('The y-coordinate '+str(y)+' is outside the range of this fitscube/velfit instance!')
3004
3005
3006
3007 ybases = ypts - (instance.dy/2.0)
3008
3009
3010 coords = np.where(ybases < y)[0]
3011 j = coords[-1]
3012 return j
3013
3015 """
3016 Find the spaxel index corresponding to the given x value.
3017
3018 Inputs:
3019 instance - The nifscube/velfit instance to be used.
3020 x - The x-coordinate that we wish to find the spaxel index for
3021 useoffsets - Specifies whether to use the raw ypts array in the instance,
3022 or the offset values (fitscube instance only). Defaults to
3023 False. An error will be thrown if useoffsets=True is passed
3024 in concert with a velfit instance.
3025
3026 Returns:
3027 i - The spaxel x-index corresponding to the passed x-coordinate.
3028 """
3029
3030 if useoffsets:
3031 xpts = instance.xptsoff
3032 else:
3033 xpts = instance.xpts
3034
3035
3036 if (x < (xpts[0] - instance.dx/2.0)) or (x > (xpts[-1] + instance.dx/2.0)):
3037 raise ValueError('The x-coordinate '+str(x)+' is outside the range of this fitscube/velfit instance!')
3038
3039
3040
3041 xbases = xpts - (instance.dx/2.0)
3042
3043
3044 coords = np.where(xbases < x)[0]
3045 i = coords[-1]
3046 return i
3047
3048 -def findIndex(instance, x, y, useoffsets=False):
3049 """
3050 Find the spaxel index corresponding to the given x & y values.
3051
3052 Inputs:
3053 instance - The nifscube/velfit instance to be used.
3054 x - The x-coordinate that we wish to find the spaxel index for
3055 y - The xycoordinate that we wish to find the spaxel index for
3056 useoffsets - Specifies whether to use the raw ypts array in the instance,
3057 or the offset values (fitscube instance only). Defaults to
3058 False. An error will be thrown if useoffsets=True is passed
3059 in concert with a velfit instance.
3060
3061 Returns:
3062 i - The spaxel x-index corresponding to the passed x-coordinate.
3063 j - The spaxel y-index corresponding to the passed y-coordinate.
3064 """
3065
3066 i = findXIndex(instance, x, useoffsets=useoffsets)
3067 j = findYIndex(instance, y, useoffsets=useoffsets)
3068
3069 return i, j
3070
3072 """
3073 Computes the distance between two points.
3074
3075 Inputs:
3076 x1, y1 - First point
3077 x2, y2 - Second point
3078
3079 Returns:
3080 dist - The distance between the two points.
3081 """
3082 xdiff = x2 - x1
3083 ydiff = y2 - y1
3084 dist = np.sqrt((xdiff*xdiff) + (ydiff*ydiff))
3085 return dist
3086
3087
3088
3089
3090
3092 """
3093 Compute the electron density equivalent to the 1.553/1.644 um [Fe II]
3094 flux ratio, as described by Pesenti et al. 2003. Result corresponds to a
3095 gas temperature of 10,000 K. Formula taken from Agra-Amboage et al. 2011.
3096
3097 Inputs:
3098 R - The 1.533/1.644 um [Fe II] line flux ratio.
3099
3100 Returns:
3101 ne - The computed electron number density. Returns units of cm^{-3}.
3102 """
3103
3104 numer = R - 0.035
3105 denom = 0.395 - R
3106 ne = 1.2e4 * (numer/denom)
3107 if (R < 0.035):
3108 return 1e2
3109 elif (R > 0.395):
3110 return 1e6
3111 else:
3112 return(ne)
3113
3114 -def lpAnderson(winf, vp, vphi, Mstar, doerrs=False, winferr=0.0, vperr=0.0,
3115 vphierr=0.0, Mstarerr=0.0):
3116 """
3117 Returns the footprint radius (and associated uncertainty, if required) of
3118 an MHD wind, as per Anderson et al. 2004.
3119
3120 Inputs:
3121 winf - The radius (i.e. distance from the outflow center) at which the
3122 poloidal and rotational velocity measurements are made, in units
3123 of AU.
3124 vp - The poloidal (i.e. downstream) outflow velocity at the point
3125 specified, in units of km s^-1.
3126 vphi - The rotational outflow velocity at the point specified, in units
3127 of km s^-1.
3128 Mstar - The mass of the central star, in units of M_sun.
3129 doerrs - A Boolean value, specifying whether to compute and return the
3130 uncertainty on the footprint radius (True) or not (False).
3131 Defaults to False.
3132 winferr - The uncertainty of the radius winf, in units of AU. Defaults to
3133 0.0.
3134 vperr - The uncertainty of poloidal velocity vp, in units km s^-1.
3135 Defaults to 0.0.
3136 vphierr - The uncertainty on the rotational velocity vphi, in units
3137 km s^-1. Defaults to 0.0.
3138 Mstar - The uncertainty on the stellar mass Mstar, in units M_Sun.
3139 Defaults to 0.0.
3140 """
3141
3142
3143
3144 winf = abs(winf)
3145 vp = abs(vp)
3146 vphi = abs(vphi)
3147 Mstar = abs(Mstar)
3148
3149
3150 w0 = 0.7*((winf/10.0)**(2.0/3.0) * (vphi/10.0)**(2.0/3.0)
3151 * (vp/100.0)**(-4.0/3.0) * (Mstar)**(1.0/3.0))
3152
3153
3154
3155 if doerrs:
3156 w0err = w0 * ((2.0/3.0)*(1.0/winf)*winferr + (2.0/3.0)*(1/vphi)*vphierr
3157 + (-4.0/3.0)*(1.0/vp)*vperr
3158 + (1.0/3.0)*(1.0/Mstar)*Mstarerr)
3159 return w0, w0err
3160 else:
3161 return w0
3162
3163
3164
3166 """
3167 Generates a colormap based on an underlying cubic function (i.e. plotting
3168 the variation in color as a function of x from 0 to 1 would appear to be
3169 cubic). Note this function is designed to deal with smoothly varying
3170 colormaps (i.e. those for which every (x,y0,y1) colormap tuple has y0 = y1)
3171 - I have no idea what will happen if a segmented colormap is passed.
3172
3173 Inputs:
3174 startcolor, - The start and stop colours of the colormap. Should be given
3175 stopcolor as RGB 3-tuples.
3176 interval - The number of segments to break the colormap into.
3177 More segments will make a more 'accurate' conversion away
3178 from the linear colormap, but will slow down computation of
3179 the new colormap and make it take up more memory. It is
3180 recommended that trial-and-error may be useful for
3181 determining an adequate value for interval. Defaults to 10.
3182
3183 Returns:
3184 cmap_cubic - The cubic-based colormap. This will be a smoothly varying
3185 colormap, regardless of the colormap that's put in.
3186 """
3187
3188 def cubic(x):
3189 c = 4*(x-0.5)**3 + 0.5
3190 return c
3191
3192 colors = ['red', 'green', 'blue']
3193 cmap_cubicdata = {}
3194 for color in range(0,len(colors)):
3195 colorstart = float(startcolor[color])
3196 colorstop = float(stopcolor[color])
3197 cmap_cubicdata[colors[color]] = []
3198 for i in range(0, interval+1):
3199 x = float(i)/float(interval)
3200 colorval = cubic(x) * (colorstop - colorstart) + colorstart
3201 cmap_cubicdata[colors[color]].append((x, colorval, colorval))
3202
3203 return matplotlib.colors.LinearSegmentedColormap(matplotlib.cm.autumn, cmap_cubicdata)
3204
3206 """
3207 Generates a colormap based on an underlying quintic function (i.e. plotting
3208 the variation in color as a function of x from 0 to 1 would appear to be
3209 cubic). Note this function is designed to deal with smoothly varying
3210 colormaps (i.e. those for which every (x,y0,y1) colormap tuple has y0 = y1)
3211 - I have no idea what will happen if a segmented colormap is passed.
3212
3213 Inputs:
3214 startcolor, - The start and stop colours of the colormap. Should be given
3215 stopcolor as RGB 3-tuples.
3216 interval - The number of segments to break the colormap into.
3217 More segments will make a more 'accurate' conversion away
3218 from the linear colormap, but will slow down computation of
3219 the new colormap and make it take up more memory. It is
3220 recommended that trial-and-error may be useful for
3221 determining an adequate value for interval. Defaults to 10.
3222
3223 Returns:
3224 cmap_cubic - The cubic-based colormap. This will be a smoothly varying
3225 colormap, regardless of the colormap that's put in.
3226 """
3227
3228 def quintic(x):
3229 c = 16*(x-0.5)**5 + 0.5
3230 return c
3231
3232 colors = ['red', 'green', 'blue']
3233 cmap_cubicdata = {}
3234 for color in range(0,len(colors)):
3235 colorstart = float(startcolor[color])
3236 colorstop = float(stopcolor[color])
3237 cmap_cubicdata[colors[color]] = []
3238 for i in range(0, interval+1):
3239 x = float(i)/float(interval)
3240 colorval = quintic(x) * (colorstop - colorstart) + colorstart
3241 cmap_cubicdata[colors[color]].append((x, colorval, colorval))
3242
3243 return matplotlib.colors.LinearSegmentedColormap(matplotlib.cm.autumn, cmap_cubicdata)
3244
3245
3247 """
3248 Generates a colormap based on an underlying heptic function (i.e. plotting
3249 the variation in color as a function of x from 0 to 1 would appear to be
3250 cubic). Note this function is designed to deal with smoothly varying
3251 colormaps (i.e. those for which every (x,y0,y1) colormap tuple has y0 = y1)
3252 - I have no idea what will happen if a segmented colormap is passed.
3253
3254 Inputs:
3255 startcolor, - The start and stop colours of the colormap. Should be given
3256 stopcolor as RGB 3-tuples.
3257 interval - The number of segments to break the colormap into.
3258 More segments will make a more 'accurate' conversion away
3259 from the linear colormap, but will slow down computation of
3260 the new colormap and make it take up more memory. It is
3261 recommended that trial-and-error may be useful for
3262 determining an adequate value for interval. Defaults to 10.
3263
3264 Returns:
3265 cmap_cubic - The cubic-based colormap. This will be a smoothly varying
3266 colormap, regardless of the colormap that's put in.
3267 """
3268
3269 def heptic(x):
3270 c = 16*(x-0.5)**5 + 0.5
3271 return c
3272
3273 colors = ['red', 'green', 'blue']
3274 cmap_cubicdata = {}
3275 for color in range(0,len(colors)):
3276 colorstart = float(startcolor[color])
3277 colorstop = float(stopcolor[color])
3278 cmap_cubicdata[colors[color]] = []
3279 for i in range(0, interval+1):
3280 x = float(i)/float(interval)
3281 colorval = heptic(x) * (colorstop - colorstart) + colorstart
3282 cmap_cubicdata[colors[color]].append((x, colorval, colorval))
3283
3284 return matplotlib.colors.LinearSegmentedColormap(matplotlib.cm.autumn, cmap_cubicdata)
3285
3287 """
3288 Generates a linearly interpolated colormap, with a 'break' in the linear
3289 profile. The break occurs breakpc along the colormap, at a height
3290 of breakheight.
3291
3292 Inputs:
3293 startcolor - The starting color of the colormap. Should be an RGB
3294 3-tuple.
3295 stopcolor - The end color of the colormap. Should be an RGB 3-tuple.
3296 breakpc - The percentage along the x-axis that the linear break
3297 occurs. Must be between 0 and 1.
3298 breakheight - The 'height' along the color axis where the break occurs.
3299 Must be between 0 and 1.
3300
3301 Returns:
3302 break_cmap - The colormap instance.
3303 """
3304
3305 breakpc = float(breakpc)
3306 if (breakpc < 0.0) or (breakpc > 1.0):
3307 raise ValueError('linearBreakColormap: Invalid breakpc value')
3308 breakheight = float(breakheight)
3309 if (breakheight < 0.0) or (breakheight > 1.0):
3310 raise ValueError('linearBreakColormap: Invalid breakheight value')
3311
3312 colors = ['red', 'green', 'blue']
3313 cmap_breakdata = {}
3314 for color in range(0,len(colors)):
3315 colorstart = float(startcolor[color])
3316 colorstop = float(stopcolor[color])
3317 cmap_breakdata[colors[color]] = [(0.0, colorstart, colorstart),
3318 (breakpc, breakheight*(colorstop
3319 -colorstart),
3320 breakheight*(colorstop-colorstart)),
3321 (1.0, colorstop, colorstop)]
3322
3323 return matplotlib.colors.LinearSegmentedColormap(matplotlib.cm.autumn,
3324 cmap_breakdata)
3325
3327 """
3328 Returns a modified version of the RdYlBu colormap, where the center of the
3329 'yellow' appears at xpt.
3330
3331 Inputs:
3332 xpt - The position in the colormap where yellow appears. Must be between
3333 0 and 1.
3334
3335 Returns:
3336 modRdYlBu - The modified colormap.
3337 """
3338
3339 xpt = float(xpt)
3340 if (xpt < 0.0) or (xpt > 1.0):
3341 raise ValueError('modRdYlBu: invalid xpoint passed')
3342
3343 colordict = {'red': [(0.0, 0.6470588445663452, 0.6470588445663452),
3344 (xpt, 1.0, 1.0),
3345 (1.0, 0.1921568661928177, 0.1921568661928177)],
3346 'green': [(0.0, 0.0, 0.0),
3347 (xpt, 1.0, 1.0),
3348 (1.0, 0.21176470816135406,0.21176470816135406)],
3349 'blue': [(0.0, 0.14901961386203766, 0.14901961386203766),
3350 (xpt, 0.7490196228027344, 0.7490196228027344),
3351 (1.0, 0.5843137502670288, 0.5843137502670288)]}
3352
3353 return matplotlib.colors.LinearSegmentedColormap('modRdYlBu', colordict)
3354
3356 """
3357 Returns a RdOrYlBu (red-orange-yellow-blue) colormap.
3358
3359 Inputs:
3360 orpt - The point between 0 and 1 where the colormap returns orange.
3361 ylpt - The point between 0 and 1 where the colormap returns yellow.
3362 Must be greater than orpt.
3363
3364 Returns:
3365 RdOrYlBu - The modified colormap.
3366 """
3367
3368 orpt = float(orpt)
3369 ylpt = float(ylpt)
3370 if (orpt < 0.0) or (orpt > 1.0):
3371 raise ValueError('RdOrYlBu: invalid orange point passed')
3372 if (ylpt < 0.0) or (ylpt > 1.0):
3373 raise ValueError('RdOrYlBu: invalid orange point passed')
3374 if (orpt >= ylpt):
3375 raise ValueError('RdOrYlBu: orange point must be less than yellow point')
3376
3377 colordict = {'red': [(0.0, 0.6470588445663452, 0.6470588445663452),
3378 (orpt, 0.9960784316062927, 0.9960784316062927),
3379 (ylpt, 1.0, 1.0),
3380 (1.0, 0.1921568661928177, 0.1921568661928177)],
3381 'green': [(0.0, 0.0, 0.0),
3382 (orpt, 0.6000000238418579, 0.6000000238418579),
3383 (ylpt, 1.0, 1.0),
3384 (1.0, 0.21176470816135406,0.21176470816135406)],
3385 'blue': [(0.0, 0.14901961386203766, 0.14901961386203766),
3386 (orpt, 0.16078431904315948, 0.16078431904315948),
3387 (ylpt, 0.7490196228027344, 0.7490196228027344),
3388 (1.0, 0.5843137502670288, 0.5843137502670288)]}
3389
3390 return matplotlib.colors.LinearSegmentedColormap('modRdYlBu', colordict)
3391
3392
3393
3394
3395
3396
3397
3398
3399
3400
3401
3402
3403
3404
3405
3406
3407
3408
3409
3410
3411
3412
3413
3414
3415
3416
3417
3418
3419 -def cmap3Color(color1, color2, color3, colorpt=0.5):
3420 """
3421 Returns a three-color colormap.
3422
3423 Inputs:
3424 color1, color2, color3 - The three colors to be used in the colormap.
3425 Should be passed as 3-tuples.
3426 colorpt - The point on the x-scale where color2 is used.
3427 Must be between 0 and 1. Defaults to 0.5.
3428
3429 Returns:
3430 colormap - The requested three-color colormap.
3431 """
3432
3433 colorpt = float(colorpt)
3434 if (colorpt < 0.0) or (colorpt > 1.0):
3435 return ValueError('cmap3Color: invalid middle color point passed')
3436
3437 colordict = {'red': [(0.0, color1[0], color1[0]),
3438 (colorpt, color2[0], color2[0]),
3439 (1.0, color3[0], color3[0])],
3440 'green': [(0.0, color1[1], color1[1]),
3441 (colorpt, color2[1], color2[1]),
3442 (1.0, color3[1], color3[1])],
3443 'blue': [(0.0, color1[2], color1[2]),
3444 (colorpt, color2[2], color2[2]),
3445 (1.0, color3[2], color3[2])]}
3446
3447 return matplotlib.colors.LinearSegmentedColormap('modRdYlBu', colordict)
3448
3450 """
3451 Define a colormap with an arbitrary number of colors.
3452
3453 Inputs:
3454 colors - A list of RGB 3-tuples to be used a colors.
3455 turnpts - A list of turning points where colors are to be used on the
3456 x (0 to 1) scale. The length of this list must be two less
3457 than the length of the colors list (as the first and last
3458 colors are automatically placed at 0.0 and 1.0). If no
3459 turning points are passed, they are automatically generated
3460 in an evenly spaced fashion. Defaults to None.
3461
3462 Returns:
3463 colormap - The requested colormap.
3464 """
3465
3466 if (turnpts != None) and (len(colors) - 2 != len(turnpts)):
3467 raise ValueError('cmapNColor: you have not specified the correct number ('+str(len(colors)-2)+') of turning points')
3468
3469 if turnpts == None:
3470 turnpts = np.arange(0.0, 1.05, 1.0/len(colors)).tolist()
3471 else:
3472 turnpts.insert(0, 0.0)
3473 turnpts.append(1.0)
3474
3475 colordict = {}
3476 colordict['red'] = []
3477 colordict['green'] = []
3478 colordict['blue'] = []
3479
3480 for i in range(len(turnpts)):
3481 colordict['red'].append((turnpts[i], colors[i][0], colors[i][0]))
3482 colordict['green'].append((turnpts[i], colors[i][1], colors[i][1]))
3483 colordict['blue'].append((turnpts[i], colors[i][2], colors[i][2]))
3484
3485 return matplotlib.colors.LinearSegmentedColormap('pyNifsNColor', colordict)
3486
3488 """
3489 Returns a three-color yellow-orange-red colormap.
3490
3491 Inputs:
3492 colorpt - The point on the color axis (0 to 1) where orange is used.
3493 Defaults to 0.5.
3494
3495 Returns:
3496 colormap - The colormap.
3497 """
3498
3499
3500 red = (0.5, 0.1, 0.15)
3501 orange = (0.996, 0.6, 0.161)
3502 yellow = (1.0, 1.0, 0.75)
3503
3504 colormap = cmap3Color(red, orange, yellow, colorpt)
3505 return colormap
3506
3508 """
3509 Returns a rainbow color map desgined by Peter McGregor.
3510
3511 Inputs:
3512 Nil.
3513
3514 Returns:
3515 RainbowPJM - A rainbow-style colormap, with only bright colors, and
3516 black and white removed.
3517 """
3518
3519 colordict = {'blue': ((0.0,0.3,0.3),
3520 (0.125,0.6,0.6),
3521 (0.25,1.0,1.0),
3522 (0.375,0.6,0.6),
3523 (0.5,0.3,0.3),
3524 (0.625,0.0,0.0),
3525 (1.0,0.0,0.0)),
3526 'green': ((0.0,0.0,0.0),
3527 (0.125,0.0,0.0),
3528 (0.25,0.6,0.6),
3529 (0.375,1.0,1.0),
3530 (0.625,1.0,1.0),
3531 (0.75, 0.6,0.6),
3532 (0.875,0.0,0.0),
3533 (1.0,0.0,0.0)),
3534 'red': ((0.0,0.0,0.0),
3535 (0.25,0.0,0.0),
3536 (0.375,0.3,0.3),
3537 (0.5,0.6,0.6),
3538 (0.625,1.0,1.0),
3539 (1.0,1.0,1.0))}
3540
3541 return matplotlib.colors.LinearSegmentedColormap('RainbowPJM',
3542 colordict)
3543
3544 -def xPV(nifscube, ypos, restwavl, vstart, vstop,
3545 x0=0.0, xstart=None, xstop=None,
3546 useoffsets=True):
3547 """
3548 Returns data for an PV diagram in the x-direction of a data cube.
3549
3550 Inputs:
3551 nifscube - The data cube to be operated on
3552 y - The y-position along which to make the cut
3553 restwavl - Rest wavelength used to determine line velocities. Must be
3554 specified in units of m.
3555 vstart - The start velocity of the PV diagram
3556 vstop - The end velocity of the PV diagram
3557 x0 - Optional. The x-position to treat as the origin of the
3558 position axis. Defaults to 0.
3559 xstart - Optional. The x-position to treat as the start of the position
3560 axis, or None to use to start of the coordinates in the
3561 nifscube. Defaults to None.
3562 xstop - Optional. The x-position to treat as the end of the position
3563 axis, or None to use to end of the coordinates in the
3564 nifscube. Defaults to None.
3565 useoffsets - Optional. Boolean value specifying whether to use offset
3566 spatial coordinates or original coordinates. Defaults to
3567 True (offset coordinates).
3568
3569 Returns:
3570 pos - The array of positions for plotting
3571 vel - The array of velocities for plotting
3572 PVdata - The data the be plotted onto the PV diagram. Conforms to
3573 the Python standard for having point (x,y) located in the
3574 array at PVdata[y,x].
3575 """
3576
3577
3578 kstart = findZIndex(nifscube, velToWavl(vstart, restwavl))
3579 kstop = findZIndex(nifscube, velToWavl(vstop, restwavl))
3580
3581 j = findYIndex(nifscube, ypos, useoffsets=useoffsets)
3582 if xstart == None:
3583 istart = 0
3584 else:
3585 istart = findXIndex(nifscube, xstart, useoffsets=useoffsets)
3586 if xstop == None:
3587 istop = len(nifscube.xpts) - 1
3588 else:
3589 istop = findXIndex(nifscube, xstop, useoffsets=useoffsets)
3590
3591
3592 wavls = nifscube.zpts[kstart:kstop+1]
3593 if nifscube.zu == "A":
3594 wavls = wavls * 1e-10
3595 elif nifscube.zu == "um":
3596 wavls = wavls * 1e-6
3597 print wavls
3598 vel = []
3599 for w in wavls:
3600 vel.append(wavlToVel(w,restwavl))
3601 vel = np.asarray(vel)
3602 if useoffsets:
3603 pos = nifscube.xptsoff[istart:istop+1]
3604 else:
3605 pos = nifscube.xpts[istart:istop+1]
3606
3607 pos = pos - x0
3608 PVdata = nifscube.scidata[kstart:kstop+1,
3609 j,
3610 istart:istop+1]
3611
3612 return pos, vel, PVdata
3613
3614 -def yPV(nifscube, xpos, restwavl, vstart, vstop,
3615 y0=0.0, ystart=None, ystop=None,
3616 useoffsets=True):
3617 """
3618 Returns data for an PV diagram in the x-direction of a data cube.
3619
3620 Inputs:
3621 nifscube - The data cube to be operated on
3622 xpos - The y-position along which to make the cut
3623 restwavl - Rest wavelength used to determine line velocities. Must be
3624 specified in units of m.
3625 vstart - The start velocity of the PV diagram
3626 vstop - The end velocity of the PV diagram
3627 y0 - Optional. The x-position to treat as the origin of the
3628 position axis. Defaults to 0.
3629 ystart - Optional. The x-position to treat as the start of the position
3630 axis, or None to use to start of the coordinates in the
3631 nifscube. Defaults to None.
3632 ystop - Optional. The x-position to treat as the end of the position
3633 axis, or None to use to end of the coordinates in the
3634 nifscube. Defaults to None.
3635 useoffsets - Optional. Boolean value specifying whether to use offset
3636 spatial coordinates or original coordinates. Defaults to
3637 True (offset coordinates).
3638
3639 Returns:
3640 pos - The array of positions for plotting
3641 vel - The array of velocities for plotting
3642 PVdata - The data the be plotted onto the PV diagram. Conforms to
3643 the Python standard for having point (x,y) located in the
3644 array at PVdata[y,x].
3645 """
3646
3647
3648 kstart = findZIndex(nifscube, velToWavl(vstart, restwavl))
3649 kstop = findZIndex(nifscube, velToWavl(vstop, restwavl))
3650
3651 i = findXIndex(nifscube, xpos, useoffsets=useoffsets)
3652 if ystart == None:
3653 jstart = 0
3654 else:
3655 jstart = findYIndex(nifscube, ystart, useoffsets=useoffsets)
3656 if ystop == None:
3657 jstop = len(nifscube.ypts) - 1
3658 else:
3659 jstop = findYIndex(nifscube, ystop, useoffsets=useoffsets)
3660
3661
3662 wavls = nifscube.zpts[kstart:kstop+1]
3663 if nifscube.zu == "A":
3664 wavls = wavls * 1e-10
3665 elif nifscube.zu == "um":
3666 wavls = wavls * 1e-6
3667 print wavls
3668 vel = []
3669 for w in wavls:
3670 vel.append(wavlToVel(w,restwavl))
3671 vel = np.asarray(vel)
3672 if useoffsets:
3673 pos = nifscube.yptsoff[jstart:jstop+1]
3674 else:
3675 pos = nifscube.ypts[jstart,jstop+1]
3676
3677 pos = pos - y0
3678 PVdata = nifscube.scidata[kstart:kstop+1,
3679 jstart:jstop+1,
3680 i]
3681
3682 return pos, vel, PVdata
3683
3684 -def xPVfitted(velfit, comp, ypos, vstart, vstop,
3685 x0=0.0, xstart=None, xstop=None):
3686 """
3687 Returns data for an PV diagram in the x-direction of a data cube.
3688
3689 Inputs:
3690 velfit - The velocity fit to be operated on
3691 comp - The component of the velocity fit to be displayed
3692 y - The y-position along which to make the cut
3693 vstart - The start velocity of the PV diagram
3694 vstop - The end velocity of the PV diagram
3695 x0 - Optional. The x-position to treat as the origin of the
3696 position axis. Defaults to 0.
3697 xstart - Optional. The x-position to treat as the start of the position
3698 axis, or None to use to start of the coordinates in the
3699 velfit. Defaults to None.
3700 xstop - Optional. The x-position to treat as the end of the position
3701 axis, or None to use to end of the coordinates in the
3702 velfit. Defaults to None.
3703
3704 Returns:
3705 pos - The array of positions for plotting
3706 vel - The array of velocities for plotting
3707 PVdata - The data the be plotted onto the PV diagram. Conforms to
3708 the Python standard for having point (x,y) located in the
3709 array at PVdata[y,x].
3710 """
3711
3712
3713
3714 comp = int(comp)
3715 if comp >= velfit.fitno.max():
3716 raise ValueError('This velfit instance does not contain that many components!')
3717 try:
3718 kstart = velfit.findVIndex(vstart)
3719 except ValueError:
3720 kstart = 0
3721 try:
3722 kstop = velfit.findVIndex(vstop)
3723 except ValueError:
3724 kstop = len(velfit.vpts) - 1
3725 j = findYIndex(velfit, ypos)
3726 if xstart == None:
3727 istart = 0
3728 else:
3729 istart = findXIndex(velfit, xstart)
3730 if xstop == None:
3731 istop = len(velfit.xpts) - 1
3732 else:
3733 istop = findXIndex(velfit, xstop)
3734
3735 pos = velfit.xpts[istart:istop+1]
3736 vel = velfit.vpts[kstart:kstop+1]
3737 PVdata = np.ma.empty([vel.shape[0], pos.shape[0]])
3738
3739 pos = pos - x0
3740
3741
3742 for k in range(len(vel)):
3743 for i in range(len(pos)):
3744 if velfit.fitno[j,istart+i] != velfit.fitno.max():
3745 PVdata[k,i] = np.ma.masked
3746 else:
3747 parmin = 3*comp + 1
3748 parmax = 3*comp + 4
3749 fitpars = ([velfit.fits[j,istart+i][velfit.fitno[j,istart+i]-1][0]] +
3750 velfit.fits[j,istart+i][velfit.fitno[j,istart+i]-1][parmin:parmax].tolist())
3751 PVdata[k,i] = fitter.n_gaussian(pars=fitpars)(vel[k])
3752
3753 return pos, vel, PVdata
3754
3755 -def yPVfitted(velfit, comp, xpos, vstart, vstop,
3756 y0=0.0, ystart=None, ystop=None):
3757 """
3758 Returns data for an PV diagram in the y-direction of a data cube.
3759
3760 Inputs:
3761 velfit - The velocity fit to be operated on
3762 comp - The component of the velocity fit to be displayed
3763 x - The x-position along which to make the cut
3764 vstart - The start velocity of the PV diagram
3765 vstop - The end velocity of the PV diagram
3766 x0 - Optional. The y-position to treat as the origin of the
3767 position axis. Defaults to 0.
3768 ystart - Optional. The y-position to treat as the start of the position
3769 axis, or None to use to start of the coordinates in the
3770 velfit. Defaults to None.
3771 ystop - Optional. The y-position to treat as the end of the position
3772 axis, or None to use to end of the coordinates in the
3773 velfit. Defaults to None.
3774
3775 Returns:
3776 pos - The array of positions for plotting
3777 vel - The array of velocities for plotting
3778 PVdata - The data the be plotted onto the PV diagram. Conforms to
3779 the Python standard for having point (x,y) located in the
3780 array at PVdata[y,x].
3781 """
3782
3783
3784
3785 comp = int(comp)
3786 if comp >= velfit.fitno.max():
3787 raise ValueError('This velfit instance does not contain that many components!')
3788 try:
3789 kstart = velfit.findVIndex(vstart)
3790 except ValueError:
3791 kstart = 0
3792 try:
3793 kstop = velfit.findVIndex(vstop)
3794 except ValueError:
3795 kstop = len(velfit.vpts) - 1
3796 i = findXIndex(velfit, xpos)
3797 if ystart == None:
3798 jstart = 0
3799 else:
3800 jstart = findYIndex(velfit, ystart)
3801 if ystop == None:
3802 jstop = len(velfit.ypts) - 1
3803 else:
3804 jstop = findYIndex(velfit, ystop)
3805
3806 pos = velfit.ypts[jstart:jstop+1]
3807 vel = velfit.vpts[kstart:kstop+1]
3808 PVdata = np.ma.empty([vel.shape[0], pos.shape[0]])
3809
3810 pos = pos - y0
3811
3812
3813 for k in range(len(vel)):
3814 for j in range(len(pos)):
3815 if velfit.fitno[jstart+j,i] != velfit.fitno.max():
3816 PVdata[k,j] = np.ma.masked
3817 else:
3818 parmin = 3*comp + 1
3819 parmax = 3*comp + 4
3820 fitpars = ([velfit.fits[jstart+j,i][velfit.fitno[jstart+j,i]-1][0]] +
3821 velfit.fits[jstart+j,i][velfit.fitno[jstart+j,i]-1][parmin:parmax].tolist())
3822 PVdata[k,j] = fitter.n_gaussian(pars=fitpars)(vel[k])
3823
3824 return pos, vel, PVdata
3825