trafo_Descriptive.py
4.7 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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
title = "Descriptive <-> Absolute"
tip = "regrids a map of projection and coordinate system into a different one"
onein = True
import numpy as np
from kapteyn import wcs
from kapteyn.interpolation import map_coordinates
from nodmath import map_interpolate
from nodmath import extract
from nodastro import nodtrafo
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 _
def nint(x):
if x > 0: return int(x+0.5)
else: return int(x-0.5)
nt = nodtrafo()
class NOD3_App:
def __init__(self, parent):
self.parent = parent
self.parent.activateWindow()
def compute_app(self, **args):
class FuncParam(DataSet):
system = ChoiceItem("Zero point", (("absolute", "absoulte"),
("descriptive", "descriptive")))
name = title.replace(" ", "")
if args == {}:
param = FuncParam(_(title), "Your choise:")
else:
param = self.parent.ScriptParameter(name, args)
#param = None
self.parent.compute_11(name, lambda m, p: self.function(m, p), param, onein)
def Error(self, msg):
QMessageBox.critical(self.parent.parent(), title,
_(u"Error:")+"\n%s" % str(msg))
def automatic(self, header, mat):
L, B = [0,0,0,0], [0,0,0,0]
X, Y = (0, header['NAXIS1'], 0, header['NAXIS1']), (0, 0, header['NAXIS2'], header['NAXIS2'])
for i in range(4):
x = header['CRVAL1'] + (X[i] - header['CRPIX1'])*header['CDELT1']
y = header['CRVAL2'] + (Y[i] - header['CRPIX2'])*header['CDELT2']
L[i], B[i] = nt.getCoord(x, y, mat)
c = np.cos(np.pi/180.0 * (min(B) + (max(B)-min(B))*0.5))
ndim = int(c*abs((max(L)-min(L))/header['CDELT1'])+1.5)
ldim = int(abs((max(L)-min(L))/header['CDELT1'])+1.5)
mdim = int(abs((max(B)-min(B))/header['CDELT2'])+1.5)
if header['CDELT1'] < 0:
l = max(L) + 0.5*ldim*header['CDELT1']
else:
l = min(L) + 0.5*ndim*header['CDELT1']
l = (l+360.0) % 360.0
b = min(B) + 0.5*mdim*header['CDELT2']
return (mdim, ndim), 0.5*ndim, 0.5*mdim, l, b
def function(self, m, p):
naxis = m.header['NAXIS']
header = m.header.copy()
header['CTYPE1'] = header['CTYPE1'].replace("DES", "CAR")
header['CTYPE2'] = header['CTYPE2'].replace("DES", "CAR")
if header['CTYPE1'].find("CAR") < 0:
self.Error(_("this app works on CAR projection only."))
return [], p
header['NAXIS'] = 2
proj1 = wcs.Projection(header)
newheader = header.copy()
if 'BLON' in header and 'BLAT' in header:
blon, blat = header['BLON'], header['BLAT']
else:
blon = header['CRVAL1']
blat = header['CRVAL2']
newheader['BLON'] = blon
newheader['BLAT'] = blat
if p.system.lower() == "descriptive":
newheader['CRVAL1'], newheader['CRVAL2'] = blon, blat
newheader['CRPIX1'], newheader['CRPIX2'] = proj1.topixel((blon, blat))
if abs(newheader['CRPIX1']*header['CDELT1']) > 180.0:
newheader['CRPIX1'] += 360./header['CDELT1']
else:
sfl = np.cos(header['CRVAL2']*np.pi/180.0)
newheader['NAXIS1'] = int(header['NAXIS1'] / sfl)
pix, piy = (newheader['NAXIS1']+1)/2.0, (newheader['NAXIS2']+1)/2.0
newheader['CRVAL1'], newheader['CRVAL2'] = 0.0, 0.0
newheader['CRPIX1'] = header['CRPIX1']/sfl - header['CRVAL1']/header['CDELT1']
newheader['CRPIX2'] -= header['CRVAL2']/header['CDELT2']
#newheader['CRPIX1'], newheader['CRPIX2'] = proj1.topixel((0.0, 0.0))
proj2 = wcs.Projection(newheader)
x, y = wcs.coordmap(proj1, proj2)
if abs(np.median(y)*header['CDELT1']) > 180.0:
y -= 360./header['CDELT2']
#print newheader['CRPIX1'], header['CRPIX1'],
#print newheader['CRPIX2'], header['CRPIX2']
#print "================================================="
#print x
#print "-------------------------------------------------"
#print y
#print "================================================="
if p.system.lower() == "descriptive":
newheader['CTYPE1'] = newheader['CTYPE1'].replace("CAR", "DES")
newheader['CTYPE2'] = newheader['CTYPE2'].replace("CAR", "DES")
newdata = map_interpolate(m.data, y, x)
m.header = newheader
m.data = newdata
m = self.parent.removeNANedges(m)
m.header['NAXIS'] = naxis
return m, p