source: trunk/DataCheck/Processing/fill_sqm_data_into_db.py@ 19043

#!/usr/bin/env python2 
# coding: utf-8
"""
authors: Max Ahnen, Dominik Neise
----------------------------------------------------------------------------
"THE BEER-WARE LICENSE" (Revision 42):
Max Ahnen and Dominik Neise wrote this file. As long as you retain this notice you
can do whatever you want with this stuff. If we meet some day, and you think
this stuff is worth it, you can buy us a beer.
----------------------------------------------------------------------------

This script calculates the mean of the SQM magnitude for a given run
(run is given by its start_ and stop_time)
and the p-value for a linear fit.

We hope that this will be used to separate good data, where 
the linear fit is perfect (p-value > 1e-3?) and bad fits, 
where clouds let the magnitude brightness fluctuate stronger,
so the fit worsens.

Returns: a string
  "results {mean magnitude:f} {p-value:f}"   
"""

from astropy.io import fits
import numpy as np
import scipy as sp
import sys
import ROOT 

import pandas as pd
from sqlalchemy import create_engine

import glob

from calendar import timegm
import time

database = {
    'user': 'factread',
    'password': 'r3adfac!',
    'host': '129.194.168.95',
    'table': 'factdata',
}

root_database = {
    'user': 'root',
    'password': '1440Gapd',
    'host': 'localhost',
    'table': 'factdata',
}


db_string = '{user}:{password}@{host}/{table}'

try:
    factdb = create_engine('mysql+mysqldb://'+ db_string.format(**database))
    factdb_root = create_engine('mysql+mysqldb://'+ db_string.format(**root_database))
except ImportError:
    factdb = create_engine('mysql+pymysql://'+db_string.format(**database))
    factdb_root = create_engine('mysql+pymysql://'+db_string.format(**root_database))


def get_list_of_SQM_files(base_path='/daq/aux'):
    return sorted(glob.glob(base_path+'/*/*/*/*.SQM_CONTROL_DATA.fits'))

def get_y_m_d(file_path):
    s = file_path.split('/')[-1].split('.')[0]
    return s[0:4], s[4:6], s[6:8]

def get_night(file_path):
    return file_path.split('/')[-1].split('.')[0]

def mag_mean_p_value(fits_file_path, start_time, stop_time):

        d = fits.open(fits_file_path)[1].data
        d = d[(d['Time'] > start_time) * (d['Time'] < stop_time)]
        if len(d)==0:
            return None

        x =  d['Time'].copy() - d['Time'][0]
        y = d['Mag'].astype(np.float64).copy()
        start_time, stop_time = x[0], x[-1]

        sigma  = 0.025 / np.sqrt(12.)
        g = ROOT.TGraphErrors(len(x), x, y, np.zeros(len(x), dtype=np.float64), np.ones(len(x), dtype=np.float64)*sigma)
        f = ROOT.TF1("linfit", "pol1", start_time, stop_time)
        g.Fit(f, "E")
        
        function = f.Eval
        y_fit = map(function,  x)
 
        result = {
            'fSqmMagMean' : d['Mag'].mean(),
            'fSqmMagLinFitPValue' : f.GetProb(),
            'fSqmMagLinFitChi2' : f.GetChisquare(),
            'fSqmMagLinFitNdf' : f.GetNDF(),
            'fSqmMagLinFitSlope' : f.GetParameters()[1],
        }

        return result

def timestamp_from_string(stri):
    return timegm(time.strptime(stri.replace('Z', 'UTC'),'%Y-%m-%d %H:%M:%S%Z'))

def update_dict_in_database(some_dict, db, primary_keys=("fNight", "fRunID")):
    commands = []
    commands.append('BEGIN;')

    # "UPDATE RunInfo SET weight = 160, desiredWeight = 145 WHERE id = 1;"
    update_string = "UPDATE RunInfo SET "
    first = True
    for k in some_dict:
        if k not in primary_keys:
            if not first:
                update_string += ', '
            else:
                first = False

            update_string += "{0} = {1}".format(k, some_dict[k])
    update_string += " WHERE "
    first = True
    for k in some_dict:
        if k in primary_keys:
            if not first:
                update_string += ' AND '
            else:
                first = False
            update_string += "{0} = {1}".format(k, some_dict[k])
    update_string += ';'

    commands.append(update_string)
    commands.append('COMMIT;')

    # print commands
    for com in commands:
        db.engine.execute(com)

    
if __name__ == "__main__":
    aux_files = get_list_of_SQM_files()
    for aux_file in aux_files:
      night = get_night(aux_file)

      query = ("SELECT fRunID, fRunStart, fRunStop from RunInfo"
               " WHERE fNight={0}").format(night)

      df = pd.read_sql_query(query, factdb)

      for i in range(df.shape[0]):
          row = df.iloc[i]

          try:
              run_start = timestamp_from_string(str(row['fRunStart'])+'Z')/(24.*3600.)
              run_stop = timestamp_from_string(str(row['fRunStop'])+'Z')/(24.*3600.)
              run_id = row['fRunID']
          except ValueError as e:
              print e
              print row['fRunStart'], row['fRunStop']
              continue

          result= mag_mean_p_value(aux_file, run_start, run_stop)
          if result is None:
              continue
          result['fNight'] = int(night)
          result['fRunID'] = run_id

          update_dict_in_database(result, factdb_root)

          print time.asctime(), aux_file, run_start, run_stop, result