source: trunk/src/useCorrTb_daily.py @ 705

Last change on this file since 705 was 705, checked in by pinsard, 8 years ago

no more /homedata/mleduc/CorrTb/, needed on ciclad
File size: 10.0 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3
4# import
5import logging
6
7from pylab import *
8from ReadRawFile import *
9from datetime import date
10import netCDF4
11from scipy import spatial
12from mpl_toolkits.basemap import pyproj
13import itertools as it
14import time
15
16##-------------------------------##
17#     Class for satellite name    #
18##-------------------------------##
19class SatName: 
20    """ ### Defined a number for each satellite """
21    def __init__(self):
22        self.satname = array(['AMSUBN15','AMSUBN16','MHSN18','MHSN19','MHSM01','MHSM02'])
23        self.satnum  = arange(size(self.satname))+1
24
25    def get_satnum(self,name):
26        return concatenate([self.satnum[nn == self.satname] for nn in name])
27   
28    def get_satname(self,num):
29        return concatenate([self.satname[nn == self.satnum] for nn in num])
30
31##-------------------------------##
32#     Correction for one orbit    #
33##-------------------------------##
34def applycorr2orb(fname,lonvec,latvec,ch=['All'],area=[-35.,35.,-180.,180.]):
35    """ ### Apply correction to AMSU or MHS orbit
36    ## Inputs
37    # fname  = file name
38    # tstep  = time resolution for new grid in minutes
39    # lonvec = longitude vecteur of new grid
40    # latvec = latitude vecteur of new grid
41    # ch     = list of extrated channel ('B1','B2','B3','B4','B5' or 'All'), default 'All'
42    # area   = limits [LS,LN,LW,LE] for extracted area, default is world (ie [-35.,35.,-180.,180.])
43    #
44    # Outputs
45    # tbcorrm = brightness temperature of orbit in new grid
46    # trefm   = ordinal time
47    #
48    """
49
50    t=time.time()
51
52    logger = logging.getLogger()
53
54    # ## Read orbit file
55    # print 'Read ', fname
56    lon, lat, dt, tb, fov = load_amsub_mhs(fname,ch,area)
57
58    try:
59        # ##  Class of each orbit points
60        # print 'Classification'
61        clz = find_class(lat,lon)
62        nozero = (clz != 0)
63
64    except ValueError:
65        logger.error('Orbit is out of the area: {0}'.format(fname))
66        return [], [], []
67
68    else:
69        # ## Coef of correction as a function of classification, fov and day of year
70        # print 'Correction'
71        coefz = find_corrcoef(dt[nozero],fov[nozero],clz[nozero],ch=ch)
72        tbcorr = array(tb)[:,nozero] - coefz
73        # ## Average on new grid
74        tbcorrm, tobsm = newgrid(lonvec,latvec,lon[nozero],lat[nozero],dt[nozero],tbcorr)
75   
76    logger.info('Deal {0} in: {1}'.format(fname, time.time()-t))
77
78    return tbcorrm, tobsm, fname.split('/')[3]
79
80##----------------------------------##
81#      Read AMSU-B and MHS files     #
82##----------------------------------##
83def load_amsub_mhs(fname,ch,area):
84    """ ### Read AMSU-B and MHS files
85    ## Inputs
86    # fname = file name
87    # ch    = list of extrated channel ('B1','B2','B3','B4','B5' or 'All'), default 'All'
88    # area  = limits [LS,LN,LW,LE] for extracted area, default is world (ie [-35,35,-180.,180])
89    #
90    ## Outputs
91    # lon  = Longitude
92    # lat  = Latitude
93    # yday = Day of year
94    # tms  = UTC time in milliseconds
95    # tb   = Brightness temperature for each channel
96    # err  = 1 if reading error, 0 else
97    #
98    """
99
100    # Index of channel
101    ch=[cch.lower() for cch in ch]
102    if (ch[0] == 'all'):
103        ich = arange(5)
104    else:
105        ich = array([int(cch[1])-1  for cch in ch])
106
107    project = os.environ['PROJECT']
108    hdr_file = '{0}/src/external_python/ReadRawFile/data/amsub-mhs_1c_header'.format(project)
109    rec_file = '{0}/src/external_python/ReadRawFile/data/amsub-mhs_1c_record'.format(project)
110    # Variables
111    dset_names = ['amb1c_latlon','amb1c_btemps','amb1c_scnlindy','amb1c_scnlintime','amb1c_scnlinyr']
112    try:
113        headers, datasets = ReadRawFile(fname, hdr_file, rec_file, False, dset_names)   
114    except:
115        logger.error('Error with {0}'.format(fname))
116        # continue
117        return 
118
119    # Latitude
120    lat = float32(datasets['amb1c_latlon'][:,:,0]/1.E4)
121    # Longitude
122    lon = float32(datasets['amb1c_latlon'][:,:,1]/1.E4)
123    # Year
124    yyyy= datasets['amb1c_scnlinyr']
125    # Day of Year
126    yday= datasets['amb1c_scnlindy']
127    # Time in milliseconds
128    tms = datasets['amb1c_scnlintime']
129    # Tb
130    tb  = float32(datasets['amb1c_btemps'][:,:,ich]/1.E2)
131    # Remplace fill values -9999.99 by NaN
132    tb[(tb == -9999.99)] = NaN
133
134    ### Extract Area
135    # Index in the defined area
136    LS, LN, LW, LE = area
137    zone = (lon>=LW) & (lon<=LE) & (lat>=LS) & (lat<=LN)
138    nscn, nfov, nch = tb.shape
139    # Extraction of Tb for each channel
140    tb = [(tb[:,:,ii])[zone] for ii in xrange(nch)]
141    # day of year and time
142    yyyy = transpose([yyyy]*nfov)[zone]
143    yday = transpose([yday]*nfov)[zone]
144    tms  = transpose([tms]*nfov)[zone]
145    # Find location in fov
146    fov = array(range(nfov)*nscn).reshape(nscn,nfov)[zone]+1
147
148    ### Convert in python date/time
149    dd = num2date(yday)
150    tt = [ms2time(ms) for ms in tms]
151    dt = array([d.replace(year=y,hour=h,minute=m,second=s) for d,y,(h,m,s) in zip(dd,yyyy,tt)])
152
153    return lon[zone], lat[zone], dt, tb, fov
154
155##-----------------------------------------------------##
156#      Convert milliseconds in hour, minute, second     #
157##-----------------------------------------------------##
158def ms2time(ms):
159    """ ### Convert milliseconds in hour, minute, second
160    # Warning 1: lost millisecond because second is integer
161    # Warning 2: no day (but hour might be greater than 24)
162
163    ## Inputs
164    # ms  = millisecond
165    #
166    ## Output
167    # h = hour
168    # m = minute
169    # s = second
170    #
171    >>> # one hour
172    >>> ms = 1000 * 60 * 60
173    >>> (h, m, s) = ms2time(ms)
174    >>> print ('{0} h {1} m {2} s'.format(h, m, s))
175    1 h 0 m 0 s
176
177    >>> # one hour + 1 ms
178    >>> ms = (1000 * 60 * 60) + 1
179    >>> (h, m, s) = ms2time(ms)
180    >>> print ('{0} h {1} m {2} s'.format(h, m, s))
181    1 h 0 m 0 s
182
183    >>> # one day + one hour
184    >>> ms = (1000 * 60 * 60 * 24) + (1000 * 60 * 60)
185    >>> (h, m, s) = ms2time(ms)
186    >>> print ('{0} h {1} m {2} s'.format(h, m, s))
187    25 h 0 m 0 s
188    """
189
190    s    = ms/1000
191    m, s = divmod(s,60)
192    h, m = divmod(m,60)
193    #d, h = divmod(h,24)
194    return h, m, s
195
196##---------------------##
197#      Find classes     #
198##---------------------##
199def find_class(lat,lon):
200    """ ### Find class of each given (lat,lon)
201    ## Inputs
202    # lat  = Latitude (degrees)
203    # lon  = Longitude (degrees)
204    #
205    ## Output
206    # clpt = class of each point (same shape as lat and lon)
207    #
208    """
209    project = os.environ['PROJECT']
210    # Read the classification file
211    fic_class = '{0}/data/classifERAI_Trop.nc'.format(project)
212    fileID = netCDF4.Dataset(fic_class)
213    lat_class = fileID.variables['lat'][:]
214    lon_class = fileID.variables['lon'][:]
215    classif   = int16(fileID.variables['classif_tot'][:])
216    fileID.close()
217
218    # Latitude limits of classification (ie tropical band)
219    latstep = abs(diff(lat_class))[0]
220    latmin = lat_class.min()-latstep/2.
221    latmax = lat_class.max()+latstep/2.
222    # Only points in the area of classification
223    inclass = logical_and(lat>=latmin,lat<=latmax)
224
225    # Geographic projection (lon/lat in degrees to x/y in m)
226    projG = pyproj.Proj("+proj=merc +ellps=WGS84 +datum=WGS84")
227    longd_class,latgd_class = meshgrid(lon_class,lat_class)
228    xclass,yclass = projG(longd_class,latgd_class)
229    xx,yy = projG(lon[inclass],lat[inclass]) 
230
231    # Find the nearest neighbour (euclidian distance) with a kd-tree
232    tree = spatial.cKDTree(zip(ravel(xclass),ravel(yclass)))
233    _, ind = tree.query(zip(xx,yy))
234    clpt = zeros(shape(lat),dtype=int16)
235    clpt[inclass] = classif.flatten()[ind]
236
237    return clpt
238
239##-------------------------------------------##
240#      Find the coefficient of correction     #
241##-------------------------------------------##
242def find_corrcoef(dt,fov,cl,ch=['All']):
243    """ ### Find the coefficient of correction
244    ## Inputs
245    # dt  = Python date
246    # fov = Location in the swath
247    # cl  = Class
248    # ch  = list of extrated channel ('B1','B2','B3','B4','B5' or 'All'), default 'All'
249    #
250    ## Outputs
251    # coef = Coefficient of correction
252    #
253    """
254    # Index of channel
255    ch=[cch.lower() for cch in ch]
256    if (ch[0] == 'all'):
257        ich = arange(5)
258    else:
259        ich = array([int(cch[1])-1  for cch in ch])
260    nch = ich.size
261   
262    # Read the classification file
263    project = os.environ['PROJECT']
264    fic_corr = '{0}/data/corrTb_AMSUB_MHS_L1C.nc'.format(project)
265    fileID = netCDF4.Dataset(fic_corr)
266    day_corr  = fileID.variables['day'][:]
267    coef_corr = fileID.variables['coefcorr'][:,:,:,ich]#[:]
268    fileID.close()
269
270    # Find the day of year
271    yday = array([d.timetuple().tm_yday for d in dt])
272    # Find the day of correction
273    nday = day_corr.size
274    daystep = diff(day_corr)[0]
275    icorr = list(flatten([[ii]*daystep for ii in xrange(nday)])) + [0]*6
276    iday = array(icorr)[yday-1]
277
278    # Find correction
279    noclass = (cl == 0)
280    icl  = array(cl-1,dtype=int)
281    ifov = array(fov-1,dtype=int)
282    coef = []
283    for ii in ich:
284        dd = zeros(cl.shape)+NaN
285        dd[~noclass] = coef_corr[iday[~noclass],icl[~noclass],ifov[~noclass],ii]
286        coef.append(dd)
287
288    return float32(coef)
289
290##---------------------------##
291#      New grid and time      #
292##---------------------------##
293def newgrid(lonvec,latvec,lon,lat,dt,tb):
294    """ ### Find the coefficient of correction
295    ## Inputs
296    #
297    #
298    ## Outputs
299    #
300    #
301    """
302
303    # from date to ordinal
304    dtord = date2num(dt)
305
306    ### Find index of each pixel in new grid 
307    # Geographic projection (lon/lat in degrees to x/y in m)
308    projG = pyproj.Proj("+proj=merc +ellps=WGS84 +datum=WGS84")
309    longd, latgd = meshgrid(lonvec,latvec)
310    xgd, ygd = projG(longd,latgd)
311    xx, yy = projG(lon,lat) 
312    # Find the nearest neighbour (euclidian distance) with a kd-tree
313    tree = spatial.cKDTree(zip(ravel(xgd),ravel(ygd)))
314    _, indll = tree.query(zip(xx,yy))
315
316    # Mean
317    nch = tb.shape[0]
318    tbm = zeros((nch,lonvec.size*latvec.size),dtype=float32)+NaN
319    dtm = zeros((lonvec.size*latvec.size))+NaN
320    for ii in unique(indll):
321        tbm[:,ii] = nanmean(tb[:,(indll==ii)],axis=1)
322        dtm[ii]  = nanmean(dtord[(indll==ii)])
323
324    return tbm, dtm
Note: See TracBrowser for help on using the repository browser.