Changeset 44


Ignore:
Timestamp:
07/23/14 19:15:22 (10 years ago)
Author:
lahlod
Message:

modifs

Location:
trunk/src/scripts_Laura/ARCTIC/Travail_CEN
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/calcul_spec_lamb_nadir.py

    r41 r44  
    2626for imo in range (0, M): 
    2727    ####### open .dat file to stack data ####### 
    28     data_spec_lamb = open ('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA23.dat', 'a') 
     28    data_spec_lamb = open ('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA89.dat', 'a') 
    2929    ############################################ 
    3030    print 'read file' + month[imo] 
     
    6969    fichier.close() 
    7070    for i in range (0, nbtotal): 
    71         if (emis[0, i] < 0.): 
    72             emis[0, i] = NaN  
    73         if (emis[0, i] > 1.): 
    74             emis[0, i] = 0.99   
     71        if (emis[3, i] < 0.): 
     72            emis[3, i] = NaN  
     73        if (emis[3, i] > 1.): 
     74            emis[3, i] = 0.99   
    7575    ############################ 
    7676    # definition of parameters # 
     
    8383    zen_angle = zen[bblat][bbzen] # zenith angle 
    8484    field_view = fov[bblat][bbzen] # field of view 
    85     e = emis[0, :][bblat][bbzen] # emissivity 
     85    e = emis[3, :][bblat][bbzen] # emissivity 
    8686    t_surf = ts[bblat][bbzen] # surface temperature 
    87     t_up = tup[0, :][bblat][bbzen] # tup 
    88     t_dn = tdn[0, :][bblat][bbzen] # tdown 
    89     transm = trans[0, :][bblat][bbzen] # transmissivity 
    90     t_b = tb[0, :][bblat][bbzen] # brightness temperature 
     87    t_up = tup[3, :][bblat][bbzen] # tup 
     88    t_dn = tdn[3, :][bblat][bbzen] # tdown 
     89    transm = trans[3, :][bblat][bbzen] # transmissivity 
     90    t_b = tb[3, :][bblat][bbzen] # brightness temperature 
    9191    L = len(llon) 
    9292    ####################################### 
  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/discrim_rate_espec-elamb_espec.py

    r40 r44  
    1414import map_ffgrid 
    1515from matplotlib import colors 
     16import map_cartesian_grid 
    1617 
    1718     
     
    2425 
    2526 
    26 ################### 
    27 # grid parameters # 
    28 ################### 
    29 dx=1. 
    30 dy=0.5 
    31 x0, x1 = -180., 180. 
    32 y0, y1 = 30., 90. 
    33 xvec = np.arange(x0, x1 + dx, dx) 
    34 yvec = np.arange(y0, y1 + dy, dy) 
    35 nx = len(xvec) 
     27 
     28######################## 
     29# grid characteristics # 
     30######################## 
     31x0 = -3000. # min limit of grid 
     32x1 = 2500. # max limit of grid 
     33dx = 40. 
     34xvec = np.arange(x0, x1+dx, dx) 
     35nx = len(xvec)  
     36y0 = -3000. # min limit of grid 
     37y1 = 3000. # max limit of grid 
     38dy = 40. 
     39yvec = np.arange(y0, y1+dy, dy) 
    3640ny = len(yvec) 
    37  
    38  
    3941 
    4042 
     
    4446############################################ 
    4547sl_rate = np.zeros([M, 31, ny, nx], float) 
    46 type_map = np.zeros([M, ny, nx], float) 
    4748for imo in range (0, M): 
    4849    print 'month: ' + month[imo] 
    49     fichier = Dataset('/mma/hermozol/Documents/Data/daily_OSISAF_SPEC_LAMB/daily_OSISAF_emis_spec_lamb_AMSUB_ffgrid_' + str(month[imo]) + '_2009.nc', 'r', format='NETCDF3_CLASSIC') 
    50     lon = fichier.variables['longitude'][:] 
    51     lat = fichier.variables['latitude'][:] 
     50    ### AMSUA23 ### 
     51    print 'open data file' 
     52    fichier = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA89_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
     53    xdist = fichier.variables['longitude'][:] 
     54    ydist = fichier.variables['latitude'][:] 
    5255    day = fichier.variables['days'][:] 
    53     osi_type = fichier.variables['osi_ice_type'][:] 
    54     type_map[imo, :, :] = osi_type[0, :, :] 
    55     amsub_spec = fichier.variables['spec_emis'][:] 
    56     amsub_lamb = fichier.variables['lamb_emis'][:] 
     56    emis_spec = fichier.variables['e_spec'][:] 
     57    emis_lamb = fichier.variables['e_lamb'][:] 
    5758    fichier.close() 
    5859    print 'calculate rate' 
    5960    for ijr in range (0, month_day[imo]): 
    60         for ilat in range (0, len(lat)): 
    61             for ilon in range (0, len(lon)): 
    62                 if (amsub_spec[ijr, ilat, ilon] == 0.): 
    63                     amsub_spec[ijr, ilat, ilon] = NaN 
    64                 if (amsub_lamb[ijr, ilat, ilon] == 0.): 
    65                     amsub_lamb[ijr, ilat, ilon] = NaN 
    66                 sl_rate[imo, ijr, ilat, ilon] = ((amsub_lamb[ijr, ilat, ilon] - amsub_spec[ijr, ilat, ilon]) / amsub_spec[ijr, ilat, ilon]) * 100. 
     61        for ilat in range (0, len(ydist)): 
     62            for ilon in range (0, len(xdist)): 
     63                sl_rate[imo, ijr, ilat, ilon] = ((emis_lamb[ilat, ilon, ijr] - emis_spec[ilat, ilon, ijr]) / emis_spec[ilat, ilon, ijr]) * 100. 
    6764     
    6865 
     
    7269# compute monthly means of emissivity ratio # 
    7370############################################# 
     71print 'compute monthly mean of emissivity ratio' 
    7472monthly_ratio = np.zeros([M, ny, nx], float) 
    7573for imo in range (0, M): 
    76     for ilon in range (0, len(lon)): 
    77         for ilat in range (0, len(lat)): 
     74    print 'month ' + str(month[imo]) 
     75    for ilon in range (0, len(xdist)): 
     76        for ilat in range (0, len(ydist)): 
    7877            monthly_ratio[imo, ilat, ilon] = mean(sl_rate[imo, 0 : month_day[imo], ilat, ilon][nonzero(isnan(sl_rate[imo, 0 : month_day[imo], ilat, ilon]) == False)]) 
    7978     
    8079 
     80     
     81 
     82     
    8183 
    8284######################################### 
    8385# map monthly means of emissivity ratio # 
    8486######################################### 
     87print 'start mapping' 
     88## parameters of coast coordinates for mapping ## 
     89plt.ion() 
     90x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
     91x_coast = x_ind 
     92y_coast = y_ind 
     93z_coast = z_ind 
     94colormap = cm.s3pcpn_l_r 
     95## start mapping ratio ## 
    8596for imo in range (0, M): 
    86     map_ffgrid.draw_map_l (lon, lat, 0., 10., 0.1, monthly_ratio[imo, :, :], cm.s3pcpn_l_r, 'emissivity ratio [(LAMB - SPEC)/SPEC]*100') 
    87     title(month[imo] + ' 2009 - AMSUB') 
    88     plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/space_evolution/EMIS/emis_rate_vs_ice_type/emis_rate_spec-lamb_AMSUB_' + month[imo] + '2009.png') 
    89     #map_ffgrid.draw_map_l (lon, lat, 0, 5, 1, type_map[imo, :, :], colors.ListedColormap(['1.', '0.25', '0.50', '0.75']), 'Ice Type') 
    90     #title(month[imo] + ' 2009 - OSISAF') 
    91     #plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/space_evolution/EMIS/emis_rate_vs_ice_type/ice_type_OSISAF_01' + month[imo] + '2009.png') 
     97    print 'map month ' +str(month[imo]) 
     98    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, monthly_ratio[imo, :, :], 0., 10., 0.1, colormap, 'emissivity ratio [(LAMB - SPEC)/SPEC]*100') 
     99    title(month[imo] + ' 2009 - AMSUA 89GHz') 
     100    plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_ratio/emis_ratio_spec-lamb_AMSUA89_' + month[imo] + '2009.png') 
     101 
     102 
     103 
    92104 
    93105 
     
    96108# stack monthly means of emissivity ratio in netcdf files # 
    97109########################################################### 
     110print 'start stacking' 
     111## AMSUA ratio ## 
    98112for imo in range (0, M): 
    99     print 'file for ' + month[imo] 
    100     rootgrp = Dataset('/mma/hermozol/Documents/Data/monthly_GLACE/gridded_data/grid_monthly_lamb-spec_ratio_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
    101     rootgrp.createDimension('longitude', len(lon)) 
    102     rootgrp.createDimension('latitude', len(lat)) 
     113    print 'stack in file month ' + str(month[imo]) 
     114    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
     115    rootgrp.createDimension('longitude', len(xdist)) 
     116    rootgrp.createDimension('latitude', len(ydist)) 
    103117    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
    104118    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
    105119    nc_ratio = rootgrp.createVariable('emis_ratio', 'f', ('latitude', 'longitude')) 
    106     nc_lon[:] = lon 
    107     nc_lat[:] = lat 
     120    nc_lon[:] = xdist 
     121    nc_lat[:] = ydist 
    108122    nc_ratio[:] = monthly_ratio[imo, :, :] 
    109123    rootgrp.close() 
    110  
    111  
    112  
    113  
    114  
    115  
    116  
    117  
    118  
    119  
    120  
    121  
    122  
    123 ############################# 
    124 # delimitation ice - no_ice # 
    125 ############################# 
    126 ice_zone = np.zeros([M, ny, nx], float) 
    127 for imo in range (0, M): 
    128     print 'month: ' + month[imo] 
    129     for ilat in range (0, len(lat)): 
    130         for ilon in range (0, len(lon)): 
    131             if (isnan(monthly_ratio[imo, ilat, ilon]) == True): 
    132                 ice_zone[imo, ilat, ilon] = NaN  
    133             else: 
    134                 if (monthly_ratio[imo, ilat, ilon] >= 3): 
    135                     ice_zone[imo, ilat, ilon] = 0. 
    136                 else: 
    137                     ice_zone[imo, ilat, ilon] = 1. 
    138 #    map_ffgrid.draw_map_l (lon, lat, 0., 1.5, 0.5, ice_zone[imo, :, :], colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : ratio < 3.5 %)') 
    139 #    title(month[imo] + ' 2009 - AMSUB') 
    140 #    plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/space_evolution/EMIS/AMSUB_emis_class/AMSUB_emis_classif_ratio-inf-3_' + month[imo] + '2009.png') 
    141  
    142  
    143  
    144 ########################################## 
    145 # calculation of number of pixels of ice # 
    146 ########################################## 
    147 Rt = 6372. 
    148  
    149 dist_along_lon = np.zeros([len(lat), len(lon)], float) 
    150 for ilat in range (0, len(lat) - 1): 
    151     for ilon in range (0, len(lon)): 
    152         # calculate distances along each longitude slice 
    153         arc_along_lon = distance_lon_lat.distance_on_unit_sphere(lat[ilat], lon[ilon], lat[ilat + 1], lon[ilon]) 
    154         dist_along_lon[ilat, ilon] = arc_along_lon * Rt 
    155  
    156 dist_along_lat = np.zeros([len(lat), len(lon)], float) 
    157 for ilon in range (0, len(lon) - 1): 
    158     for ilat in range (0, len(lat)): 
    159         # calculate distances along each latitude slice 
    160         arc_along_lat = distance_lon_lat.distance_on_unit_sphere(lat[ilat], lon[ilon], lat[ilat], lon[ilon + 1]) 
    161         dist_along_lat[ilat, ilon] = arc_along_lat * Rt 
    162  
    163 dist_along_lon = 55. 
    164 pix_area = np.zeros([M, len(lat), len(lon)], float) 
    165 total_ice_area = np.zeros([M], float) 
    166 for imo in range (0, M): 
    167     for ilon in range (0, len(lon)): 
    168         for ilat in range (0, len(lat)): 
    169             if (ice_zone[imo, ilat, ilon] == 0.): 
    170                 pix_area[imo, ilat, ilon] = NaN 
    171             else: 
    172                 if (isnan(ice_zone[imo, ilat, ilon]) == True): 
    173                     pix_area[imo, ilat, ilon] = NaN 
    174                 else: 
    175                     pix_area[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon] 
    176     total_ice_area[imo] = np.sum(reshape(pix_area[imo, :, :], size(pix_area[imo, :, :]))[nonzero(isnan(reshape(pix_area[imo, :, :], size(pix_area[imo, :, :]))) == False)]) 
    177  
    178  
    179 figure() 
    180 plot(total_ice_area, '.b') 
    181 map_ffgrid.draw_map_l (lon, lat, 0., 5300., 10., pix_area, cm.s3pcpn_l_r, 'dist along lon') 
    182  
    183  
    184  
    185  
    186  
    187  
    188  
    189  
    190  
    191  
    192 ######## tests ########### 
    193 np.transpose(dist_along_lon) 
    194  
    195 ice_pixels = np.zeros([M], float) 
    196 for imo in range (0, M): 
    197     ice_vec = reshape(ice_zone[imo, :, :], size(ice_zone[imo, :, :])) 
    198     ice_pixels[imo] = len(np.where(ice_vec == 1.)[0]) 
    199  
    200  
    201  
    202 radius = np.zeros([len(lat)], float) 
    203 dx_km = np.zeros([len(lat)], float) 
    204 for ilat in range (0, len(lat)): 
    205     radius[ilat] = Rt * cos(lat[ilat] * c) # calculate radius of earth at a given latitude 
    206     dx_km[ilat] = c * dx * radius[ilat] # calculate distance in km of longitude step dx (dx = 1.) 
    207  
    208  
    209  
    210  
    211  
    212  
    213  
    214  
    215  
    216  
  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/ice_delimit_AMSUB_data.py

    r40 r44  
    88from mpl_toolkits.basemap import shiftgrid, cm 
    99from netCDF4 import Dataset 
    10 #import arctic_map # function to regrid coast limits 
    11 #import cartesian_grid_test # function to convert grid from polar to cartesian 
     10import arctic_map # function to regrid coast limits 
     11import cartesian_grid_test # function to convert grid from polar to cartesian 
    1212import scipy.special 
    1313import ffgrid2 
    1414import map_ffgrid 
    1515from matplotlib import colors 
    16 import distance_lon_lat 
    17  
     16from matplotlib.font_manager import FontProperties 
     17import map_cartesian_grid 
    1818 
    1919MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) 
     
    2323 
    2424 
    25  
    26 ################### 
    27 # grid parameters # 
    28 ################### 
    29 dx=1. 
    30 dy=0.5 
    31 x0, x1 = -180., 180. 
    32 y0, y1 = 30., 90. 
    33 xvec = np.arange(x0, x1 + dx, dx) 
    34 yvec = np.arange(y0, y1 + dy, dy) 
    35 nx = len(xvec) 
     25######################## 
     26# grid characteristics # 
     27######################## 
     28x0 = -3000. # min limit of grid 
     29x1 = 2500. # max limit of grid 
     30dx = 40. 
     31xvec = np.arange(x0, x1+dx, dx) 
     32nx = len(xvec)  
     33y0 = -3000. # min limit of grid 
     34y1 = 3000. # max limit of grid 
     35dy = 40. 
     36yvec = np.arange(y0, y1+dy, dy) 
    3637ny = len(yvec) 
    3738 
     
    4142# read NETCDF files # 
    4243##################### 
    43 a = 361 * 51 # size of the grid with lats > 65 deg (lon * lat) 
    44 SP = np.zeros([M, 51, 361], float) 
    45 LB = np.zeros([M, 51, 361], float) 
    46 ER = np.zeros([M, 51, 361], float) 
    47 r_vec = np.zeros([M, a], float) 
    48 s_vec = np.zeros([M, a], float) 
    49 l_vec = np.zeros([M, a], float) 
    50 sl_vec = np.zeros([M, a], float) 
    51 for imo in range (0, M): 
    52     print 'file for ' + month[imo] 
    53     ## read netcdf file of emissivity ratio data ## 
    54     print 'ratio file' 
    55     ratio_file = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/grid_monthly_lamb-spec_ratio_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    56     lon_ratio = ratio_file.variables['longitude'][:] 
    57     lat_ratio = ratio_file.variables['latitude'][:] 
    58     emis_ratio = ratio_file.variables['emis_ratio'][:] 
    59     ratio_file.close() 
    60     ## read netcdf file of emissivity LAMB and SPEC data ## 
    61     print 'SPEC LAMB file' 
    62     spec_lamb_file = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/grid_monthly_data_lamb_spec_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    63     lon_sl = spec_lamb_file.variables['longitude'][:] 
    64     lat_sl = spec_lamb_file.variables['latitude'][:] 
    65     e_spec = spec_lamb_file.variables['e_spec'][:] 
    66     e_lamb = spec_lamb_file.variables['e_lamb'][:] 
    67     e_diff = spec_lamb_file.variables['e_spec_lamb'][:] 
    68     spec_lamb_file.close() 
    69     ## keep same dimension between ration and emissivity ## 
    70     dim_equal = np.where(lat_ratio >= 65.)[0] 
    71     ldim = len(dim_equal) 
    72     new_lat = lat_ratio[dim_equal[0] : dim_equal[ldim-1] + 1] # only keep latitude > 65 deg 
    73     ratio = emis_ratio[dim_equal[0] : dim_equal[ldim-1] + 1, :]# only keep corresponding emissivity ratio 
    74     SP[imo, :, :] = e_spec 
    75     LB[imo, :, :] = e_lamb 
    76     ER[imo, :, :] = ratio 
    77     r_vec[imo, :] = reshape(ratio, size(ratio)) 
    78     s_vec[imo, :] = reshape(e_spec, size(e_spec)) 
    79     l_vec[imo, :] = reshape(e_lamb, size(e_lamb)) 
    80     sl_vec[imo, :] = reshape(e_diff, size(e_diff)) 
    81  
    82  
    83  
    84  
    85 nan_index = np.zeros([M], float) 
    86 for imo in range (0, M): 
    87     nan_index[imo] = np.where(isnan(SP[imo, :, 180]) == True)[0][0] - 1 # select index of latitude vector from which data = NaN because of selection of near nadir angles only // above this index (excluded), all NaN values 
    88  
    89  
    90  
    91 # map test 
    92 map_ffgrid.draw_map_l (lon_sl, lat_sl, 0.6, 1.02, 0.01, SP[0, :, 180], cm.s3pcpn_l_r, 'emissivity SPEC') 
    93  
    94  
    95  
    96 ############################################################ 
    97 # distribution of emissivity rate emissivity spec and lamb # 
    98 ############################################################ 
    99 ## histogram plots 
    100  
    101 ##### emissivity rate #### 
     44emis_spec_month = np.zeros([M, ny, nx], float) 
     45emis_lamb_month = np.zeros([M, ny, nx], float) 
     46ratio_emis = np.zeros([M, ny, nx], float) 
     47ratio_emis_vec = np.zeros([M, ny * nx], float) 
     48spec_emis_vec = np.zeros([M, ny * nx], float) 
     49lamb_emis_vec = np.zeros([M, ny * nx], float) 
     50for imo in range (0, M): 
     51    print 'month: ' + month[imo] 
     52    ### emissivity ### 
     53    print 'open emissivity file' 
     54    fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
     55    xdist = fichier_emis.variables['longitude'][:] 
     56    ydist = fichier_emis.variables['latitude'][:] 
     57    day = fichier_emis.variables['days'][:] 
     58    emis_spec = fichier_emis.variables['e_spec'][:] 
     59    emis_lamb = fichier_emis.variables['e_lamb'][:] 
     60    fichier_emis.close() 
     61    ##### calculate monthly mean of emis ##### 
     62    for ilat in range(0, len(ydist)): 
     63        for ilon in range (0, len(xdist)): 
     64            emis_spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_spec[ilat, ilon, 0 : month_day[imo]]) == False)]) 
     65            emis_lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_lamb[ilat, ilon, 0 : month_day[imo]]) == False)]) 
     66    ### monthly emissivity ratio ### 
     67    print 'open emissivity ratio file' 
     68    fichier_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA23_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
     69    xdist = fichier_ratio.variables['longitude'][:] 
     70    ydist = fichier_ratio.variables['latitude'][:] 
     71    ratio_emis[imo, :, :] = fichier_ratio.variables['emis_ratio'][:] 
     72    fichier_ratio.close() 
     73    ### reshape in monthly vectors all data in 2D arrays ### 
     74    print 'reshape data' 
     75    ratio_emis_vec[imo, :] = reshape(ratio_emis[imo, :, :], size(ratio_emis[imo, :, :])) 
     76    spec_emis_vec[imo, :] = reshape(emis_spec_month[imo, :, :], size(emis_spec_month[imo, :, :])) 
     77    lamb_emis_vec[imo, :] = reshape(emis_lamb_month[imo, :, :], size(emis_lamb_month[imo, :, :])) 
     78 
     79 
     80''' 
     81# test: 
     82x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
     83x_coast = x_ind 
     84y_coast = y_ind 
     85z_coast = z_ind 
     86map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, emis_spec[:, :, 0], 0.4, 1.2, 0.01, cm.s3pcpn_l_r, 'test') 
     87''' 
     88 
     89 
    10290c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) 
     91################################### 
     92# distribution of SPEC emissivity # 
     93################################### 
     94hist_vals_spec = np.zeros([M, 50], float) 
     95corresp_emis_spec = np.zeros([M, 50], float) 
     96for imo in range (0, M): 
     97    hist_vals_spec[imo, :] = hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 
     98    for ibin in range (0, 50): 
     99        corresp_emis_spec[imo, ibin] = mean(hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 
     100 
     101## plot first six months of spec emissivity histograms ## 
     102ion() 
    103103figure() 
    104104for imo in range (0, 6): 
    105     hist(r_vec[imo, :][nonzero(isnan(r_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step', color = str(c[imo]), label = str(month[imo]))[0] 
    106  
    107 legend() 
    108 grid() 
    109 xlim(0, 10) 
    110 xticks(np.arange(0., 12., 1.), np.array(['0%', '1%', '2%', '3%', '4%', '5%', '6%', '7%', '8%', '9%', '10%'])) 
     105    plot(corresp_emis_spec[imo], hist_vals_spec[imo], c = str(c[imo]), label = str(month[imo])) 
     106 
     107grid() 
     108xlim(corresp_emis_spec.min() - 0.02, 1.02) 
     109xlabel('emissivity SPEC') 
     110ylabel('frequency of occurence') 
     111fontP = FontProperties() 
     112fontP.set_size('small') 
     113legend(loc = 2, prop = fontP) 
     114plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA23_JANUARY-JUNE_2009.png') 
     115## plot six following months of spec emissivity histograms ## 
     116figure() 
     117for imo in range (6, M): 
     118    plot(corresp_emis_spec[imo], hist_vals_spec[imo], c = str(c[imo - 6]), label = str(month[imo])) 
     119 
     120grid() 
     121xlim(corresp_emis_spec.min() - 0.02, 1.02) 
     122xlabel('emissivity SPEC') 
     123ylabel('frequency of occurence') 
     124legend(loc = 9, prop = fontP) 
     125plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA23_JULY-DECEMBER_2009.png') 
     126 
     127# find the emissivity value corresponding to the threshold of ice/no_ice limit 
     128emis_lim_spec = np.zeros([M], float) 
     129for imo in range (0, M): 
     130    print 'month ' + str(month[imo]) 
     131    bb = np.where((corresp_emis_spec[imo, :] < 0.6))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 
     132    aa = np.where(hist_vals_spec[imo, bb] < max(hist_vals_spec[imo, bb]) / 2.)[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity 
     133    cc = np.where(hist_vals_spec[imo, bb] == max(hist_vals_spec[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 
     134    dd = np.where(aa > cc[0])[0] 
     135    emis_lim_spec[imo] = corresp_emis_spec[imo, bb][aa][dd[0]] 
     136 
     137# plot monthly evolution of emissivity spec threshold used for delimitation 
     138figure() 
     139plot(np.arange(0, M, 1), emis_lim_spec, 'b-+', label = 'emis SPEC') 
     140ylabel('emissivity SPEC threshold - AMSUA 23GHz') 
     141grid() 
     142xticks(np.arange(0, M, 1), month, rotation = 25) 
     143legend(loc = 2) 
     144xlim(-1, M) 
     145plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_SPEC_AMSUA23_2009.png') 
     146 
     147 
     148################################### 
     149# distribution of LAMB emissivity # 
     150################################### 
     151hist_vals_lamb = np.zeros([M, 50], float) 
     152corresp_emis_lamb = np.zeros([M, 50], float) 
     153for imo in range (0, M): 
     154    hist_vals_lamb[imo, :] = hist(lamb_emis_vec[imo, :][nonzero(isnan(lamb_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 
     155    for ibin in range (0, 50): 
     156        corresp_emis_lamb[imo, ibin] = mean(hist(lamb_emis_vec[imo, :][nonzero(isnan(lamb_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 
     157 
     158## plot first six months of spec emissivity histograms ## 
     159figure() 
     160for imo in range (0, 6): 
     161    plot(corresp_emis_lamb[imo], hist_vals_lamb[imo], c = str(c[imo]), label = str(month[imo])) 
     162 
     163grid() 
     164xlim(corresp_emis_lamb.min() - 0.02, 1.02) 
     165xlabel('emissivity LAMB') 
     166ylabel('frequency of occurence') 
     167legend(loc = 2, prop = fontP) 
     168plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA23_JANUARY-JUNE_2009.png') 
     169## plot six following months of spec emissivity histograms ## 
     170figure() 
     171for imo in range (6, M): 
     172    plot(corresp_emis_lamb[imo], hist_vals_lamb[imo], c = str(c[imo - 6]), label = str(month[imo])) 
     173 
     174grid() 
     175xlim(corresp_emis_lamb.min() - 0.02, 1.02) 
     176xlabel('emissivity LAMB') 
     177ylabel('frequency of occurence') 
     178legend(loc = 9, prop = fontP) 
     179plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA23_JULY-DECEMBER_2009.png') 
     180 
     181# find the emissivity value corresponding to the threshold of ice/no_ice limit 
     182emis_lim_lamb = np.zeros([M], float) 
     183for imo in range (0, M): 
     184    print 'mont ' + str(month[imo]) 
     185    bb = np.where((corresp_emis_lamb[imo, :] < 0.6))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 
     186    aa = np.where(hist_vals_lamb[imo, bb] < max(hist_vals_lamb[imo, bb]) / 2.)[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity 
     187    cc = np.where(hist_vals_lamb[imo, bb] == max(hist_vals_lamb[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 
     188    dd = np.where(aa > cc[0])[0] 
     189    emis_lim_lamb[imo] = corresp_emis_lamb[imo, bb][aa][dd[0]] 
     190 
     191# plot monthly evolution of emissivity spec threshold used for delimitation 
     192figure() 
     193plot(np.arange(0, M, 1), emis_lim_lamb, 'b-+', label = 'emis LAMB') 
     194ylabel('emissivity LAMB threshold - AMSUB 23GHz') 
     195grid() 
     196xticks(np.arange(0, M, 1), month, rotation = 25) 
     197legend(loc = 2) 
     198xlim(-1, M) 
     199plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_LAMB_AMSUA23_2009.png') 
     200 
     201 
     202 
     203################################### 
     204# distribution of emissivity rate # 
     205################################### 
     206hist_vals_rate = np.zeros([M, 50], float) 
     207corresp_emis_rate = np.zeros([M, 50], float) 
     208for imo in range (0, M): 
     209    hist_vals_rate[imo, :] = hist(ratio_emis_vec[imo, :][nonzero(isnan(ratio_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 
     210    for ibin in range (0, 50): 
     211        corresp_emis_rate[imo, ibin] = mean(hist(ratio_emis_vec[imo, :][nonzero(isnan(ratio_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 
     212 
     213## plot first six months of spec emissivity histograms ## 
     214figure() 
     215for imo in range (0, 6): 
     216    plot(corresp_emis_rate[imo], hist_vals_rate[imo], c = str(c[imo]), label = str(month[imo])) 
     217 
     218grid() 
     219xlim(corresp_emis_rate.min() - 0.02, corresp_emis_rate.max() + 0.02) 
    111220xlabel('emissivity rate (%)') 
    112 ylabel('frequency of appearence') 
    113 plt.savefig('/mma/hermozol/Documents/figure_output/ice_class_AMSUB/emiss_rate_AMSUB_JANUARY-JUNE_2009.png') 
    114  
     221ylabel('frequency of occurence') 
     222legend(prop = fontP) 
     223plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA23_JANUARY-JUNE_2009.png') 
     224## plot six following months of spec emissivity histograms ## 
    115225figure() 
    116226for imo in range (6, M): 
    117     hist(r_vec[imo, :][nonzero(isnan(r_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step', color = str(c[imo - 6]), label = str(month[imo])) 
    118  
    119 legend() 
    120 grid() 
    121 xlim(0, 10) 
    122 xticks(np.arange(0., 12., 1.), np.array(['0%', '1%', '2%', '3%', '4%', '5%', '6%', '7%', '8%', '9%', '10%'])) 
     227    plot(corresp_emis_rate[imo], hist_vals_rate[imo], c = str(c[imo - 6]), label = str(month[imo])) 
     228 
     229grid() 
     230xlim(corresp_emis_rate.min() - 0.02, corresp_emis_rate.max() + 0.02) 
    123231xlabel('emissivity rate (%)') 
    124 ylabel('frequency of appearence') 
    125 plt.savefig('/mma/hermozol/Documents/figure_output/ice_class_AMSUB/emiss_rate_AMSUB_JULY-DECEMBER_2009.png') 
    126  
    127 ##### SPEC emissivity #### 
    128 hist_vals = np.zeros([M, 50], float) 
    129 corresp_emis = np.zeros([M, 50], float) 
    130 for imo in range (0, M): 
    131     hist_vals[imo, :] = hist(s_vec[imo, :][nonzero(isnan(s_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 
    132     for ibin in range (0, 50): 
    133         corresp_emis[imo, ibin] = mean(hist(s_vec[imo, :][nonzero(isnan(s_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    134  
    135 figure() 
    136 for imo in range (0, 6): 
    137     plot(corresp_emis[imo], hist_vals[imo], c = str(c[imo]), label = str(month[imo])) 
    138  
    139 grid() 
    140 xlim(0.5, 1.02) 
    141 xlabel('emissivity SPEC') 
    142 ylabel('frequency of occurence') 
     232ylabel('frequency of occurence') 
     233legend(loc = 9, prop = fontP) 
     234plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA23_JULY-DECEMBER_2009.png') 
     235 
     236# no definition of threshold for ICE / NO_ICE delimitation because no apparent distinction signal // use of this signal for internal ice classification 
     237''' 
     238# find the emissivity value corresponding to the threshold of ice/no_ice limit 
     239emis_lim_rate = np.zeros([M], float) 
     240for imo in range (0, M): 
     241    bb = np.where((corresp_emis_rate[imo, :] < 0.6))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 
     242    aa = np.where(hist_vals_rate[imo, bb] < max(hist_vals_rate[imo, bb]) / 4.)[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity 
     243    cc = np.where(hist_vals_rate[imo, bb] == max(hist_vals_rate[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 
     244    dd = np.where(aa > cc[0])[0] 
     245    emis_lim_rate[imo] = corresp_emis_rate[imo, bb][aa][dd[0]] 
     246 
     247# plot monthly evolution of emissivity spec threshold used for delimitation 
     248figure() 
     249plot(np.arange(0, M, 1), emis_lim_rate, 'b-+', label = 'emis LAMB') 
     250ylabel('emissivity rate threshold - AMSUB 89GHz') 
     251grid() 
     252xticks(np.arange(0, M, 1), month, rotation = 25) 
    143253legend(loc = 2) 
    144 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_spec_AMSUB_JANUARY-JUNE_2009.png') 
    145  
    146  
    147 figure() 
    148 for imo in range (6, M): 
    149     plot(corresp_emis[imo], hist_vals[imo], c = str(c[imo - 6]), label = str(month[imo])) 
    150  
    151 grid() 
    152 xlim(0.5, 1.02) 
    153 xlabel('emissivity SPEC') 
    154 ylabel('frequency of occurence') 
    155 legend(loc = 2) 
    156 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_spec_AMSUB_JULY-DECEMBER_2009.png') 
    157  
    158  
    159 # find the emissivity value corresponding to the threshold of ice/no_ice limit 
    160 # after many tests, we choose an emissivity SPEC value which frequency of occurence is minimal after the peak of frequency around 0.7  
    161 emis_lim = np.zeros([M], float) 
    162 for imo in range (0, M): 
    163     bb = np.where((corresp_emis[imo, :] > 0.7) & (corresp_emis[imo, :] < 0.8))[0] 
    164     aa = np.where(hist_vals[imo, bb] == min(hist_vals[imo, bb]))[0] 
    165     emis_lim[imo] = corresp_emis[imo, bb][aa] 
    166  
    167 figure() 
    168 plot(np.arange(0, M, 1), emis_lim, 'b-+', label = 'emis SPEC') 
    169 ylabel('emissivity SPEC threshold - AMSUB 89GHz') 
    170 grid() 
    171 xticks(np.arange(0, M, 1), month, rotation = 25) 
    172 legend(loc = 2) 
    173 xlim(-1, M) 
    174 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emis_lim_frequ_SPEC_AMSUB_2009.png') 
    175  
    176  
    177 ##### LAMB emissivity #### 
    178 hist_vals_l = np.zeros([M, 50], float) 
    179 corresp_emis_l = np.zeros([M, 50], float) 
    180 for imo in range (0, M): 
    181     hist_vals_l[imo, :] = hist(l_vec[imo, :][nonzero(isnan(l_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 
    182     for ibin in range (0, 50): 
    183         corresp_emis_l[imo, ibin] = mean(hist(l_vec[imo, :][nonzero(isnan(l_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    184  
    185 figure() 
    186 for imo in range (0, 6): 
    187     plot(corresp_emis_l[imo], hist_vals_l[imo], c = str(c[imo]), label = str(month[imo])) 
    188  
    189 grid() 
    190 xlim(0.5, 1.02) 
    191 xlabel('emissivity LAMB') 
    192 ylabel('frequency of occurence') 
    193 legend(loc = 2) 
    194 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_lamb_AMSUB_JANUARY-JUNE_2009.png') 
    195  
    196  
    197 figure() 
    198 for imo in range (6, M): 
    199     plot(corresp_emis_l[imo], hist_vals_l[imo], c = str(c[imo - 6]), label = str(month[imo])) 
    200  
    201 grid() 
    202 xlim(0.5, 1.02) 
    203 xlabel('emissivity LAMB') 
    204 ylabel('frequency of occurence') 
    205 legend(loc = 2) 
    206 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_lamb_AMSUB_JULY-DECEMBER_2009.png') 
    207  
    208  
    209 # find the emissivity value corresponding to the threshold of ice/no_ice limit 
    210 # after many tests, we choose an emissivity LAMB value which frequency of occurence is minimal after the peak of frequency around 0.73  
    211 emis_lim_l = np.zeros([M], float) 
    212 for imo in range (0, M): 
    213     bb = np.where((corresp_emis_l[imo, :] > 0.72) & (corresp_emis_l[imo, :] < 0.8))[0] 
    214     aa = np.where(hist_vals_l[imo, bb] == min(hist_vals_l[imo, bb]))[0] 
    215     emis_lim_l[imo] = corresp_emis_l[imo, bb][aa] 
    216  
    217 figure() 
    218 plot(np.arange(0, M, 1), emis_lim_l, 'b-+', label = 'emis LAMB') 
    219 ylabel('emissivity LAMB threshold - AMSUB 89GHz') 
    220 grid() 
    221 xticks(np.arange(0, M, 1), month, rotation = 25) 
    222 legend(loc = 2) 
    223 xlim(-1, M) 
    224 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emis_lim_frequ_LAMB_AMSUB_2009.png') 
    225  
    226  
    227  
    228 ####################################################### 
    229 # delimitation ice - no_ice with different conditions # 
    230 ####################################################### 
    231 ## condition : threshold of SPEC emissivity ('SP') = f[month] ==> take out all emissivities <= f[month] (open sea) // f[month] depends on the month, around 0.7 
    232 ice_zone_l = np.zeros([M, 51, 361], float) 
    233 for imo in range (0, M): 
    234     print 'month: ' + month[imo] 
    235     for ilat in range (0, 51): 
    236         for ilon in range (0, 361): 
    237             if (isnan(LB[imo, ilat, ilon]) == True): 
    238                 ice_zone_l[imo, ilat, ilon] = NaN  
     254xlim(-1, M) 
     255plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_rate_AMSUA89_2009.png') 
     256''' 
     257 
     258 
     259############################################################################## 
     260# plot comparison of emissivity thresholds between spec and lamb assumptions # 
     261############################################################################## 
     262figure() 
     263plot(np.arange(0, M, 1), emis_lim_spec, 'b-+', label = 'emis SPEC') 
     264plot(np.arange(0, M, 1), emis_lim_lamb, 'g-+', label = 'emis LAMB') 
     265ylabel('emissivity threshold - AMSUA 23GHz') 
     266grid() 
     267xticks(np.arange(0, M, 1), month, rotation = 25) 
     268legend(loc = 2, prop = fontP) 
     269xlim(-1, M) 
     270plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_spec_and_lamb_AMSUA23_2009.png') 
     271 
     272 
     273 
     274 
     275################################################################ 
     276# delimitation ice - no_ice with emissivity or ratio threshold # 
     277################################################################ 
     278x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
     279x_coast = x_ind 
     280y_coast = y_ind 
     281z_coast = z_ind 
     282ice_zone = np.zeros([M, 151, 139], float) 
     283for imo in range (0, M): 
     284    print 'month: ' + str(month[imo]) 
     285    for ilat in range (0, 151): 
     286        for ilon in range (0, 139): 
     287            if (isnan(emis_spec_month[imo, ilat, ilon]) == True): 
     288                ice_zone[imo, ilat, ilon] = NaN  
    239289            else: 
    240                 if (LB[imo, ilat, ilon] <= emis_lim_l[imo]): # use the monthly SPEC emissivity threshold  
    241                     ice_zone_l[imo, ilat, ilon] = 0. 
     290                if (emis_spec_month[imo, ilat, ilon] <= emis_lim_spec[imo]): # use the monthly SPEC emissivity threshold  
     291                    ice_zone[imo, ilat, ilon] = 0. 
    242292                else: 
    243                     ice_zone_l[imo, ilat, ilon] = 1.          
    244     map_ffgrid.draw_map_l (lon_sl, lat_sl, -1., 2., 1., ice_zone_l[imo, :, :], colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_LAMB > ' + str(emis_lim_l[imo]) + ')') 
    245     title(month[imo] + ' 2009 - AMSUB') 
    246     plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/AMSUB_emis_class/AMSUB_emis_classif_elamb_thresh' + '_' + month[imo] + '2009.png') 
    247  
    248  
    249  
    250  
    251  
    252 ########################################## 
    253 # calculation of number of pixels of ice # 
    254 ########################################## 
    255 Rt = 6372. 
    256  
    257 '''dist_along_lon = np.zeros([len(lat_sl), len(lon_sl)], float) 
    258 for ilat in range (0, len(lat_sl) - 1): # leave last lat index (lat = 90.) 
    259     for ilon in range (0, len(lon_sl)): 
    260         # calculate distances along each longitude slice 
    261         arc_along_lon = distance_lon_lat.distance_on_unit_sphere(lat_sl[ilat], lon_sl[ilon], lat_sl[ilat + 1], lon_sl[ilon]) 
    262         dist_along_lon[ilat, ilon] = arc_along_lon * Rt''' 
    263  
    264 dist_along_lon =  55.60618997 
    265 dist_along_lat = np.zeros([len(lat_sl), len(lon_sl)], float) 
    266 for ilon in range (0, len(lon_sl) - 1): # leave last lon index (lon = 180.) 
    267     for ilat in range (0, len(lat_sl)): 
    268         # calculate distances along each latitude slice 
    269         arc_along_lat = distance_lon_lat.distance_on_unit_sphere(lat_sl[ilat], lon_sl[ilon], lat_sl[ilat], lon_sl[ilon + 1]) 
    270         dist_along_lat[ilat, ilon] = arc_along_lat * Rt 
    271  
    272 # deal with last lon index (lon = 180.) 
    273 for ilat in range (0, len(lat_sl)): 
    274     dist_along_lat[ilat, 360] = dist_along_lat[ilat, 359] 
    275  
    276  
    277 ## calculate the area of each pixel AMSUB ## 
    278 pix_area_l = np.zeros([M, len(lat_sl), len(lon_sl)], float) 
    279 total_ice_area_l = np.zeros([M], float) 
    280 for imo in range (0, M): 
    281     for ilon in range (0, len(lon_sl)): 
    282         for ilat in range (0, len(lat_sl)): 
    283             if (ice_zone_l[imo, ilat, ilon] == 0.): 
    284                 pix_area_l[imo, ilat, ilon] = NaN 
    285             else: 
    286                 if ((lat_sl[ilat] <= 83.5) & (isnan(ice_zone_l[imo, ilat, ilon]) == True)): 
    287                     pix_area_l[imo, ilat, ilon] = NaN 
    288                 else: 
    289                     if ((lat_sl[ilat] > 83.5)): 
    290                         pix_area_l[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon] 
    291                     else: 
    292                         pix_area_l[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon] 
    293     total_ice_area_l[imo] = np.sum(reshape(pix_area_l[imo, :, :], size(pix_area_l[imo, :, :]))[nonzero(isnan(reshape(pix_area_l[imo, :, :], size(pix_area_l[imo, :, :]))) == False)]) 
    294     map_ffgrid.draw_map_l (lon_sl, lat_sl, 0., 2614., 10., pix_area_l[imo,:,:], cm.s3pcpn_l_r, 'pixel area') 
    295     title(str(month[imo]) + '2009 - AMSUB') 
    296     plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_pixel_area_AMSUB_' + month[imo] + '2009.png') 
    297  
    298  
    299  
    300  
    301  
    302  
    303 ################################################## 
    304 # Read OSISAF ice types to have total ice extent # 
    305 ################################################## 
    306 osi_type = np.zeros([M, 31, 51, 361], float) 
    307 for imo in range (0, M): 
    308     print 'month ' + month[imo] 
    309     # osisaf 
    310     print 'read OSISAF' 
    311     fichier_osi = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_ffgrid/OSISAF_ice_types_ffgrid_' + month[imo] + '_2009.nc', 'r', format='NETCDF3_CLASSIC') 
    312     lon_osi = fichier_osi.variables['longitude'][:] 
    313     lat_osi = fichier_osi.variables['latitude'][:] 
    314     days_osi = fichier_osi.variables['days'][:] 
    315     osi_ice_type = fichier_osi.variables['osi_ice_type'][:] 
    316     osi_type[imo, 0 : month_day[imo], :, :] = osi_ice_type[:, np.where(lat_osi == 65.)[0][0] : len(lat_osi),  :] # keep the same dimension in lat (start at 65 deg) 
    317     fichier_osi.close() 
    318  
    319  
    320 osi_ice = np.zeros([M, 31, 51, 361], float) 
    321 for imo in range (0, M): 
    322     for ilon in range (0, len(lon_sl)): 
    323         for ilat in range (0, len(lat_sl)): 
     293                    ice_zone[imo, ilat, ilon] = 1.    
     294    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone[imo, :, :], -1, 2., 1., colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_SPEC > ' + str(emis_lim_spec[imo])[0:6] + ')') 
     295    title(str(month[imo]) + ' 2009 - AMSUA 23GHz') 
     296    plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/AMSUA23_ice_emis_spec_thresh' + '_' + month[imo] + '2009.png') 
     297 
     298 
     299 
     300 
     301 
     302########################### 
     303# calculation of ice area # 
     304########################### 
     305# nb of pixels near pole = 926 
     306pix_area = 40. * 40. 
     307nb_pix = np.zeros([M], float) 
     308total_ice_area = np.zeros([M], float) 
     309for imo in range (0, M): 
     310    ice = reshape(ice_zone[imo, :, :], size(ice_zone[imo, :, :])) 
     311    nb_pix[imo] = len(np.where(ice == 1.)[0]) 
     312    total_ice_area[imo] = (pix_area * nb_pix[imo]) + (926. * pix_area) 
     313 
     314 
     315spec_ice_area = total_ice_area 
     316 
     317 
     318figure() 
     319plot(lamb_ice_area, 'g', label = 'lamb') 
     320plot(spec_ice_area, 'b', label = 'spec') 
     321plot(total_monthly_ice_area_osisaf, 'r', label = 'OSISAF') 
     322ylabel('total ice area (square km)') 
     323grid() 
     324xticks(np.arange(0, M, 1), month, rotation = 25) 
     325legend(loc = 3, prop = fontP) 
     326title('AMSUA 23GHz') 
     327xlim(-1, M) 
     328plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/total_ice_area_AMSUA23_OSISAF_2009.png') 
     329 
     330########################## 
     331# stack of ice area data # 
     332########################## 
     333# stack in .dat file OSISAF data 
     334print 'start stacking in .dat file' 
     335data_osisaf = open ('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_daily_ice_no-ice_2009.dat', 'a') 
     336for imo in range (0, M): 
     337    for ijr in range (0, month_day[imo]): 
     338        data_osisaf.write(('%(months)10s    %(days)2.0f    %(total_monthly_ice_area_osisaf)10.5f    %(total_ice_area_osisaf)10.5f   \n' % { 
     339                                            'months':month[imo], 
     340                                            'days':np.arange(1, month_day[imo] + 1)[ijr], 
     341                                            'total_monthly_ice_area_osisaf':total_monthly_ice_area_osisaf[imo], 
     342                                            'total_ice_area_osisaf':total_ice_area_osisaf[imo, ijr], 
     343                                            })) 
     344 
     345data_osisaf.close() 
     346 
     347 
     348# stack in .dat file  
     349print 'start stacking in .dat file' 
     350data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/2009.dat', 'a') 
     351for imo in range (0, 50): 
     352    for ii in range (0, month_day[imo]): 
     353        data_classif.write(('%(months)10s    %(hist_vals_spec)10.5f    %(corresp_emis_spec)10.5f    %(hist_vals_lamb)10.5f    %(orresp_emis_lamb)10.5f    %(hist_vals_rate)10.5f    %(corresp_emis_rate)10.5f    %(emis_lim_spec)10.5f    (emis_lim_lamb)10.5f    (spec_ice_area)10.5f    (lamb_ice_area)10.5f   \n' % { 
     354                                            'months':month[imo], 
     355                                            'hist_vals_spec':hist_vals_spec[imo, ii], 
     356                                            'corresp_emis_spec':corresp_emis_spec[imo, ii], 
     357                                            'hist_vals_lamb':hist_vals_lamb[imo, ii], 
     358                                            'corresp_emis_lamb':corresp_emis_lamb[imo, ii], 
     359                                            'hist_vals_rate':hist_vals_rate[imo, ii], 
     360                                            'emis_lim_spec':emis_lim_spec[imo], 
     361                                            'emis_lim_lamb':emis_lim_lamb[imo], 
     362                                            'spec_ice_area':spec_ice_area[imo], 
     363                                            'lamb_ice_area':lamb_ice_area[imo], 
     364                                            })) 
     365 
     366data_classif.close() 
     367 
     368# stack in netcdf file 
     369print 'start stacking' 
     370for imo in range (0, M): 
     371    print 'stack in file month ' + str(month[imo]) 
     372    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + month[imo] + '2009_AMSUA23.nc', 'w', format='NETCDF3_CLASSIC') 
     373    rootgrp.createDimension('longitude', len(xdist)) 
     374    rootgrp.createDimension('latitude', len(ydist)) 
     375    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
     376    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
     377    nc_ice = rootgrp.createVariable('ice_area', 'f', ('latitude', 'longitude')) 
     378    nc_lon[:] = xdist 
     379    nc_lat[:] = ydist 
     380    nc_ice[:] = ice_zone[imo, :, :] 
     381    rootgrp.close() 
     382 
     383 
     384############################################################################################################# 
     385 
     386######################################### 
     387# comparison with OSISAF sea ice extent # 
     388######################################### 
     389osi_ice = np.zeros([M, 31, 151, 139]) 
     390for imo in range (0, M): 
     391    print 'open OSISAF file month ' + str(month[imo]) 
     392    fichier_osisaf = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_cartesian_grid/OSISAF_ice_types_cartes_grid_' + str(month[imo]) + '_2009.nc', 'r', format='NETCDF3_CLASSIC') 
     393    xdist = fichier_osisaf.variables['longitude'][:] 
     394    ydist = fichier_osisaf.variables['latitude'][:] 
     395    day = fichier_osisaf.variables['days'][:] 
     396    osi_type = fichier_osisaf.variables['osi_ice_type'][:] 
     397    fichier_osisaf.close() 
     398    for ilon in range (0, len(xdist)): 
     399        for ilat in range (0, len(ydist)): 
    324400            for ijr in range (0, month_day[imo]): 
    325                 if (osi_type[imo, ijr, ilat, ilon] >= 2.): 
     401                if ((isnan(osi_type[ijr, ilat, ilon]) == True) & (xdist[ilon] >= -312.5) & (xdist[ilon] <= 328.5) & (ydist[ilat] >=-334.3) & (ydist[ilat] <= 338.5)): 
    326402                    osi_ice[imo, ijr, ilat, ilon] = 1. 
    327403                else: 
    328                     if ((lat_sl[ilat] <= 86.) & (isnan(osi_type[imo, ijr, ilat, ilon]) == True)): 
    329                         osi_ice[imo, ijr, ilat, ilon] = NaN 
    330                     else:  
    331                         if (lat_sl[ilat] > 86.): 
    332                             osi_ice[imo, ijr, ilat, ilon] = 1. 
    333                         else: 
     404                    if (osi_type[ijr, ilat, ilon] >= 2.): 
     405                        osi_ice[imo, ijr, ilat, ilon] = 1. 
     406                    else: 
     407                        if (isnan(osi_type[ijr, ilat, ilon]) == True): 
    334408                            osi_ice[imo, ijr, ilat, ilon] = NaN 
    335  
    336  
    337  
    338 osi_ice_month = np.zeros([M, 51, 361], float)                        
    339 for imo in range (0, M): 
    340     for ilon in range (0, len(lon_sl)): 
    341         for ilat in range (0, len(lat_sl)): 
    342             how_many_ice_pix = len(np.where(osi_ice[imo, :, ilat, ilon] == 1.)[0]) 
    343             if (how_many_ice_pix >= 15): 
    344                 osi_ice_month[imo, ilat, ilon] = 1. 
    345             else: 
    346                 osi_ice_month[imo, ilat, ilon] = NaN 
    347  
    348  
    349  
    350 ########################################### 
    351 # calculate the area of each pixel OSISAF # 
    352 ########################################### 
    353 pix_area_osi = np.zeros([M, len(lat_sl), len(lon_sl)], float) 
    354 total_ice_area_osi = np.zeros([M], float) 
    355 for imo in range (0, M): 
    356     for ilon in range (0, len(lon_sl)): 
    357         for ilat in range (0, len(lat_sl)): 
    358             if (osi_ice_month[imo, ilat, ilon] == 1.): 
    359                 pix_area_osi[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon] 
    360             else: 
    361                 pix_area_osi[imo, ilat, ilon] = NaN 
    362     total_ice_area_osi[imo] = np.sum(reshape(pix_area_osi[imo, :, :], size(pix_area_osi[imo, :, :]))[nonzero(isnan(reshape(pix_area_osi[imo, :, :], size(pix_area_osi[imo, :, :]))) == False)]) 
    363     map_ffgrid.draw_map_l (lon_sl, lat_sl, 0., 2614., 10., pix_area_osi[imo,:,:], cm.s3pcpn_l_r, 'pixel area') 
    364     title(str(month[imo]) + '2009 - OSISAF') 
    365     plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_pix_area/ice_pixel_area_OSISAF_' + month[imo] + '2009.png') 
    366  
    367  
    368  
    369  
     409                        else:  
     410                            osi_ice[imo, ijr, ilat, ilon] = 0. 
     411 
     412 
     413# test: 
     414colormap = colors.ListedColormap(['cyan','blue']) 
     415map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_ice[10, 0, :, :], -1, 2, 1, colormap, 'ice_type') 
     416 
     417 
     418################################## 
     419# calculation of OSISAF ice area # 
     420################################## 
     421pix_area = 40. * 40. 
     422nb_pix_osisaf = np.zeros([M, 31], float) 
     423total_ice_area_osisaf = np.zeros([M, 31], float) 
     424total_monthly_ice_area_osisaf = np.zeros([M], float) 
     425for imo in range (0, M): 
     426    for ijr in range (0, month_day[imo]): 
     427        ice = reshape(osi_ice[imo, ijr, :, :], size(osi_ice[imo, ijr, :, :])) 
     428        nb_pix_osisaf[imo, ijr] = len(np.where(ice == 1.)[0]) 
     429        total_ice_area_osisaf[imo, ijr] = pix_area * nb_pix_osisaf[imo, ijr] 
     430    total_monthly_ice_area_osisaf[imo] = mean(total_ice_area_osisaf[imo, 0 : month_day[imo]][nonzero(isnan(total_ice_area_osisaf[imo, 0 : month_day[imo]]) == False)]) + (926. * pix_area) 
     431 
     432 
     433 
     434 
     435 
     436################################### 
     437# plot time evolution of ice area # 
     438################################### 
    370439figure() 
    371440plot(total_ice_area, '-+b', label = 'AMSUB SPEC') 
     
    381450 
    382451 
    383 map_ffgrid.draw_map_l (lon, lat, 0., 5300., 10., pix_area, cm.s3pcpn_l_r, 'dist along lon') 
    384  
    385  
    386  
    387  
    388  
     452 
     453 
     454 
     455 
     456 
     457 
  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_cartesian_gridded_AMSUA23.py

    r42 r44  
    5353for imo in range (0, M): 
    5454    print 'open file ' + str(month[imo]) 
    55     file_amsua23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    56     lon_amsua = file_amsua23.variables['longitude'][:] 
    57     lat_amsua = file_amsua23.variables['latitude'][:] 
    58     days = file_amsua23.variables['days'][:] 
    59     es_a = file_amsua23.variables['e_spec'][:] 
    60     el_a = file_amsua23.variables['e_lamb'][:] 
    61     esl_a = file_amsua23.variables['e_spec_lamb'][:] 
    62     file_amsua23.close() 
     55    file_amsua = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
     56    lon_amsua = file_amsua.variables['longitude'][:] 
     57    lat_amsua = file_amsua.variables['latitude'][:] 
     58    days = file_amsua.variables['days'][:] 
     59    es_a = file_amsua.variables['e_spec'][:] 
     60    el_a = file_amsua.variables['e_lamb'][:] 
     61    esl_a = file_amsua.variables['e_spec_lamb'][:] 
     62    file_amsua.close() 
    6363    print 'compute monthly mean of data for map' 
    6464    for ilon in range (0, len(lon_amsua)): 
     
    6969 
    7070 
    71  
     71print 'map of data' 
     72ion() 
     73x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
     74x_coast = x_ind 
     75y_coast = y_ind 
     76z_coast = z_ind 
     77colormap = cm.s3pcpn_l_r 
    7278for imo in range (0, M): 
    73     x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    74     x_coast = x_ind 
    75     y_coast = y_ind 
    76     z_coast = z_ind 
    77     colormap = cm.s3pcpn_l_r 
    7879    # emiss SPEC 
    79     map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, es_a_month[imo, :, :], 0.45, 1.02, 0.01, colormap, 'emissivity SPEC') 
    80     title(month[imo] + ' 2009 - AMSUA 23GHz') 
    81     savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/emis_spec_AMSUA23_' + month[imo] + '2009.png') 
     80    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, es_a_month[imo, :, :], 0.6, 1.02, 0.01, colormap, 'emissivity SPEC') 
     81    title(month[imo] + ' 2009 - AMSUA 89GHz') 
     82    savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec_AMSUA89_' + month[imo] + '2009.png') 
    8283    # emis SPEC-LAMB 
    8384    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, esl_a_month[imo, :, :], 0., 0.04, 0.001, colormap, 'emissivity LAMB - SPEC') 
    84     title(month[imo] + ' 2009 - AMSUA 23GHz') 
    85     savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/emis_spec-lamb_AMSUA23_' + month[imo] + '2009.png') 
     85    title(month[imo] + ' 2009 - AMSUA 89GHz') 
     86    savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec-lamb_AMSUA89_' + month[imo] + '2009.png') 
    8687 
    8788 
    8889 
    89 ################### 
     90'''################### 
    9091# Read AMSUB data # 
    9192################### 
     
    125126    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, esl_b_month[imo, :, :], 0., 0.05, 0.001, colormap, 'emissivity LAMB - SPEC') 
    126127    title(month[imo] + ' 2009 - AMSUB 89GHz') 
    127     savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/emis_spec-lamb_AMSUB89_' + month[imo] + '2009.png') 
     128    savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/emis_spec-lamb_AMSUB89_' + month[imo] + '2009.png')''' 
    128129 
    129130 
  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_spec_lamb_nadir.py

    r42 r44  
    5959for imo in range (0, M): 
    6060    print 'month: ' + month[imo] 
    61     fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.dat','r') 
     61    fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA89.dat','r') 
    6262    numlines = 0 
    6363    for line in fichier: numlines += 1 
    6464    fichier.close 
    65     fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.dat','r')   
     65    fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA89.dat','r')   
    6666    nbtotal = numlines-1 
    6767    iligne = 0 
     
    185185    e_sl_100_day = zgrid_output 
    186186    esl100[:, :, 0 : month_day[imo], imo] = e_sl_100_day[:, :, :]''' 
    187  
    188  
    189  
    190 ############################################### 
    191 # stack gridded data spec lamb in NETCDF file # 
    192 ############################################### 
    193 print 'stacking of gridded data' 
    194 for imo in range (0, M): 
    195     print 'file for ' + month[imo] 
    196     rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
     187    ############################################### 
     188    # stack gridded data spec lamb in NETCDF file # 
     189    ############################################### 
     190    print 'stacking of gridded data' 
     191    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
    197192    rootgrp.createDimension('longitude', nx) 
    198193    rootgrp.createDimension('latitude', ny) 
  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_spec_lamb_nadir_bis.py

    r41 r44  
    4545# grid data from .dat files # 
    4646############################# 
    47 tu = np.zeros([ny, nx, 31, M], float) 
     47'''tu = np.zeros([ny, nx, 31, M], float) 
    4848td = np.zeros([ny, nx, 31, M], float) 
    4949tbs = np.zeros([ny, nx, 31, M], float) 
    50 tbl = np.zeros([ny, nx, 31, M], float) 
     50tbl = np.zeros([ny, nx, 31, M], float)''' 
    5151es = np.zeros([ny, nx, 31, M], float) 
    5252el = np.zeros([ny, nx, 31, M], float) 
    5353esl = np.zeros([ny, nx, 31, M], float) 
    54 esl00 = np.zeros([ny, nx, 31, M], float) 
     54'''esl00 = np.zeros([ny, nx, 31, M], float) 
    5555esl25 = np.zeros([ny, nx, 31, M], float) 
    5656esl50 = np.zeros([ny, nx, 31, M], float) 
    5757esl75 = np.zeros([ny, nx, 31, M], float) 
    58 esl100 = np.zeros([ny, nx, 31, M], float) 
     58esl100 = np.zeros([ny, nx, 31, M], float)''' 
    5959for imo in range (0, M): 
    6060    print 'month: ' + month[imo] 
    61     fichier=open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA23.dat','r') 
     61    fichier=open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA50.dat','r') 
    6262    numlines = 0 
    6363    for line in fichier: numlines += 1 
    6464    fichier.close 
    65     fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA23.dat','r')   
     65    fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA50.dat','r')   
    6666    nbtotal = numlines-1 
    6767    iligne = 0 
     
    7070    lon = np.zeros([nbtotal],float) 
    7171    e = np.zeros([nbtotal],float) 
    72     ts = np.zeros([nbtotal],float) 
     72    '''ts = np.zeros([nbtotal],float) 
    7373    tup = np.zeros([nbtotal],float) 
    7474    tdn = np.zeros([nbtotal],float) 
     
    7878    tdn_lamb = np.zeros([nbtotal],float) 
    7979    tb_spec = np.zeros([nbtotal],float) 
    80     tb_lamb = np.zeros([nbtotal],float) 
     80    tb_lamb = np.zeros([nbtotal],float)''' 
    8181    e_spec = np.zeros([nbtotal],float) 
    8282    e_lamb = np.zeros([nbtotal],float) 
    8383    e_spec_lamb = np.zeros([nbtotal],float) 
    84     e_sl_00 = np.zeros([nbtotal],float) 
     84    '''e_sl_00 = np.zeros([nbtotal],float) 
    8585    e_sl_25 = np.zeros([nbtotal],float) 
    8686    e_sl_50 = np.zeros([nbtotal],float) 
    8787    e_sl_75 = np.zeros([nbtotal],float) 
    88     e_sl_100 = np.zeros([nbtotal],float) 
     88    e_sl_100 = np.zeros([nbtotal],float)''' 
    8989    while (iligne < nbtotal) : 
    9090         line=fichier.readline() 
     
    9494         lon[iligne] = float(liste[1]) 
    9595         e[iligne] = float(liste[5]) 
    96          ts[iligne] = float(liste[6]) 
     96         '''ts[iligne] = float(liste[6]) 
    9797         tup[iligne] = float(liste[7]) 
    9898         tdn[iligne] = float(liste[8]) 
     
    102102         tdn_lamb[iligne] = float(liste[13]) 
    103103         tb_spec[iligne]  = float(liste[14]) 
    104          tb_lamb[iligne] = float(liste[15]) 
     104         tb_lamb[iligne] = float(liste[15])''' 
    105105         e_spec[iligne] = float(liste[16]) 
    106106         e_lamb[iligne] = float(liste[17]) 
    107107         e_spec_lamb[iligne] = float(liste[18]) 
    108          e_sl_00[iligne] = float(liste[19]) 
     108         '''e_sl_00[iligne] = float(liste[19]) 
    109109         e_sl_25[iligne] = float(liste[20]) 
    110110         e_sl_50[iligne] = float(liste[21]) 
    111111         e_sl_75[iligne] = float(liste[22]) 
    112          e_sl_100[iligne] = float(liste[23]) 
     112         e_sl_100[iligne] = float(liste[23])''' 
    113113         iligne=iligne+1 
    114114    fichier.close() 
    115     print 'tup' 
     115    '''print 'tup' 
    116116    z0 = tup.min() 
    117117    z1 = tup.max() 
     
    136136    zgrid_output, ngrid_output, z2grid_output, sigmagrid_output, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(month_day[imo], jjr, month[imo], lon, lat, tb_lamb, z0, z1, dx, dy) 
    137137    tb_lamb_day = zgrid_output 
    138     tbl[:, :, 0 : month_day[imo], imo] = tb_lamb_day[:, :, :] 
     138    tbl[:, :, 0 : month_day[imo], imo] = tb_lamb_day[:, :, :]''' 
    139139    print 'emis spec' 
    140140    z0 = e_spec.min() 
     
    155155    e_spec_lamb_day = zgrid_output 
    156156    esl[:, :, 0 : month_day[imo], imo] = e_spec_lamb_day[:, :, :] 
    157     print 'emis spec lamb 00' 
     157    '''print 'emis spec lamb 00' 
    158158    z0 = e_sl_00.min() 
    159159    z1 = e_sl_00.max() 
     
    184184    zgrid_output, ngrid_output, z2grid_output, sigmagrid_output, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(month_day[imo], jjr, month[imo], lon, lat, e_sl_100, z0, z1, dx, dy) 
    185185    e_sl_100_day = zgrid_output 
    186     esl100[:, :, 0 : month_day[imo], imo] = e_sl_100_day[:, :, :] 
     186    esl100[:, :, 0 : month_day[imo], imo] = e_sl_100_day[:, :, :]''' 
    187187    ############################################### 
    188188    # stack gridded data spec lamb in NETCDF file # 
     
    190190    print 'stacking of gridded data' 
    191191    print 'file for ' + month[imo] 
    192     rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
     192    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA50_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
    193193    rootgrp.createDimension('longitude', nx) 
    194194    rootgrp.createDimension('latitude', ny) 
     
    197197    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
    198198    nc_day = rootgrp.createVariable('days', 'f', ('days',)) 
    199     nc_tup = rootgrp.createVariable('tup', 'f', ('latitude', 'longitude', 'days')) 
     199    '''nc_tup = rootgrp.createVariable('tup', 'f', ('latitude', 'longitude', 'days')) 
    200200    nc_tdn = rootgrp.createVariable('tdn', 'f', ('latitude', 'longitude', 'days')) 
    201201    nc_tb_spec = rootgrp.createVariable('tb_spec', 'f', ('latitude', 'longitude', 'days')) 
    202     nc_tb_lamb = rootgrp.createVariable('tb_lamb', 'f', ('latitude', 'longitude', 'days')) 
     202    nc_tb_lamb = rootgrp.createVariable('tb_lamb', 'f', ('latitude', 'longitude', 'days'))''' 
    203203    nc_e_spec = rootgrp.createVariable('e_spec', 'f', ('latitude', 'longitude', 'days')) 
    204204    nc_e_lamb = rootgrp.createVariable('e_lamb', 'f', ('latitude', 'longitude', 'days')) 
    205205    nc_e_spec_lamb = rootgrp.createVariable('e_spec_lamb', 'f', ('latitude', 'longitude', 'days')) 
    206     nc_e_sl_00 = rootgrp.createVariable('e_mixed_s00', 'f', ('latitude', 'longitude', 'days')) 
     206    '''nc_e_sl_00 = rootgrp.createVariable('e_mixed_s00', 'f', ('latitude', 'longitude', 'days')) 
    207207    nc_e_sl_25 = rootgrp.createVariable('e_mixed_s25', 'f', ('latitude', 'longitude', 'days')) 
    208208    nc_e_sl_50 = rootgrp.createVariable('e_mixed_s50', 'f', ('latitude', 'longitude', 'days')) 
    209209    nc_e_sl_75 = rootgrp.createVariable('e_mixed_s75', 'f', ('latitude', 'longitude', 'days')) 
    210     nc_e_sl_100 = rootgrp.createVariable('e_mixed_s100', 'f', ('latitude', 'longitude', 'days')) 
     210    nc_e_sl_100 = rootgrp.createVariable('e_mixed_s100', 'f', ('latitude', 'longitude', 'days'))''' 
    211211    nc_lon[:] = xvec 
    212212    nc_lat[:] = yvec 
    213     nc_tup[:] = tu[:, :, 0 : month_day[imo], imo] 
     213    '''nc_tup[:] = tu[:, :, 0 : month_day[imo], imo] 
    214214    nc_tdn[:] = td[:, :, 0 : month_day[imo], imo] 
    215215    nc_tb_spec[:] = tbs[:, :, 0 : month_day[imo], imo] 
    216     nc_tb_lamb[:] = tbl[:, :, 0 : month_day[imo], imo] 
     216    nc_tb_lamb[:] = tbl[:, :, 0 : month_day[imo], imo]''' 
    217217    nc_e_spec[:] = es[:, :, 0 : month_day[imo], imo] 
    218218    nc_e_lamb[:] = el[:, :, 0 : month_day[imo], imo] 
    219219    nc_e_spec_lamb[:] = esl[:, :, 0 : month_day[imo], imo] 
    220     nc_e_sl_00[:] = esl00[:, :, 0 : month_day[imo], imo] 
     220    '''nc_e_sl_00[:] = esl00[:, :, 0 : month_day[imo], imo] 
    221221    nc_e_sl_25[:] = esl25[:, :, 0 : month_day[imo], imo] 
    222222    nc_e_sl_50[:] = esl50[:, :, 0 : month_day[imo], imo] 
    223223    nc_e_sl_75[:] = esl75[:, :, 0 : month_day[imo], imo] 
    224     nc_e_sl_100[:] = esl100[:, :, 0 : month_day[imo], imo] 
     224    nc_e_sl_100[:] = esl100[:, :, 0 : month_day[imo], imo]''' 
    225225    rootgrp.close() 
    226226 
  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/temporal_evolution_AMSUA23-AMSUB89.py

    r43 r44  
    100100 
    101101 
     102#90, 51 
     103#95, 82 
     104 
     105spec_amsua = np.append(es_a[0, 95, 82, 0 : 31], es_a[1, 95, 82, 0 : 28]) 
     106spec_amsub = np.append(es_b[0, 95, 82, 0 : 31], es_b[1, 95, 82, 0 : 28]) 
     107ambig = np.append(amb_ice[0, 95, 82, 0 : 31], amb_ice[1, 95, 82, 0 : 28]) 
     108for imo in range (2, M): 
     109    spec_amsua = np.append(spec_amsua, es_a[imo, 95, 82, 0 : month_day[imo]]) 
     110    spec_amsub = np.append(spec_amsub, es_b[imo, 95, 82, 0 : month_day[imo]]) 
     111    ambig = np.append(ambig, amb_ice[imo, 95, 82, 0 : month_day[imo]]) 
    102112 
    103113 
    104 spec_amsua = np.append(es_a[0, 119, 51, 0 : 31], es_a[1, 119, 51, 0 : 31]) 
    105 spec_amsub = np.append(es_b[0, 119, 51, 0 : 31], es_b[1, 119, 51, 0 : 31]) 
    106 ambig = np.append(amb_ice[0, 119, 51, 0 : 31], amb_ice[1, 119, 51, 0 : 31]) 
    107 for imo in range (2, M): 
    108     spec_amsua = np.append(spec_amsua, es_a[imo, 119, 51, 0 : month_day[imo]]) 
    109     spec_amsub = np.append(spec_amsub, es_b[imo, 119, 51, 0 : month_day[imo]]) 
    110     ambig = np.append(ambig, amb_ice[imo, 119, 51, 0 : month_day[imo]]) 
    111  
     114# seperately AMSUA and AMSUB 
    112115ion() 
    113116figure() 
    114 plot(spec_amsua[nonzero(spec_amsua != 0.)], 'b', label = 'AMSUA 23GHz') 
    115 plot(spec_amsub[nonzero(spec_amsub != 0.)], 'g', label = 'AMSUB 89GHz') 
    116 plot(np.arange(0, 368, 1)[nonzero(ambig == 1.)], spec_amsua[nonzero(ambig == 1.)], 'ob', label = 'obs in ambiguous area') 
    117 plot(np.arange(0, 368, 1)[nonzero(ambig == 1.)], spec_amsub[nonzero(ambig == 1.)], 'og', label = 'obs in ambiguous area') 
     117plot(np.arange(0, 365, 1), spec_amsua, 'b', label = 'AMSUA 23GHz') 
     118plot(np.arange(0, 365, 1), spec_amsub, 'g', label = 'AMSUB 89GHz') 
     119plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_amsua[nonzero(ambig == 1.)], 'ob', label = 'ambiguous area') 
     120plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_amsub[nonzero(ambig == 1.)], 'og', label = 'ambiguous area') 
    118121legend(loc = 3) 
     122xlim(-3, 368) 
     123xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 334]), month, rotation = 20) 
     124grid() 
     125 
     126# bias AMSUA-AMSUB 
     127figure() 
     128plot(np.arange(0, 365, 1), spec_amsua - spec_amsub, 'r', label = 'AMSUA23 - AMSUB89') 
     129plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], (spec_amsua - spec_amsub)[nonzero(ambig == 1.)], 'or', label = 'ambiguous area') 
     130plot(np.arange(0, 365, 1), np.zeros([365], float), 'k') 
     131legend(loc = 3) 
     132xlim(-3, 368) 
     133xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 335]), month, rotation = 20) 
     134 
     135 
     136# anomaly with a 15 day time range  
     137spec_anom = np.zeros([365], float) 
     138spec_diff = spec_amsua - spec_amsub 
     139for ii in range (15, 365 - 15): 
     140    spec_anom[ii] = (spec_diff[ii]) - mean(spec_diff[ii - 15 : ii + 15][nonzero(isnan(spec_diff[ii - 15 : ii + 15]) == False)]) 
     141 
     142figure() 
     143plot(np.arange(0, 365, 1), spec_anom, 'r') 
     144plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_anom[nonzero(ambig == 1.)], 'or', label = 'ambiguous area') 
     145plot(np.arange(0, 365, 1), np.zeros([365], float), 'k') 
     146legend(loc = 3) 
     147xlim(-3, 368) 
     148xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 335]), month, rotation = 20) 
     149 
     150 
     151 
     152 
    119153 
    120154# plot of anomaly (espec amsub - espec amsua) on a 15 day sliding base (15 day correct to simulate season transition) 
Note: See TracChangeset for help on using the changeset viewer.