filter_asmooth.py
2.38 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
title = "Gaussian Smooth"
tip = "applies Gaussian filter to data (smoothing)"
import numpy as np
import scipy.ndimage as spi
from guidata.dataset.datatypes import DataSet
from guidata.dataset.dataitems import (IntItem, FloatArrayItem, StringItem,
ChoiceItem, FloatItem, DictItem,
BoolItem)
from guiqwt.config import _
from nodmath import nan_interpolation
class NOD3_App:
def __init__(self, parent):
self.parent = parent
self.parent.activateWindow()
def compute_app(self):
class Param(DataSet):
hpbwx = FloatItem(u"HPBWx", default=1.)
hpbwy = FloatItem(u"HPBWy", default=1.)
Unit = ChoiceItem("Unit", (("Degree", "Degree"), ("Arcmin", "Arcmin"),
("Arcsec", "Arcsec")), default = "Arcmin")
#Scale = BoolItem("Scale", default=True)
param = Param(_("Gaussian filter"), "smooth data")
param.Scale = True
name = title.replace(" ", "")
self.parent.compute_11(name, lambda m, p: self.function(m, p), param)
def function(self, m, p):
if p.Unit == "Degree": u = 1.0
elif p.Unit == "Arcmin": u = 1.0/60.0
elif p.Unit == "Arcsec": u = 1.0/3600.0
f = 2.0*np.sqrt(2.0*np.log(2.0))
if "BMAJ" in m.header:
hpbwx = m.header["BMAJ"]
if "BMIN" in m.header:
hpbwy = m.header["BMIN"]
else:
hpbwy = m.header["BMAJ"]
else:
fghz = m.header['CRVAL3']*1.e-9
diam = m.header['DIAMETER']
clam = 0.299792458 / fghz
hpbw = 180.0/np.pi * 1.22 * clam / diam
hpbwx, hpbwy = hpbw, hpbw
delx = abs(m.header['CDELT1'])
dely = abs(m.header['CDELT2'])
sigmax = np.sqrt((u*p.hpbwx/delx)**2 - (hpbwx/delx)**2) / f
sigmay = np.sqrt((u*p.hpbwy/dely)**2 - (hpbwy/dely)**2) / f
m.header["BMAJ"] = max(u*p.hpbwx, u*p.hpbwy)
m.header["BMIN"] = min(u*p.hpbwx, u*p.hpbwy)
m.header["BPA"] = 0.0
if p.Scale:
f = u**2*p.hpbwx*p.hpbwy / (hpbwx*hpbwy)
#print f
else:
f = 1
mask = np.isnan(m.data)
m.data = nan_interpolation(m.data)
m.data = f*spi.gaussian_filter(m.data, (sigmay, sigmax))
m.data = np.where(mask, np.nan, m.data)
return m, p