source: trunk/src/file_Tbcorr_daily.py @ 698

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

# import

import argparse

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)

    ### 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())
    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)
    
    # number of processes
    npr = 10
    
    for day in rrule(DAILY,dtstart=dmin,until=dmax):
        print 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()
    
        print '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
        print '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()
            print output+' created'


if __name__ == '__main__':
    file_Tbcorr_daily()