proc_spectindex.py
2.07 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
title = "Spectral Index"
tip = "spectral index"
onein = 2
import numpy as np
from guidata.qt.QtGui import QMessageBox
from guidata.dataset.datatypes import DataSet
from guidata.dataset.dataitems import (IntItem, StringItem, ChoiceItem, FloatItem, BoolItem)
from guiqwt.config import _
from nodfitting import correlate
class NOD3_App:
def __init__(self, parent):
self.parent = parent
self.parent.activateWindow()
def Error(self, msg):
QMessageBox.critical(self.parent.parent(), title,
_(u"Error:")+"\n%s" % str(msg))
def compute_app(self, **args):
units = ('Pixel', 'Degree', 'Arcmin', 'Arcsec')
class FuncParam(DataSet):
rms = FloatItem('RMS', default=0.)
name = title.replace(" ", "")
if args == {}:
param = FuncParam(_(title), "cuts values at 3*RMS")
else:
param = self.parent.ScriptParameter(name, args)
# if no parameter needed set param to None. activate next line
#param = None
self.parent.compute_11(name, lambda m, p: self.function(m, p), param, onein)
def function(self, ms, p):
m1 = ms[0]
m2 = ms[1]
#a, b = correlate(m1.data, m2.data, fmax=p.rms, sort=False)
#m1.data += a
freq1 = m1.header['CRVAL3']
freq2 = m2.header['CRVAL3']
mask1 = ~np.isnan(m1.data)
mask2 = ~np.isnan(m2.data)
mask1 *= m1.data > 3*p.rms
mask2 *= m2.data > 3*p.rms
f = np.log10(freq1/freq2)
dm1 = +1.0/(np.log(10.0)*f*m1.data) # partial derivation m1.data
dm2 = -1.0/(np.log(10.0)*f*m2.data) # partial derivation m2.data
m1.data = np.where(mask1*mask2, np.log10(m1.data/m2.data) / np.log10(freq1/freq2), np.nan)
m2.data = p.rms*np.sqrt(dm1**2 + dm2**2)
m2.data = np.where(mask1*mask2, m2.data, np.nan)
del m1.header['CRVAL3']
del m2.header['CRVAL3']
m1.header['MAPTYPE'] = "SI"
m2.header['MAPTYPE'] = "Error"
m1.header['BUNIT'] = ""
m2.header['BUNIT'] = ""
return [m2, m1], p