source: mire/process.py @ 32

#!/bin/env python
# coding: utf-8
# Un programme pour exploiter les données des mires
# -- Frédéric Meynadier

import string
from numpy import *

from functions import *

from datetime import *


from optparse import OptionParser
parser = OptionParser()
parser.add_option("-c","--catalog",dest="catalogfile",
                  help="Name of the catalog file")
parser.add_option("-r","--reference",dest="reference",
                  help="Name of the file with reference coordinates")
parser.add_option("-o","--outputfits",dest="outputfits",
                  default="output.fits",
                  help="Name of the output fits file")
parser.add_option("-m","--measurement",dest="measurement",
                  default="output.mes",
                  help="Name of the output measurement file")
parser.add_option("--ring",dest="ring",
                  default="output.ring",
                  help="Name of the output measurement file")
parser.add_option("-s","--square",dest="square",
                  default="output.square",
                  help="Name of the output measurement file")
parser.add_option("-t","--timing",dest="timing",
                  help="Name of the input timing file")


(options, args) = parser.parse_args ()

# Recherche de la date dans un fichier des timings
# (généré par recup_date, soit :
# ls -l $1| awk '{print$8 "\t" $7 "\t" $6}' | sort -k2 | sed -e 's/*//'
# où $1 est le répertoire contenant les .dat non modifiés.)
# format attendu : nom_fichier phase_cycle date_isoformat

for line in file(options.timing):
    line_s = string.split(line)
    if (len(line_s) > 2):
        nom = line_s[0] [:-4]
        if (nom == options.catalogfile[:-6]):
            dt = datetime.strptime(line_s[2],
                                   "%Y-%m-%dT%H:%M:%S")
            phase = float(line_s[1])

# Donner les coordonnées approximatives des 4 repÚres

rep = zeros((4,2))
i = 0
for line in open(options.reference):
    if ((line[0] != "#") and (line != '\n')):
        spl = string.split(line)
        rep[i][0] = float(spl[0])
        rep[i][1] = float(spl[1])
        i+=1

# Lecture du fichier catalogue
fields = ["FLUX_MAX",
          "X_IMAGE",
          "Y_IMAGE",
          "THETA_IMAGE",
          "ELONGATION",
          "ELLIPTICITY",
          "FWHM_IMAGE",
          "ERRX2_IMAGE",
          "ERRY2_IMAGE"]

catalog = []
for line in open(options.catalogfile):
    if ((line[0] != "#") and (line != '\n')):
        tmp = {}
        spl = string.split(line)
        for i in range(len(fields)):
            tmp[fields[i]] = float(spl[i])
        tmp['MEAN_DIST_TO_NGBRS']=0
        catalog.append(tmp)


mode_degrade = 0
# Recherche de la position des repÚres selon catalogue - attention aux doublons
real_rep = zeros((4,2))
for i in range(4):
    detections = 0
    for source in catalog:
        pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])
        max_dist = 5 # pixels
        if ((modulus(rep[i,:] - pos_source)) < max_dist):
            max_dist = modulus(rep[i,:] - pos_source)
            detections +=1
            real_rep[i,:] = pos_source
    if (detections > 1):
        print "Attention ! ambiguïté pour repÚre ",i,\
              ", nb de détections :",detections
        
    if (detections ==0):
        print "Erreur ! Pas de détection pour repÚre ",i
        real_rep[i,:]=rep[i,:]
        mode_degrade = 1

if (mode_degrade==1):
    print "Passage en mode dégradé, sorties réduites"

# Détermination du centre : barycentre des 4 repÚres
center = array([sum(real_rep[:,0])/4,sum(real_rep[:,1])/4])

# Vecteurs unitaires de référence
vec_i = ((real_rep[0,:] - real_rep[1,:])/20 \
         + (real_rep[3,:] - real_rep[2,:])/20) / 2 
vec_j = ((real_rep[1,:] - real_rep[2,:])/20 \
         + (real_rep[0,:] - real_rep[3,:])/20) / 2 

# Détection des noeuds de grille les plus proches de chaque source
vec_i_2 = dot(vec_i,vec_i)
vec_j_2 = dot(vec_j,vec_j)
for source in catalog:
    pos_source = array([source['X_IMAGE'],source['Y_IMAGE']])-center
    # dans le référentiel grille (en projetant sur les vecteurs unitaires) :
    pos_source_grid = array([dot(pos_source,vec_i)/vec_i_2,
                             dot(pos_source,vec_j)/vec_j_2])
    # on arrondit au noeud le plus proche
    knot = pos_source_grid.round()
    # et on calcule la distance entre les deux
    r = modulus(pos_source_grid - knot)
    # enregistrement des résultats
    source['I_GRID'] = knot[0]
    source['J_GRID'] = knot[1]
    source['KNOT_DIST'] = r
    #if (modulus(knot) <1):
    #    print source['I_GRID'], source['J_GRID'], source['KNOT_DIST'], pos_source_grid, pos_source,source['X_IMAGE'], source['Y_IMAGE']

# Chasse aux doublons
# tableau 81*81 contenant :
# champ 0 : nb de sources rattachée au noeud
# champ 1 : distance actuelle entre le noeud et la source
# champ 2 : index courant de la source dans catalog
# champ 3 : distance moyenne aux voisins (calculé plus tard)
lookup_sources = zeros((81,81,4))
for si in range(len(catalog)):
    source = catalog[si]
    i = source['I_GRID']+40
    j = source['J_GRID']+40
    if ((lookup_sources[i,j,0] < 1)
        or (lookup_sources[i,j,1] > source['KNOT_DIST'])):
        lookup_sources[i,j,2] = si 
        lookup_sources[i,j,1] = source['KNOT_DIST']
    lookup_sources[i,j,0] += 1

# On élimine le noeud (1,1), trop prÚs du "P"
lookup_sources[41,41,0] = 0

# Calcul de la distance moyenne entre noeuds voisins

fen = 2

shift_list = []

for i in range(-fen,fen+1):
    for j in range(-fen,fen+1):
        if (i!=0 or j!=0):
            shift_list.append([i,j])

for i in range(fen,81-fen):
    for j in range(fen,81-fen):
        if (lookup_sources[i,j,0] > 0):
            knot_idx = int(lookup_sources[i,j,2])
            catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = 0
            nb_voisin = 0
            distance = 0
            weights = 0
            pos_source = array([catalog[knot_idx]['X_IMAGE'],
                                catalog[knot_idx]['Y_IMAGE']])
            for shifts in shift_list :
                i_voisin = i + shifts[0]
                j_voisin = j + shifts[1]
                if (lookup_sources[i_voisin,j_voisin,0] > 0):
                    neigh_idx = int(lookup_sources[i_voisin,j_voisin,2])
                    nb_voisin +=1
                    pos_voisin = array([catalog[neigh_idx]['X_IMAGE'],
                                        catalog[neigh_idx]['Y_IMAGE']])
                    distance += modulus(pos_source-pos_voisin)
                    weights += modulus(shifts)
            if (nb_voisin > 0):
                catalog[knot_idx]['MEAN_DIST_TO_NGBRS'] = float(distance)/weights
                lookup_sources[i,j,3] = (float(distance)/weights)
        
cleaned_catalog = []
for i in range(81):
    for j in range(81):
        if (lookup_sources[i,j,0] > 0):
            cleaned_catalog.append(catalog[int(lookup_sources[i,j,2])])

# Sortie du catalogue nettoyé des doublons

of = open('clean.cat','w')

fields = ["FLUX_MAX",
          "X_IMAGE",
          "Y_IMAGE",
          "THETA_IMAGE",
          "ELONGATION",
          "ELLIPTICITY",
          "FWHM_IMAGE",
          "ERRX2_IMAGE",
          "ERRY2_IMAGE",
          "I_GRID",
          "J_GRID",
          "KNOT_DIST",
          "MEAN_DIST_TO_NGBRS"]

fmts = ["%10.3f ",
        "%8.3f ",
        "%8.3f ",
        "%8.3f ",
        "%8.3f ",
        "%8.3f ",
        "%8.3f ",
        "%8.5f ",
        "%8.5f ",
        "%5d ",
        "%5d ",
        "%5.3f ",
        "%8.3f"]

for source in cleaned_catalog:
    outs = ""
    for i in range(len(fields)):
        outs += fmts[i] % source[fields[i]]
    outs += "\n"
    of.write(outs)

of.close()

# Génération d'un fichier fits contenant les résultats


import pyfits

#a = zeros((2048,2048))
#
# patches width in pixel:
## width = int(modulus(vec_i))

## for i in range(len(cleaned_catalog)):
##     print i
##     source = cleaned_catalog[i]
##     pos_source = array([int(source['X_IMAGE']),int(source['Y_IMAGE'])])
##     pos_mean_dist = source['MEAN_DIST_TO_NGBRS']
##     for rel_x in range(-width,width):
##         for rel_y in range(-width,width):
##             rel_vec = array([rel_x,rel_y])
##             # Forme du patch : carré
##             if ( (dot(rel_vec,vec_i) < vec_i_2/2) and
##                  (dot(rel_vec,vec_j) < vec_j_2/2)):
##                 (x,y) = pos_source+rel_vec
##                 if ( x>0 and x<2048 and y>0 and y<2048):
##                     a[x,y] = pos_mean_dist



a = lookup_sources[:,:,3]

for i in range(81):
    for j in range(81):
        if (a[i,j] > 0):
            a[i,j] = 60.0/a[i,j]
        else :
            a[i,j] = 1.06

# Interpolation de la case du "P"

a[41,41] = (a[42,41] + a[40,41] + a[41,42] + a[41,40]) /4.

hdu = pyfits.PrimaryHDU(a)
hdulist = pyfits.HDUList([hdu])

hdu.writeto(options.outputfits,clobber=True)


# Facteur d'échelle moyen

tmp = []
for i in range(81):
    for j in range(81):
        if (lookup_sources[i,j,0] == 1):
            tmp.append(lookup_sources[i,j,3])


meas = open(options.measurement,'w')

meas.write("# Moyenne des mesures /  Module vec_i /  Module vec_j / deltaT (min) / date\n")
# transfo pixels / arcminutes -> mas / pixel
moyenne_mesures = 60 * (1. / (sum(tmp)/len(tmp)))
module_veci = 60 * (1. /  modulus(vec_i))
module_vecj = 60 * (1. /  modulus(vec_j))

outs = "%10.7f %10.7f %10.7f %7.3f %s\n" % (moyenne_mesures,
                                            module_veci,
                                            module_vecj,
                                            phase,
                                            dt.isoformat())
meas.write(outs)

meas.close()

## Facteur d'échelle par rapport au centre du champ sur une couronne
## représentative du diamÚtre solaire

# Quel noeud est (à peu prÚs) au centre ?


dist_mini = 2048
i_centre = 0
for i in range(len(cleaned_catalog)):
    source = cleaned_catalog[i]
    dist = modulus([source["X_IMAGE"] - 1024, source["Y_IMAGE"] - 1024])
    if (dist < dist_mini):
        dist_mini = dist
        i_centre = i

#print "Source au centre : ",cleaned_catalog[i_centre]

# On parcourt une couronne de noeuds et on calcule le facteur d'échelle pour chacun

rmin = 14 # en inter-noeud
rmax = 18 # en inter-noeud

source_centrale = cleaned_catalog[i_centre]

couronne=[]

for source in cleaned_catalog:
    tmp = {}
    dist_knots = modulus([source["I_GRID"] - source_centrale["I_GRID"],
                          source["J_GRID"] - source_centrale["J_GRID"]])
    x_rel = source["X_IMAGE"] - source_centrale["X_IMAGE"]
    y_rel = source["Y_IMAGE"] - source_centrale["Y_IMAGE"]
    dist_pix = modulus([x_rel,y_rel])
    if ((dist_knots >= rmin) and (dist_knots <= rmax)):
        tmp['ScFact'] = dist_knots * 60. / dist_pix # arcsec / pix
        tmp['r'] = dist_pix
        tmp['theta'] = math.atan(y_rel/x_rel)
        if (x_rel < 0):
            tmp['theta']+=math.pi
        if (tmp['theta'] < 0):
            tmp['theta'] += 2*math.pi
        if (tmp['theta'] >= 2*math.pi):
            tmp['theta'] -= 2*math.pi
        couronne.append(tmp)

output_couronne = file(options.ring,'w')

outs = "# phase: %f date: %s\n" % (phase,
                                 dt.isoformat())
output_couronne.write(outs)

for source in couronne:
    outs = "%f %f %10.7f\n" % (source['r'],
                               math.degrees(source['theta']),
                               source['ScFact'])
    output_couronne.write(outs)



output_couronne.close()
        
# On parcourt un carré de surface équivalente à la couronne définie ci dessus
# Le facteur d'échelle est calculé par rapport au point opposé du carré
# soit un cÎté de 28 noeuds


# Le centre a été trouvé pour la couronne
i_cen = cleaned_catalog[i_centre]["I_GRID"]
j_cen = cleaned_catalog[i_centre]["J_GRID"]

# On se refait un "lookup" simplifié

synthese = zeros((81,81,2))

for source in cleaned_catalog:
    synthese[source['I_GRID'],source['J_GRID'],0] = source['X_IMAGE']
    synthese[source['I_GRID'],source['J_GRID'],1] = source['Y_IMAGE']

valeurs =  [] 
for i in range(-14,15):
    vecteur = synthese[i_cen+14,j_cen+i,:] - synthese[i_cen-14,j_cen+i,:]
    valeurs.append(modulus(vecteur))

valeurs_orth = []
for i in range(-14,15):
    vecteur = synthese[i_cen+i,j_cen+14,:] - synthese[i_cen+i,j_cen-14,:]
    valeurs_orth.append(modulus(vecteur))

for i in range(len(valeurs)):
    if valeurs[i] > 0:
        valeurs[i] = 28*60./valeurs[i]

for i in range(len(valeurs_orth)):
    if valeurs_orth[i] > 0:
        valeurs_orth[i] = 28*60./valeurs_orth[i]


# on vire les valeurs aberrantes
def f(x):
    return ((x < 1.07) and (x > 1.05))

valeurs = filter(f, valeurs)
valeurs_orth = filter(f, valeurs_orth)

output_square = open(options.square,'w')
outs = "# moyennes: %f / %f phase: %f date: %s\n" % (mean(valeurs),
                                                     mean(valeurs_orth),
                                                     phase,
                                                     dt.isoformat())
output_square.write(outs)

for valeur in valeurs:
    outs = "%f\n" % (valeur)
    output_square.write(outs)

output_square.write('orth \n')

for valeur in valeurs_orth:
    outs = "%f\n" % (valeur)
    output_square.write(outs)

output_square.close()