base_refadjust.py
3.15 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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
title = "RefAdjust"
tip = "tool tip"
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 nodmath import nan_interpolation, subpixel_shift, same_size
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):
#s = StringItem('s', default="string")
#i = IntItem('i', default=0, max=100, min=0)
#a = FloatItem('a', default=1.)
#b = BoolItem("bool", default=True)
#c = ChoiceItem(_(u"Mode"), zip(units, units), default=units[0])
clip = FloatItem('Cutoff', default=0.1, min=0.0)
name = title.replace(" ", "")
if args == {}:
param = FuncParam(_(title), "RefAdjust parameter (Refmap first):")
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 ref_adjust(self, rdata, data, clip=None, order=1):
diff = rdata - data
if clip != None:
diff = self.sig_clip(diff, clip, 3)
rows, cols = diff.shape
x = np.arange(cols)
y = np.arange(rows)
for row in range(rows):
mask = ~np.isnan(diff[row])
if len(x[mask]) > order:
poly = np.poly1d(np.polyfit(x[mask], diff[row][mask], order))
data[row] += poly(x)
diff = rdata - data
if clip != None:
diff = self.sig_clip(diff, clip, 3)
rows, cols = diff.shape
for col in range(cols):
mask = ~np.isnan(diff[:,col])
if len(y[mask]) > order:
poly = np.poly1d(np.polyfit(y[mask], diff[:,col][mask], order))
data[:,col] += poly(y)
return data
def sig_clip(self, a, clip, n):
if clip > 0:
for i in range(n):
c = clip*np.std(a)
a = np.clip(a, -c, c)
return a
def function(self, ms, p):
out, ms = same_size([ms[1], ms[0]], True)
#return ms, p
if out > 0:
self.Error("one or more maps have different scales or coordinate systems")
return [], p
m = ms[0]
mr = ms[1]
mask = ~np.isnan(m.data)
m.data = nan_interpolation(m.data)
mr.data = nan_interpolation(mr.data)
dx, dy, mr.data = subpixel_shift(m.data, mr.data, clip=p.clip)
a, b = correlate(mr.data, m.data, fmax=m.data.min(), out=False, sort=False)
mr.data = a + b*mr.data
m.data = self.ref_adjust(mr.data, m.data, clip=3.0)
m.data = np.where(mask, m.data, np.nan)
return m, p