source: trunk/src/file_Tbcorr_daily.py @ 694

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

# import
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

rep_fic='/homedata/mleduc/CorrTb/'
### Satellite, Instrument, Produit Level, Channel
sat  = ['M02','M01','N19','N18','N16','N15']
inst = ['MHS','MHS','MHS','MHS','AMSUB','AMSUB']
lvl  = 'L1C'
ch   = ['B1','B2','B3','B4','B5']
### Period
dmin = date(2008,5,1)
dmax = date(2008,5,3)
# Zone
#LS = -1.
#LN = 4.
#LS = 7.
#LN = 12.
#LW = -10.
#LE = 5.

# New Grid
resol = .5
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',))
        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'