source: mire/process.py @ 18

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

mises au point
  • Property svn:executable set to *
File size: 9.0 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
160fen = 2
161
162shift_list = []
163
164for i in range(-fen,fen+1):
165    for j in range(-fen,fen+1):
166        if (i!=0 or j!=0):
167            shift_list.append([i,j])
168
169for i in range(fen,81-fen):
170    for j in range(fen,81-fen):
171        if (lookup_sources[i,j,0] > 0):
172            knot_idx = int(lookup_sources[i,j,2])
173            catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = 0
174            nb_voisin = 0
175            distance = 0
176            weights = 0
177            pos_source = array([catalog[knot_idx]['X_IMAGE'],
178                                catalog[knot_idx]['Y_IMAGE']])
179            for shifts in shift_list :
180                i_voisin = i + shifts[0]
181                j_voisin = j + shifts[1]
182                if (lookup_sources[i_voisin,j_voisin,0] > 0):
183                    neigh_idx = int(lookup_sources[i_voisin,j_voisin,2])
184                    nb_voisin +=1
185                    pos_voisin = array([catalog[neigh_idx]['X_IMAGE'],
186                                        catalog[neigh_idx]['Y_IMAGE']])
187                    distance += modulus(pos_source-pos_voisin)
188                    weights += modulus(shifts)
189            if (nb_voisin > 0):
190                catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = float(distance)/weights
191                lookup_sources[i,j,3] = float(distance)/weights
192       
193cleaned_catalog = []
194for i in range(81):
195    for j in range(81):
196        if (lookup_sources[i,j,0] > 0):
197            cleaned_catalog.append(catalog[int(lookup_sources[i,j,2])])
198
199# Sortie du catalogue nettoyé des doublons
200
201of = open('clean.cat','w')
202
203fields = ["FLUX_MAX",
204          "X_IMAGE",
205          "Y_IMAGE",
206          "THETA_IMAGE",
207          "ELONGATION",
208          "ELLIPTICITY",
209          "FWHM_IMAGE",
210          "ERRX2_IMAGE",
211          "ERRY2_IMAGE",
212          "I_GRID",
213          "J_GRID",
214          "KNOT_DIST",
215          "MEAN_DIST_TO_NGBRS"]
216
217fmts = ["%10.3f ",
218        "%8.3f ",
219        "%8.3f ",
220        "%8.3f ",
221        "%8.3f ",
222        "%8.3f ",
223        "%8.3f ",
224        "%8.5f ",
225        "%8.5f ",
226        "%5d ",
227        "%5d ",
228        "%5.3f ",
229        "%8.3f"]
230
231for source in cleaned_catalog:
232    outs = ""
233    for i in range(len(fields)):
234        outs += fmts[i] % source[fields[i]]
235    outs += "\n"
236    of.write(outs)
237
238of.close()
239
240# Génération d'un fichier fits contenant les résultats
241
242
243import pyfits
244
245#a = zeros((2048,2048))
246#
247# patches width in pixel:
248## width = int(modulus(vec_i))
249
250## for i in range(len(cleaned_catalog)):
251##     print i
252##     source = cleaned_catalog[i]
253##     pos_source = array([int(source['X_IMAGE']),int(source['Y_IMAGE'])])
254##     pos_mean_dist = source['MEAN_DIST_TO_NGBRS']
255##     for rel_x in range(-width,width):
256##         for rel_y in range(-width,width):
257##             rel_vec = array([rel_x,rel_y])
258##             # Forme du patch : carré
259##             if ( (dot(rel_vec,vec_i) < vec_i_2/2) and
260##                  (dot(rel_vec,vec_j) < vec_j_2/2)):
261##                 (x,y) = pos_source+rel_vec
262##                 if ( x>0 and x<2048 and y>0 and y<2048):
263##                     a[x,y] = pos_mean_dist
264
265
266
267a = lookup_sources[:,:,3]
268           
269hdu = pyfits.PrimaryHDU(a)
270hdulist = pyfits.HDUList([hdu])
271
272hdu.writeto(options.outputfits,clobber=True)
273
274# Facteur d'échelle moyen
275
276tmp = []
277for i in range(81):
278    for j in range(81):
279        if (lookup_sources[i,j,0] == 1):
280            tmp.append(lookup_sources[i,j,3])
281
282
283meas = open(options.measurement,'w')
284
285meas.write("# Moyenne des mesures /  Module vec_i /  Module vec_j / deltaT (min) / date\n")
286outs = "%7.4f %7.4f %7.4f %7.3f %s\n" % (sum(tmp)/len(tmp),
287                                         modulus(vec_i),
288                                         modulus(vec_j),
289                                         delta_dt.seconds/60.,
290                                         dt.isoformat())
291meas.write(outs)
292
293meas.close()
294
Note: See TracBrowser for help on using the repository browser.