Changeset 708 for trunk

Timestamp:

04/27/16 01:37:34 (8 years ago)

Author:

llmlod

Message:

Correction with spatialy smoothed coef (to be checked)

Location:

trunk/src

Files:

• file_Tbcorr_daily.py (modified) (8 diffs)
• useCorrTb_daily.py (modified) (5 diffs)

trunk/src/file_Tbcorr_daily.py

@@ -13 +13 @@
 import os
 from pylab import date2num
-from numpy import arange, size, hstack, dstack, dtype, argsort, nanmean, isnan
+from numpy import arange, size, hstack, dstack, dtype, argsort, mean, nanmean, isnan
 from dateutil.rrule import rrule, DAILY
 from useCorrTb_daily import applycorr2orb, SatName
@@ -85 +85 @@
                         type=valid_date,
                         metavar='<yyyymmdd>',
-                        default=datetime.date(2008, 5, 1),
+                        default=datetime.date(2011, 1, 3),
                         help=('first date yyyymmdd'))
 
@@ -93 +93 @@
                         type=valid_date,
                         metavar='<yyyymmdd>',
-                        default=datetime.date(2008, 5, 3),
+                        default=datetime.date(2011, 1, 3),
                         help=('last date yyyymmdd'))
 
@@ -150 +150 @@
 
     # number of processes
-    npr = 10
-
+    npr = 30
     for day in rrule(DAILY, dtstart=dmin, until=dmax):
         logger.info('{0}'.format(day))
@@ -159 +158 @@
         repsat = [ii+ss for ii, ss in zip(inst, sat)]
         fname = hstack([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 = multiprocessing.Pool(processes=npr) 
@@ -183 +182 @@
 
         # Fillvalue
-        tbcorr[isnan(tbcorr)] = -99.99
+        fvalue = -9999
+        tbcorr[isnan(tbcorr)] = fvalue
 
         # Save in netcdf file
         logger.info('Save in nc')
-        nsat = satvec.size
         nlon = lonvec.size
         nlat = latvec.size
         torig = datetime.date(1970, 1, 1)
         for cc, ich in enumerate(ch):
+            # Index to remove orbites completly 
+            iorb = (mean(tbcorr[cc,:,:].squeeze(),axis=0) != fvalue)
+            nsat = satvec[iorb].size
             # Output file name
             output = '%sTbCorr_AMSUMHS_%s_%3.2fdeg_%s.nc' % (rep_fic, ich, resol, day.strftime('%Y%m%d'))
@@ -217 +219 @@
             sats.long_name = 'Satellite'
             sats.description = [str(num)+' = '+ns +', ' for num, ns in zip(satdef.satnum, satdef.satname)]
-            sats[:] = satnum
+            sats[:] = satnum[iorb]
 
             # Creation variable
@@ -223 +225 @@
             dates.long_name = 'Time'
             dates.units = 'hours since %s'%torig.strftime('%Y-%m-%d %H:%M:%S')
-            dates[:] = ((tobs-date2num(torig))*24.).reshape(nlat, nlon, nsat)
+            dates[:] = ((tobs[:,iorb]-date2num(torig))*24.).reshape(nlat, nlon, nsat)
 
             sclfactor = .01
-            tbfinal = ncfile.createVariable('tbcorr', dtype('int16').char, ('lat', 'lon', 'satellite', ), fill_value=-9999.)
+            tbfinal = ncfile.createVariable('tbcorr', dtype('int16').char, ('lat', 'lon', 'satellite', ), fill_value=fvalue)
             tbfinal.description = 'Corrected Tb of channel %s'%ich
             tbfinal.units = 'Kelvin'
             tbfinal.scale_factor = sclfactor
-            tbfinal[:] = tbcorr[cc, :, :].reshape(nlat, nlon, nsat)
+            tbfinal[:] = tbcorr[cc, :, iorb].reshape(nlat, nlon, nsat)
 
             # Fermeture du fichier

trunk/src/useCorrTb_daily.py

@@ -8 +8 @@
 import time
 import netCDF4
-from pylab import date2num, num2date, flatten
-from numpy import array, arange, float32, NaN, transpose, int16, diff, logical_and, meshgrid, ravel, zeros, unique, nanmean
+from pylab import date2num, num2date#, flatten
+from numpy import array, arange, float32, NaN, transpose, int16, diff, logical_and, meshgrid, ravel, zeros, unique, nanmean, isnan, hstack
 from ReadRawFile import ReadRawFile
 from scipy import spatial
@@ -50 +50 @@
     logger = logging.getLogger()
 
+    if ch[0].lower() == 'all':
+        ch = ['B1', 'B2', 'B3', 'B4', 'B5']
+
     # Read orbit file
     lon, lat, dt, tb, fov = load_amsub_mhs(fname, ch, area)
 
     try:
-        # Class of each orbit points
-        clz = find_class(lat, lon)
-        nozero = (clz != 0)
-
-    except ValueError:
+        # Coef of correction as a function of classification, fov and day of year
+        coefz = find_corrcoef(dt, lon, lat, fov, ch)
+        nonan = ~isnan(coefz[0,:])
+
+    except StandardError:
         logger.error('Orbit is out of the area: {0}'.format(fname))
         return [], [], []
 
     else:
-        # Coef of correction as a function of classification, fov and day of year
-        coefz = find_corrcoef(dt[nozero], fov[nozero], clz[nozero], ch=ch)
-        tbcorr = array(tb)[:, nozero] - coefz
+        tbcorr = array(tb)[:,nonan] - coefz[:,nonan]
         # Average on new grid
-        tbcorrm, tobsm = newgrid(lonvec, latvec, lon[nozero], lat[nozero], dt[nozero], tbcorr)
+        tbcorrm, tobsm = newgrid(lonvec, latvec, lon[nonan], lat[nonan], dt[nonan], tbcorr)
 
     logger.info('Deal {0} in: {1}'.format(fname, time.time()-t))
@@ -78 +79 @@
     *Inputs:
     fname = file name
-    ch    = list of extrated channel ('B1', 'B2', 'B3', 'B4', 'B5' or 'All'), default 'All'
-    area  = limits [LS, LN, LW, LE] for extracted area, default is world (ie [-35, 35, -180., 180])
+    ch    = list of extrated channel ('B1', 'B2', 'B3', 'B4', 'B5')
+    area  = limits [LS, LN, LW, LE] for extracted area, default is Tropics (ie [-35, 35, -180., 180])
 
     *Outputs:
@@ -90 +91 @@
     """
 
-    # Index of channel
-    ch = [cch.lower() for cch in ch]
-    if ch[0] == 'all':
-        ich = arange(5)
-    else:
-        ich = array([int(cch[1])-1 for cch in ch])
+    # Index channel to extract
+    ich = array([int(cch[1])-1 for cch in ch])
 
     project = os.environ['PROJECT']
@@ -188 +185 @@
     return h, m, s
 
-# ---Find classes--- #
-def find_class(lat, lon):
-    """Find class of each given (lat, lon)
-    *Inputs:
+# ---Find coefficient of correction--- #
+def find_corrcoef(dt, lon, lat, fov, ch):
+    """Find the coefficient of correction
+    *Inputs:
+    dt  = Python date
     lat  = Latitude (degrees)
     lon  = Longitude (degrees)
-
-    *Outputs:
-    clpt = class of each point (same shape as lat and lon)
-    """
-    project = os.environ['PROJECT']
-    # Read the classification file
-    fic_class = '{0}/data/classifERAI_Trop.nc'.format(project)
-    fileID = netCDF4.Dataset(fic_class)
-    lat_class = fileID.variables['lat'][:]
-    lon_class = fileID.variables['lon'][:]
-    classif = int16(fileID.variables['classif_tot'][:])
-    fileID.close()
-
-    # Latitude limits of classification (ie tropical band)
-    latstep = abs(diff(lat_class))[0]
-    latmin = lat_class.min()-latstep/2.
-    latmax = lat_class.max()+latstep/2.
-    # Only points in the area of classification
-    inclass = logical_and(lat >= latmin, lat <= latmax)
-
-    # Geographic projection (lon/lat in degrees to x/y in m)
-    projG = pyproj.Proj("+proj=merc +ellps=WGS84 +datum=WGS84")
-    longd_class, latgd_class = meshgrid(lon_class, lat_class)
-    xclass, yclass = projG(longd_class, latgd_class)
-    xx, yy = projG(lon[inclass], lat[inclass])
-
-    # Find the nearest neighbour (euclidian distance) with a kd-tree
-    tree = spatial.cKDTree(zip(ravel(xclass), ravel(yclass)))
-    _, ind = tree.query(zip(xx, yy))
-    clpt = zeros(lat.shape, dtype=int16)
-    clpt[inclass] = classif.flatten()[ind]
-
-    return clpt
-
-# ---Find the coefficient of correction--- #
-def find_corrcoef(dt, fov, cl, ch=['All']):
-    """Find the coefficient of correction
-    *Inputs:
-    dt  = Python date
     fov = Location in the swath
-    cl  = Class
-    ch  = list of extrated channel ('B1', 'B2', 'B3', 'B4', 'B5' or 'All'), default 'All'
+    ch  = list of extrated channel ('B1', 'B2', 'B3', 'B4', 'B5')
 
     *Outputs:
     coef = Coefficient of correction
     """
-    # Index of channel
-    ch = [cch.lower() for cch in ch]
-    if ch[0] == 'all':
-        ich = arange(5)
-    else:
-        ich = array([int(cch[1])-1 for cch in ch])
-
-    # Read the classification file
-    project = os.environ['PROJECT']
-    fic_corr = '{0}/data/corrTb_AMSUB_MHS_L1C.nc'.format(project)
-    fileID = netCDF4.Dataset(fic_corr)
-    day_corr = fileID.variables['day'][:]
-    coef_corr = fileID.variables['coefcorr'][:, :, :, ich]
-    fileID.close()
 
     # Find the day of year
     yday = array([d.timetuple().tm_yday for d in dt])
-    # Find the day of correction
-    nday = day_corr.size
-    daystep = diff(day_corr)[0]
-    icorr = list(flatten([[ii]*daystep for ii in xrange(nday)])) + [0]*6
-    iday = array(icorr)[yday-1]
-
-    # Find correction
-    noclass = (cl == 0)
-    icl = array(cl-1, dtype=int)
-    ifov = array(fov-1, dtype=int)
+
     coef = []
-    for ii in ich:
-        dd = zeros(cl.shape)+NaN
-        dd[~noclass] = coef_corr[iday[~noclass], icl[~noclass], ifov[~noclass], ii]
-        coef.append(dd)
+    # Read the coefficient correction file for each channel 
+    for ich in ch:
+        fic_corr = '/homedata/mleduc/CorrTb/corrTb_AMSUB_MHS_L1C_smooth_%s.nc' % ich
+        fileID = netCDF4.Dataset(fic_corr)
+        coef_corr = fileID.variables['corrTb'][:]
+        if ich == ch[0]:
+            lat_corr = fileID.variables['lat'][:]
+            lon_corr = fileID.variables['lon'][:]
+            day_corr = fileID.variables['day'][:]
+
+            # Latitude limits of classification (ie tropical band)
+            latstep = abs(diff(lat_corr))[0]
+            latmin = lat_corr.min()-latstep/2.
+            latmax = lat_corr.max()+latstep/2.
+
+            # Geographic projection (lon/lat in degrees to x/y in m)
+            projG = pyproj.Proj("+proj=merc +ellps=WGS84 +datum=WGS84")
+            longd_corr, latgd_corr = meshgrid(lon_corr, lat_corr)
+            xcorr, ycorr = projG(longd_corr, latgd_corr)
+
+            # Kd-tree to find the nearest neighbour (euclidian distance)
+            tree = spatial.cKDTree(zip(ravel(xcorr), ravel(ycorr)))
+
+            # Find the day of correction
+            nday = day_corr.size
+            daystep = int(abs(diff(day_corr))[0])
+            icorr = [[ii]*daystep for ii in xrange(nday)] + [0]*6
+
+        fileID.close()
+        # Only points in the area of classification
+        inclass = logical_and(lat >= latmin, lat <= latmax)
+        # index for days
+        iday = hstack(icorr)[yday[inclass]-1]
+        # index for location in the swath
+        ifov = array(fov[inclass]-1, dtype=int)
+        # index for lon/lat
+        xx, yy = projG(lon[inclass], lat[inclass])
+        _, ind = tree.query(zip(xx, yy))
+        # Find coefficient of correction
+        coefch = zeros(lat.shape, dtype=float32)+NaN
+        coefch[inclass] = [coef_corr[dd, :, :, ff].flatten()[ii] for dd, ff, ii in zip(iday, ifov, ind)]
+        # Concatenate channel
+        coef.append(coefch)
 
     return float32(coef)