source: mire/process.py @ 17

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

Intégration timing

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