source: trunk/src/file_Tbcorr_daily.py @ 699

#!/usr/bin/env python
# -*- coding: utf-8 -*-

# import

import argparse

import inspect
from datetime import datetime
import logging

#from pylab import *
import glob
#from datetime import date
from useCorrTb_daily import *
import itertools as it
from multiprocessing import Pool
from functools import partial
import matplotlib.gridspec as gridspec
import time


def file_Tbcorr_daily():

    """
    correct AMSU, MHS and METOP files
    """
    project_id_env = os.environ['PROJECT_ID']

    rep_fic=project_id_env

    parser = argparse.ArgumentParser(description=__doc__,
                                     formatter_class=argparse.RawDescriptionHelpFormatter)

    parser.add_argument("--verbose", 
                        action='store_true',
                        default=False,
                        help=("run in verbose mode."))

    ### Satellite, Instrument, Produit Level, Channel
    parser.add_argument('--sat',
                        dest='sat',
                        nargs="*",
                        choices=['M02','M01','N19','N18','N16','N15'],
                        metavar="<satellite>",
                        default=['M02','M01','N19','N18','N16','N15'],
                        help=("satellite"))

    parser.add_argument('--inst',
                        dest='inst',
                        nargs="*",
                        choices=['MHS','MHS','MHS','MHS','AMSUB','AMSUB'],
                        metavar="<instrument>",
                        default=['MHS','MHS','MHS','MHS','AMSUB','AMSUB'],
                        help=("instrument"))

    parser.add_argument('--level', 
                        dest='lvl',
                        action='store',
                        metavar='<level>',
                        default='L1C',
                        help=("level"))

    parser.add_argument('--ch',
                        dest='ch',
                        nargs='*',
                        choices=['B1','B2','B3','B4','B5'],
                        metavar='<channel>',
                        default=['B1','B2','B3','B4','B5'],
                        help=("channel"))

    ### Period
    def valid_date(datestring):
        """
        used to parse yyyymmdd argument in CLI
        """
        
        return datetime.datetime.strptime(datestring, '%Y%m%d')

    parser.add_argument('--dmin',
                        dest='dmin',
                        action='store',
                        type=valid_date,
                        metavar='<yyyymmdd>',
                        default=datetime.date(2008,5,1),
                        help=('first date yyyymmdd'))

    parser.add_argument('--dmax',
                        dest='dmax',
                        action='store',
                        type=valid_date,
                        metavar='<yyyymmdd>',
                        default=datetime.date(2008,5,3),
                        help=('last date yyyymmdd'))

    # Zone
    #LS = -1.
    #LN = 4.
    #LS = 7.
    #LN = 12.
    #LW = -10.
    #LE = 5.
    
    # New Grid
    resol = .5

    opts = vars(parser.parse_args())
    is_verbose = opts['verbose']
    lvl = opts['lvl']
    sat = opts['sat']
    inst = opts['inst']
    ch = opts['ch']
    dmin = opts['dmin']
    dmax = opts['dmax']

    latvec = arange(-34.5,34.6,resol)
    lonvec = arange(-180.,179.76,resol)

    # define log file
    command = '{0}'.format(inspect.stack()[0][3])
    log_date = datetime.datetime.utcnow().strftime('%Y%m%dT%H%M%SZ')
    log = '{0}/{1}.log.{2}'.format(os.environ['PROJECT_LOG'], command, log_date)

    logger = logging.getLogger()
    handler = logging.FileHandler(log)
    logger.addHandler(handler)
    if is_verbose:
        logger.setLevel(logging.DEBUG)
    else:
        logger.setLevel(logging.INFO)
    
    logger = logging.getLogger()
    logger.info('[Context]')
    logger.info('PID={0}'.format(os.getpid()))
    logger.info('command={0}'.format(command))
    import socket
    logger.info('hostname={0}'.format(socket.getfqdn()))
    logger.info('runtime={0}'.format(log_date))
    logger.info('log={0}'.format(log))

    logger.info('[Parameters]')
    logger.info('lvl={0}'.format(lvl))
    logger.info('sat={0}'.format(sat))
    logger.info('inst={0}'.format(inst))
    logger.info('ch={0}'.format(ch))
    logger.info('dmin={0}'.format(dmin))
    logger.info('dmax={0}'.format(dmax))

    # number of processes
    npr = 10
    
    for day in rrule(DAILY,dtstart=dmin,until=dmax):
        logger.info('{0}'.format(day))
        ### File list
        repday = [day.strftime('%Y_%m_%d')]#[ff.strftime('%Y_%m_%d') for ff in rrule(DAILY,dtstart=day+relativedelta(days=-1),until=day+relativedelta(days=+1))]
        repsat = [ii+ss for ii,ss in zip(inst,sat)]
        fname = hstack([glob.glob('/bdd/AMSU-1C/%s/L1C/*/%s/*.L1C'%(reps,repd)) for reps,repd in it.product(repsat,repday)])
    
        ### Read files (in parallel with Pool)
        pool = Pool(processes=npr) 
        res = pool.map(partial(applycorr2orb,lonvec=lonvec,latvec=latvec,ch=ch),fname)
        pool.close()
        pool.join()
    
        logger.info('Join')

        tbcorr, tobs, satvec = [dstack([proc[var] for proc in res if size(proc[var])]) for var in xrange(3)]
        del res
    
        # Sort ??? sur le temps ???
        isort = argsort(nanmean(tobs.squeeze(),axis=0))
        tbcorr = tbcorr[:,:,isort]
        tobs = tobs.squeeze()[:,isort]
        satvec = satvec.squeeze()[isort]
    
        # Satellite/Instrument identification
        satdef = SatName()
        satnum = satdef.get_satnum(satvec.squeeze())
    
        ### Fillvalue
        tbcorr[isnan(tbcorr)] = -99.99
    
        ### Save in netcdf file
        logger.info('Save in nc')
        nsat = satvec.size
        nlon = lonvec.size
        nlat = latvec.size
        torig = date(day.year,1,1)
        for cc, ich in enumerate(ch):
            # Output file name
            output = '%sTbCorr_AMSUMHS_%s_%3.2fdeg_%s.nc'%(rep_fic,ich,resol,day.strftime('%Y%m%d'))
            # Creation du fichier    
            ncfile = netCDF4.Dataset(output,'w',format='NETCDF4_CLASSIC')
            # Definition des dimensions
            ncfile.createDimension('lon',nlon)
            ncfile.createDimension('lat',nlat)
            ncfile.createDimension('satellite',nsat)
            
            # Creation des variables dimensions
            lons = ncfile.createVariable('lon',dtype('float32').char,('lon',))
            lons.long_name    = 'Longitude'
            lons.units        = 'degrees_east'
            lons.actual_range = '%4.2f,%4.2f'%(lonvec.min(),lonvec.max())
            lons[:]           = lonvec
        
            lats = ncfile.createVariable('lat',dtype('float32').char,('lat',))
            lats.long_name    = 'Latitude'
            lats.units        = 'degrees_north'
            lats.actual_range = '%4.2f,%4.2f'%(latvec.min(),latvec.max())
            lats[:]           = latvec
    
            sats = ncfile.createVariable('satellite',dtype('int16').char,('satellite',))
            sats.long_name    = 'Satellite'
            sats.description  = [str(num)+' = '+ns +', ' for num,ns in zip(satdef.satnum,satdef.satname)]
            sats[:]           = satnum
    
            # Creation variable
            dates = ncfile.createVariable('time',dtype('float').char,('lat','lon','satellite',))
            dates.long_name   = 'Time'
            dates.units       = 'hours since %s'%torig.strftime('%Y-%m-%d %H:%M:%S')
            #dates[:]          = (fix((tobs-date2num(torig))*24.*100)/100.).reshape(nlat,nlon,nsat)
            dates[:]          = ((tobs-date2num(torig))*24.).reshape(nlat,nlon,nsat)
    
            sclfactor = .01
            tbfinal = ncfile.createVariable('tbcorr',dtype('int16').char,('lat','lon','satellite',),fill_value=-9999.)
            tbfinal.description  = 'Corrected Tb of channel %s'%ich
            tbfinal.units        = 'Kelvin'
            tbfinal.scale_factor = sclfactor
            #tbfinal._FillValue   = -9999.
            tbfinal[:]           = tbcorr[cc,:,:].reshape(nlat,nlon,nsat)
    
            # Fermeture du fichier
            ncfile.close()
            logger.info('{0} created'.format(output))


if __name__ == '__main__':
    file_Tbcorr_daily()