source: mire/process.py

Last change on this file was 34, checked in by meynadie, 15 years ago

Ajout 00README

  • Property svn:executable set to *
File size: 13.3 KB
Line 
1#!/bin/env python
2# coding: utf-8
3# Un programme pour exploiter les données des mires
4# -- Frédéric Meynadier
5
6import string
7from numpy import *
8
9from functions import *
10
11from datetime import *
12
13
14from optparse import OptionParser
15parser = OptionParser()
16parser.add_option("-c","--catalog",dest="catalogfile",
17                  help="Name of the catalog file")
18parser.add_option("-r","--reference",dest="reference",
19                  help="Name of the file with reference coordinates")
20parser.add_option("-o","--outputfits",dest="outputfits",
21                  default="output.fits",
22                  help="Name of the output fits file")
23parser.add_option("-m","--measurement",dest="measurement",
24                  default="output.mes",
25                  help="Name of the output measurement file")
26parser.add_option("--ring",dest="ring",
27                  default="output.ring",
28                  help="Name of the output measurement file")
29parser.add_option("-s","--square",dest="square",
30                  default="output.square",
31                  help="Name of the output measurement file")
32parser.add_option("-t","--timing",dest="timing",
33                  help="Name of the input timing file")
34
35
36(options, args) = parser.parse_args ()
37
38# Recherche de la date dans un fichier des timings
39# (généré par recup_date, soit :
40# ls -l $1| awk '{print$8 "\t" $7 "\t" $6}' | sort -k2 | sed -e 's/*//'
41# où $1 est le répertoire contenant les .dat non modifiés.)
42# format attendu : nom_fichier phase_cycle date_isoformat
43
44for line in file(options.timing):
45    line_s = string.split(line)
46    if (len(line_s) > 2):
47        nom = line_s[0] [:-4]
48        if (nom == options.catalogfile[:-6]):
49            dt = datetime.strptime(line_s[2],
50                                   "%Y-%m-%dT%H:%M:%S")
51            phase = float(line_s[1])
52
53# Donner les coordonnées approximatives des 4 repÚres
54
55rep = zeros((4,2))
56i = 0
57for line in open(options.reference):
58    if ((line[0] != "#") and (line != '\n')):
59        spl = string.split(line)
60        rep[i][0] = float(spl[0])
61        rep[i][1] = float(spl[1])
62        i+=1
63
64# Lecture du fichier catalogue
65fields = ["FLUX_MAX",
66          "X_IMAGE",
67          "Y_IMAGE",
68          "THETA_IMAGE",
69          "ELONGATION",
70          "ELLIPTICITY",
71          "FWHM_IMAGE",
72          "ERRX2_IMAGE",
73          "ERRY2_IMAGE"]
74
75catalog = []
76for line in open(options.catalogfile):
77    if ((line[0] != "#") and (line != '\n')):
78        tmp = {}
79        spl = string.split(line)
80        for i in range(len(fields)):
81            tmp[fields[i]] = float(spl[i])
82        tmp['MEAN_DIST_TO_NGBRS']=0
83        catalog.append(tmp)
84
85
86mode_degrade = 0
87# Recherche de la position des repÚres selon catalogue - attention aux doublons
88real_rep = zeros((4,2))
89for i in range(4):
90    detections = 0
91    for source in catalog:
92        pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])
93        max_dist = 5 # pixels
94        if ((modulus(rep[i,:] - pos_source)) < max_dist):
95            max_dist = modulus(rep[i,:] - pos_source)
96            detections +=1
97            real_rep[i,:] = pos_source
98    if (detections > 1):
99        print "Attention ! ambiguïté pour repÚre ",i,\
100              ", nb de détections :",detections
101       
102    if (detections ==0):
103        print "Erreur ! Pas de détection pour repÚre ",i
104        real_rep[i,:]=rep[i,:]
105        mode_degrade = 1
106
107if (mode_degrade==1):
108    print "Passage en mode dégradé, sorties réduites"
109    exit(2)
110
111# Détermination du centre : barycentre des 4 repÚres
112center = array([sum(real_rep[:,0])/4,sum(real_rep[:,1])/4])
113
114# Vecteurs unitaires de référence
115vec_i = ((real_rep[0,:] - real_rep[1,:])/20 \
116         + (real_rep[3,:] - real_rep[2,:])/20) / 2 
117vec_j = ((real_rep[1,:] - real_rep[2,:])/20 \
118         + (real_rep[0,:] - real_rep[3,:])/20) / 2 
119
120# Détection des noeuds de grille les plus proches de chaque source
121vec_i_2 = dot(vec_i,vec_i)
122vec_j_2 = dot(vec_j,vec_j)
123for source in catalog:
124    pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])-center
125    # dans le référentiel grille (en projetant sur les vecteurs unitaires) :
126    pos_source_grid = array([dot(pos_source,vec_i)/vec_i_2,
127                             dot(pos_source,vec_j)/vec_j_2])
128    # on arrondit au noeud le plus proche
129    knot = pos_source_grid.round()
130    # et on calcule la distance entre les deux
131    r = modulus(pos_source_grid - knot)
132    # enregistrement des résultats
133    source['I_GRID'] = knot[0]
134    source['J_GRID'] = knot[1]
135    source['KNOT_DIST'] = r
136    #if (modulus(knot) <1):
137    #    print source['I_GRID'], source['J_GRID'], source['KNOT_DIST'], pos_source_grid, pos_source,source['X_IMAGE'], source['Y_IMAGE']
138
139# Chasse aux doublons
140# tableau 81*81 contenant :
141# champ 0 : nb de sources rattachée au noeud
142# champ 1 : distance actuelle entre le noeud et la source
143# champ 2 : index courant de la source dans catalog
144# champ 3 : distance moyenne aux voisins (calculé plus tard)
145lookup_sources = zeros((81,81,4))
146for si in range(len(catalog)):
147    source = catalog[si]
148    i = source['I_GRID']+40
149    j = source['J_GRID']+40
150    if ((lookup_sources[i,j,0] < 1)
151        or (lookup_sources[i,j,1] > source['KNOT_DIST'])):
152        lookup_sources[i,j,2] = si
153        lookup_sources[i,j,1] = source['KNOT_DIST']
154    lookup_sources[i,j,0] += 1
155
156# On élimine le noeud (1,1), trop prÚs du "P"
157lookup_sources[41,41,0] = 0
158
159# Calcul de la distance moyenne entre noeuds voisins
160
161fen = 2
162
163shift_list = []
164
165for i in range(-fen,fen+1):
166    for j in range(-fen,fen+1):
167        if (i!=0 or j!=0):
168            shift_list.append([i,j])
169
170for i in range(fen,81-fen):
171    for j in range(fen,81-fen):
172        if (lookup_sources[i,j,0] > 0):
173            knot_idx = int(lookup_sources[i,j,2])
174            catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = 0
175            nb_voisin = 0
176            distance = 0
177            weights = 0
178            pos_source = array([catalog[knot_idx]['X_IMAGE'],
179                                catalog[knot_idx]['Y_IMAGE']])
180            for shifts in shift_list :
181                i_voisin = i + shifts[0]
182                j_voisin = j + shifts[1]
183                if (lookup_sources[i_voisin,j_voisin,0] > 0):
184                    neigh_idx = int(lookup_sources[i_voisin,j_voisin,2])
185                    nb_voisin +=1
186                    pos_voisin = array([catalog[neigh_idx]['X_IMAGE'],
187                                        catalog[neigh_idx]['Y_IMAGE']])
188                    distance += modulus(pos_source-pos_voisin)
189                    weights += modulus(shifts)
190            if (nb_voisin > 0):
191                catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = float(distance)/weights
192                lookup_sources[i,j,3] = (float(distance)/weights)
193       
194cleaned_catalog = []
195for i in range(81):
196    for j in range(81):
197        if (lookup_sources[i,j,0] > 0):
198            cleaned_catalog.append(catalog[int(lookup_sources[i,j,2])])
199
200# Sortie du catalogue nettoyé des doublons
201
202of = open('clean.cat','w')
203
204fields = ["FLUX_MAX",
205          "X_IMAGE",
206          "Y_IMAGE",
207          "THETA_IMAGE",
208          "ELONGATION",
209          "ELLIPTICITY",
210          "FWHM_IMAGE",
211          "ERRX2_IMAGE",
212          "ERRY2_IMAGE",
213          "I_GRID",
214          "J_GRID",
215          "KNOT_DIST",
216          "MEAN_DIST_TO_NGBRS"]
217
218fmts = ["%10.3f ",
219        "%8.3f ",
220        "%8.3f ",
221        "%8.3f ",
222        "%8.3f ",
223        "%8.3f ",
224        "%8.3f ",
225        "%8.5f ",
226        "%8.5f ",
227        "%5d ",
228        "%5d ",
229        "%5.3f ",
230        "%8.3f"]
231
232for source in cleaned_catalog:
233    outs = ""
234    for i in range(len(fields)):
235        outs += fmts[i] % source[fields[i]]
236    outs += "\n"
237    of.write(outs)
238
239of.close()
240
241# Génération d'un fichier fits contenant les résultats
242
243
244import pyfits
245
246#a = zeros((2048,2048))
247#
248# patches width in pixel:
249## width = int(modulus(vec_i))
250
251## for i in range(len(cleaned_catalog)):
252##     print i
253##     source = cleaned_catalog[i]
254##     pos_source = array([int(source['X_IMAGE']),int(source['Y_IMAGE'])])
255##     pos_mean_dist = source['MEAN_DIST_TO_NGBRS']
256##     for rel_x in range(-width,width):
257##         for rel_y in range(-width,width):
258##             rel_vec = array([rel_x,rel_y])
259##             # Forme du patch : carré
260##             if ( (dot(rel_vec,vec_i) < vec_i_2/2) and
261##                  (dot(rel_vec,vec_j) < vec_j_2/2)):
262##                 (x,y) = pos_source+rel_vec
263##                 if ( x>0 and x<2048 and y>0 and y<2048):
264##                     a[x,y] = pos_mean_dist
265
266
267
268a = lookup_sources[:,:,3]
269
270for i in range(81):
271    for j in range(81):
272        if (a[i,j] > 0):
273            a[i,j] = 60.0/a[i,j]
274        else :
275            a[i,j] = 1.06
276
277# Interpolation de la case du "P"
278
279a[41,41] = (a[42,41] + a[40,41] + a[41,42] + a[41,40]) /4.
280
281hdu = pyfits.PrimaryHDU(a)
282hdulist = pyfits.HDUList([hdu])
283
284hdu.writeto(options.outputfits,clobber=True)
285
286
287# Facteur d'échelle moyen
288
289tmp = []
290for i in range(81):
291    for j in range(81):
292        if (lookup_sources[i,j,0] == 1):
293            tmp.append(lookup_sources[i,j,3])
294
295
296meas = open(options.measurement,'w')
297
298meas.write("# Moyenne des mesures /  Module vec_i /  Module vec_j / deltaT (min) / date\n")
299# transfo pixels / arcminutes -> mas / pixel
300moyenne_mesures = 60 * (1. / (sum(tmp)/len(tmp)))
301module_veci = 60 * (1. /  modulus(vec_i))
302module_vecj = 60 * (1. /  modulus(vec_j))
303
304outs = "%10.7f %10.7f %10.7f %7.3f %s\n" % (moyenne_mesures,
305                                            module_veci,
306                                            module_vecj,
307                                            phase,
308                                            dt.isoformat())
309meas.write(outs)
310
311meas.close()
312
313## Facteur d'échelle par rapport au centre du champ sur une couronne
314## représentative du diamÚtre solaire
315
316# Quel noeud est (à peu prÚs) au centre ?
317
318
319dist_mini = 2048
320i_centre = 0
321for i in range(len(cleaned_catalog)):
322    source = cleaned_catalog[i]
323    dist = modulus([source["X_IMAGE"] - 1024, source["Y_IMAGE"] - 1024])
324    if (dist < dist_mini):
325        dist_mini = dist
326        i_centre = i
327
328#print "Source au centre : ",cleaned_catalog[i_centre]
329
330# On parcourt une couronne de noeuds et on calcule le facteur d'échelle pour chacun
331
332rmin = 14 # en inter-noeud
333rmax = 18 # en inter-noeud
334
335source_centrale = cleaned_catalog[i_centre]
336
337couronne=[]
338
339for source in cleaned_catalog:
340    tmp = {}
341    dist_knots = modulus([source["I_GRID"] - source_centrale["I_GRID"],
342                          source["J_GRID"] - source_centrale["J_GRID"]])
343    x_rel = source["X_IMAGE"] - source_centrale["X_IMAGE"]
344    y_rel = source["Y_IMAGE"] - source_centrale["Y_IMAGE"]
345    dist_pix = modulus([x_rel,y_rel])
346    if ((dist_knots >= rmin) and (dist_knots <= rmax)):
347        tmp['ScFact'] = dist_knots * 60. / dist_pix # arcsec / pix
348        tmp['r'] = dist_pix
349        tmp['theta'] = math.atan(y_rel/x_rel)
350        if (x_rel < 0):
351            tmp['theta']+=math.pi
352        # Offset de pi pour se mettre en convention "astro"
353        # theta = 0 -> P droit quand nord vers le haut
354        tmp['theta']+=math.pi
355        while (tmp['theta'] < 0):
356            tmp['theta'] += 2*math.pi
357        while (tmp['theta'] >= 2*math.pi):
358            tmp['theta'] -= 2*math.pi
359        # On n'écrit pas les valeurs aberrantes
360        if ((tmp['ScFact'] > 1.05) and (tmp['ScFact'] < 1.07) or 1):
361            couronne.append(tmp)
362
363output_couronne = file(options.ring,'w')
364
365outs = "# phase: %f date: %s\n" % (phase,
366                                 dt.isoformat())
367output_couronne.write(outs)
368
369for source in couronne:
370    outs = "%f %f %10.7f\n" % (source['r'],
371                               math.degrees(source['theta']),
372                               source['ScFact'])
373    output_couronne.write(outs)
374
375
376
377output_couronne.close()
378       
379# On parcourt un carré de surface équivalente à la couronne définie ci dessus
380# Le facteur d'échelle est calculé par rapport au point opposé du carré
381# soit un cÎté de 28 noeuds
382
383
384# Le centre a été trouvé pour la couronne
385i_cen = cleaned_catalog[i_centre]["I_GRID"]
386j_cen = cleaned_catalog[i_centre]["J_GRID"]
387
388# On se refait un "lookup" simplifié
389
390synthese = zeros((81,81,2))
391
392for source in cleaned_catalog:
393    synthese[source['I_GRID'],source['J_GRID'],0] = source['X_IMAGE']
394    synthese[source['I_GRID'],source['J_GRID'],1] = source['Y_IMAGE']
395
396valeurs =  [] 
397for i in range(-14,15):
398    vecteur = synthese[i_cen+14,j_cen+i,:] - synthese[i_cen-14,j_cen+i,:]
399    valeurs.append(modulus(vecteur))
400
401valeurs_orth = []
402for i in range(-14,15):
403    vecteur = synthese[i_cen+i,j_cen+14,:] - synthese[i_cen+i,j_cen-14,:]
404    valeurs_orth.append(modulus(vecteur))
405
406for i in range(len(valeurs)):
407    if valeurs[i] > 0:
408        valeurs[i] = 28*60./valeurs[i]
409
410for i in range(len(valeurs_orth)):
411    if valeurs_orth[i] > 0:
412        valeurs_orth[i] = 28*60./valeurs_orth[i]
413
414
415# on vire les valeurs aberrantes
416def f(x):
417    return ((x < 1.07) and (x > 1.05))
418
419valeurs = filter(f, valeurs)
420valeurs_orth = filter(f, valeurs_orth)
421
422output_square = open(options.square,'w')
423outs = "# moyennes: %f / %f phase: %f date: %s\n" % (mean(valeurs),
424                                                     mean(valeurs_orth),
425                                                     phase,
426                                                     dt.isoformat())
427output_square.write(outs)
428
429for valeur in valeurs:
430    outs = "%f\n" % (valeur)
431    output_square.write(outs)
432
433output_square.write('orth \n')
434
435for valeur in valeurs_orth:
436    outs = "%f\n" % (valeur)
437    output_square.write(outs)
438
439output_square.close()
440
Note: See TracBrowser for help on using the repository browser.