source: mire/process.py @ 27

Last change on this file since 27 was 27, checked in by meynadie, 16 years ago

passage en timing_cycle

  • Property svn:executable set to *
File size: 8.4 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("-t","--timing",dest="timing",
27                  help="Name of the input timing file")
28
29
30(options, args) = parser.parse_args ()
31
32# Recherche de la date dans un fichier des timings
33# (généré par recup_date, soit :
34# ls -l $1| awk '{print$8 "\t" $7 "\t" $6}' | sort -k2 | sed -e 's/*//'
35# où $1 est le répertoire contenant les .dat non modifiés.)
36# format attendu : nom_fichier phase_cycle date_isoformat
37
38for line in file(options.timing):
39    line_s = string.split(line)
40    if (len(line_s) > 2):
41        nom = line_s[0] [:-4]
42        if (nom == options.catalogfile[:-6]):
43            dt = datetime.strptime(line_s[2],
44                                   "%Y-%m-%dT%H:%M:%S")
45            phase = float(line_s[1])
46
47# Donner les coordonnées approximatives des 4 repÚres
48
49rep = zeros((4,2))
50i = 0
51for line in open(options.reference):
52    if ((line[0] != "#") and (line != '\n')):
53        spl = string.split(line)
54        rep[i][0] = float(spl[0])
55        rep[i][1] = float(spl[1])
56        i+=1
57
58# Lecture du fichier catalogue
59fields = ["FLUX_MAX",
60          "X_IMAGE",
61          "Y_IMAGE",
62          "THETA_IMAGE",
63          "ELONGATION",
64          "ELLIPTICITY",
65          "FWHM_IMAGE",
66          "ERRX2_IMAGE",
67          "ERRY2_IMAGE"]
68
69catalog = []
70for line in open(options.catalogfile):
71    if ((line[0] != "#") and (line != '\n')):
72        tmp = {}
73        spl = string.split(line)
74        for i in range(len(fields)):
75            tmp[fields[i]] = float(spl[i])
76        tmp['MEAN_DIST_TO_NGBRS']=0
77        catalog.append(tmp)
78
79
80# Recherche de la position des repÚres selon catalogue - attention aux doublons
81real_rep = zeros((4,2))
82for i in range(4):
83    detections = 0
84    for source in catalog:
85        pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])
86        if ((modulus(rep[i,:] - pos_source)) < 5):
87            detections +=1
88            real_rep[i,:] = pos_source
89    if (detections !=1):
90        print "Attention ! ambiguïté pour repÚre ",i,\
91              ", nb de détections :",detections
92       
93
94# Détermination du centre : barycentre des 4 repÚres
95center = array([sum(real_rep[:,0])/4,sum(real_rep[:,1])/4])
96
97# Vecteurs unitaires de référence
98vec_i = ((real_rep[0,:] - real_rep[1,:])/20 \
99         + (real_rep[3,:] - real_rep[2,:])/20) / 2 
100vec_j = ((real_rep[1,:] - real_rep[2,:])/20 \
101         + (real_rep[0,:] - real_rep[3,:])/20) / 2 
102
103# Détection des noeuds de grille les plus proches de chaque source
104vec_i_2 = dot(vec_i,vec_i)
105vec_j_2 = dot(vec_j,vec_j)
106for source in catalog:
107    pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])-center
108    # dans le référentiel grille (en projetant sur les vecteurs unitaires) :
109    pos_source_grid = array([dot(pos_source,vec_i)/vec_i_2,
110                             dot(pos_source,vec_j)/vec_j_2])
111    # on arrondit au noeud le plus proche
112    knot = pos_source_grid.round()
113    # et on calcule la distance entre les deux
114    r = modulus(pos_source_grid - knot)
115    # enregistrement des résultats
116    source['I_GRID'] = knot[0]
117    source['J_GRID'] = knot[1]
118    source['KNOT_DIST'] = r
119    #if (modulus(knot) <1):
120    #    print source['I_GRID'], source['J_GRID'], source['KNOT_DIST'], pos_source_grid, pos_source,source['X_IMAGE'], source['Y_IMAGE']
121
122# Chasse aux doublons
123# tableau 81*81 contenant :
124# champ 0 : nb de sources rattachée au noeud
125# champ 1 : distance actuelle entre le noeud et la source
126# champ 2 : index courant de la source dans catalog
127# champ 3 : distance moyenne aux voisins (calculé plus tard)
128lookup_sources = zeros((81,81,4))
129for si in range(len(catalog)):
130    source = catalog[si]
131    i = source['I_GRID']+40
132    j = source['J_GRID']+40
133    if ((lookup_sources[i,j,0] < 1)
134        or (lookup_sources[i,j,1] > source['KNOT_DIST'])):
135        lookup_sources[i,j,2] = si
136        lookup_sources[i,j,1] = source['KNOT_DIST']
137    lookup_sources[i,j,0] += 1
138
139# On élimine le noeud (1,1), trop prÚs du "P"
140lookup_sources[41,41,0] = 0
141
142# Calcul de la distance moyenne entre noeuds voisins
143
144fen = 2
145
146shift_list = []
147
148for i in range(-fen,fen+1):
149    for j in range(-fen,fen+1):
150        if (i!=0 or j!=0):
151            shift_list.append([i,j])
152
153for i in range(fen,81-fen):
154    for j in range(fen,81-fen):
155        if (lookup_sources[i,j,0] > 0):
156            knot_idx = int(lookup_sources[i,j,2])
157            catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = 0
158            nb_voisin = 0
159            distance = 0
160            weights = 0
161            pos_source = array([catalog[knot_idx]['X_IMAGE'],
162                                catalog[knot_idx]['Y_IMAGE']])
163            for shifts in shift_list :
164                i_voisin = i + shifts[0]
165                j_voisin = j + shifts[1]
166                if (lookup_sources[i_voisin,j_voisin,0] > 0):
167                    neigh_idx = int(lookup_sources[i_voisin,j_voisin,2])
168                    nb_voisin +=1
169                    pos_voisin = array([catalog[neigh_idx]['X_IMAGE'],
170                                        catalog[neigh_idx]['Y_IMAGE']])
171                    distance += modulus(pos_source-pos_voisin)
172                    weights += modulus(shifts)
173            if (nb_voisin > 0):
174                catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = float(distance)/weights
175                lookup_sources[i,j,3] = float(distance)/weights
176       
177cleaned_catalog = []
178for i in range(81):
179    for j in range(81):
180        if (lookup_sources[i,j,0] > 0):
181            cleaned_catalog.append(catalog[int(lookup_sources[i,j,2])])
182
183# Sortie du catalogue nettoyé des doublons
184
185of = open('clean.cat','w')
186
187fields = ["FLUX_MAX",
188          "X_IMAGE",
189          "Y_IMAGE",
190          "THETA_IMAGE",
191          "ELONGATION",
192          "ELLIPTICITY",
193          "FWHM_IMAGE",
194          "ERRX2_IMAGE",
195          "ERRY2_IMAGE",
196          "I_GRID",
197          "J_GRID",
198          "KNOT_DIST",
199          "MEAN_DIST_TO_NGBRS"]
200
201fmts = ["%10.3f ",
202        "%8.3f ",
203        "%8.3f ",
204        "%8.3f ",
205        "%8.3f ",
206        "%8.3f ",
207        "%8.3f ",
208        "%8.5f ",
209        "%8.5f ",
210        "%5d ",
211        "%5d ",
212        "%5.3f ",
213        "%8.3f"]
214
215for source in cleaned_catalog:
216    outs = ""
217    for i in range(len(fields)):
218        outs += fmts[i] % source[fields[i]]
219    outs += "\n"
220    of.write(outs)
221
222of.close()
223
224# Génération d'un fichier fits contenant les résultats
225
226
227import pyfits
228
229#a = zeros((2048,2048))
230#
231# patches width in pixel:
232## width = int(modulus(vec_i))
233
234## for i in range(len(cleaned_catalog)):
235##     print i
236##     source = cleaned_catalog[i]
237##     pos_source = array([int(source['X_IMAGE']),int(source['Y_IMAGE'])])
238##     pos_mean_dist = source['MEAN_DIST_TO_NGBRS']
239##     for rel_x in range(-width,width):
240##         for rel_y in range(-width,width):
241##             rel_vec = array([rel_x,rel_y])
242##             # Forme du patch : carré
243##             if ( (dot(rel_vec,vec_i) < vec_i_2/2) and
244##                  (dot(rel_vec,vec_j) < vec_j_2/2)):
245##                 (x,y) = pos_source+rel_vec
246##                 if ( x>0 and x<2048 and y>0 and y<2048):
247##                     a[x,y] = pos_mean_dist
248
249
250
251a = lookup_sources[:,:,3]
252           
253hdu = pyfits.PrimaryHDU(a)
254hdulist = pyfits.HDUList([hdu])
255
256hdu.writeto(options.outputfits,clobber=True)
257
258# Facteur d'échelle moyen
259
260tmp = []
261for i in range(81):
262    for j in range(81):
263        if (lookup_sources[i,j,0] == 1):
264            tmp.append(lookup_sources[i,j,3])
265
266
267meas = open(options.measurement,'w')
268
269meas.write("# Moyenne des mesures /  Module vec_i /  Module vec_j / deltaT (min) / date\n")
270outs = "%7.4f %7.4f %7.4f %7.3f %s\n" % (sum(tmp)/len(tmp),
271                                         modulus(vec_i),
272                                         modulus(vec_j),
273                                         phase,
274                                         dt.isoformat())
275meas.write(outs)
276
277meas.close()
278
Note: See TracBrowser for help on using the repository browser.