Changeset 45


Ignore:
Timestamp:
07/24/14 17:12:02 (10 years ago)
Author:
lahlod
Message:

new_scripts

Location:
trunk/src/scripts_Laura/ARCTIC/Travail_CEN
Files:
2 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/discrim_rate_espec-elamb_espec.py

    r44 r45  
    5050    ### AMSUA23 ### 
    5151    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') 
     52    fichier = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA50_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    5353    xdist = fichier.variables['longitude'][:] 
    5454    ydist = fichier.variables['latitude'][:] 
    … …  
    9696for imo in range (0, M): 
    9797    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') 
     98    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, monthly_ratio[imo, :, :], 0., 20., 0.1, colormap, 'emissivity ratio [(LAMB - SPEC)/SPEC]*100') 
     99    title(month[imo] + ' 2009 - AMSUA 50GHz') 
     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_AMSUA50_' + month[imo] + '2009.png') 
    101101 
    102102 
    … …  
    112112for imo in range (0, M): 
    113113    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') 
     114    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA50_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
    115115    rootgrp.createDimension('longitude', len(xdist)) 
    116116    rootgrp.createDimension('latitude', len(ydist)) 

  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/ice_delimit_AMSUB_data.py

    r44 r45  
    1717import map_cartesian_grid 
    1818 
     19 
     20 
     21 
    1922MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) 
    2023month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER']) 
    … …  
    4851spec_emis_vec = np.zeros([M, ny * nx], float) 
    4952lamb_emis_vec = np.zeros([M, ny * nx], float) 
     53print 'read data from files' 
    5054for imo in range (0, M): 
    5155    print 'month: ' + month[imo] 
    5256    ### emissivity ### 
    5357    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') 
     58    fichier_emis = 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') 
    5559    xdist = fichier_emis.variables['longitude'][:] 
    5660    ydist = fichier_emis.variables['latitude'][:] 
    … …  
    6670    ### monthly emissivity ratio ### 
    6771    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') 
     72    fichier_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA89_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    6973    xdist = fichier_ratio.variables['longitude'][:] 
    7074    ydist = fichier_ratio.variables['latitude'][:] 
    … …  
    8993 
    9094c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) 
     95limit_coef_spec = np.array([0.6, 0.7, 0.8]) # i1 = 23GHz / i2 = 50GHz / i3 = 89GHz 
     96limit_coef_lamb = np.array([0.6, 0.8, 0.8]) # i1 = 23GHz / i2 = 50GHz / i3 = 89GHz 
     97idata = 2 
     98print 'distribution of emissivity values' 
    9199################################### 
    92100# distribution of SPEC emissivity # 
    93101################################### 
     102print 'spec' 
    94103hist_vals_spec = np.zeros([M, 50], float) 
    95104corresp_emis_spec = np.zeros([M, 50], float) 
    … …  
    112121fontP.set_size('small') 
    113122legend(loc = 2, prop = fontP) 
    114 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA23_JANUARY-JUNE_2009.png') 
     123#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA50_JANUARY-JUNE_2009.png') 
    115124## plot six following months of spec emissivity histograms ## 
    116125figure() 
    … …  
    123132ylabel('frequency of occurence') 
    124133legend(loc = 9, prop = fontP) 
    125 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA23_JULY-DECEMBER_2009.png') 
     134#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA50_JULY-DECEMBER_2009.png') 
    126135 
    127136# find the emissivity value corresponding to the threshold of ice/no_ice limit 
    … …  
    129138for imo in range (0, M): 
    130139    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 
     140    bb = np.where((corresp_emis_spec[imo, :] < limit_coef_spec[idata]))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 
    132141    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 
    133142    cc = np.where(hist_vals_spec[imo, bb] == max(hist_vals_spec[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 
    … …  
    138147figure() 
    139148plot(np.arange(0, M, 1), emis_lim_spec, 'b-+', label = 'emis SPEC') 
    140 ylabel('emissivity SPEC threshold - AMSUA 23GHz') 
     149ylabel('emissivity SPEC threshold') 
    141150grid() 
    142151xticks(np.arange(0, M, 1), month, rotation = 25) 
    143152legend(loc = 2) 
    144153xlim(-1, M) 
    145 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_SPEC_AMSUA23_2009.png') 
     154#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_SPEC_AMSUA50_2009.png') 
    146155 
    147156 
    … …  
    149158# distribution of LAMB emissivity # 
    150159################################### 
     160print 'lamb' 
    151161hist_vals_lamb = np.zeros([M, 50], float) 
    152162corresp_emis_lamb = np.zeros([M, 50], float) 
    … …  
    166176ylabel('frequency of occurence') 
    167177legend(loc = 2, prop = fontP) 
    168 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA23_JANUARY-JUNE_2009.png') 
     178plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA50_JANUARY-JUNE_2009.png') 
    169179## plot six following months of spec emissivity histograms ## 
    170180figure() 
    … …  
    177187ylabel('frequency of occurence') 
    178188legend(loc = 9, prop = fontP) 
    179 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA23_JULY-DECEMBER_2009.png') 
     189#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA50_JULY-DECEMBER_2009.png') 
    180190 
    181191# find the emissivity value corresponding to the threshold of ice/no_ice limit 
    … …  
    183193for imo in range (0, M): 
    184194    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 
     195    bb = np.where((corresp_emis_lamb[imo, :] < limit_coef_lamb[idata]))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 
    186196    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 
    187197    cc = np.where(hist_vals_lamb[imo, bb] == max(hist_vals_lamb[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 
    … …  
    192202figure() 
    193203plot(np.arange(0, M, 1), emis_lim_lamb, 'b-+', label = 'emis LAMB') 
    194 ylabel('emissivity LAMB threshold - AMSUB 23GHz') 
     204ylabel('emissivity LAMB threshold') 
    195205grid() 
    196206xticks(np.arange(0, M, 1), month, rotation = 25) 
    197207legend(loc = 2) 
    198208xlim(-1, M) 
    199 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_LAMB_AMSUA23_2009.png') 
     209#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_LAMB_AMSUA50_2009.png') 
    200210 
    201211 
    … …  
    204214# distribution of emissivity rate # 
    205215################################### 
     216print 'rate' 
    206217hist_vals_rate = np.zeros([M, 50], float) 
    207218corresp_emis_rate = np.zeros([M, 50], float) 
    … …  
    221232ylabel('frequency of occurence') 
    222233legend(prop = fontP) 
    223 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA23_JANUARY-JUNE_2009.png') 
     234#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA50_JANUARY-JUNE_2009.png') 
    224235## plot six following months of spec emissivity histograms ## 
    225236figure() 
    … …  
    232243ylabel('frequency of occurence') 
    233244legend(loc = 9, prop = fontP) 
    234 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA23_JULY-DECEMBER_2009.png') 
     245#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA50_JULY-DECEMBER_2009.png') 
    235246 
    236247# no definition of threshold for ICE / NO_ICE delimitation because no apparent distinction signal // use of this signal for internal ice classification 
    … …  
    263274plot(np.arange(0, M, 1), emis_lim_spec, 'b-+', label = 'emis SPEC') 
    264275plot(np.arange(0, M, 1), emis_lim_lamb, 'g-+', label = 'emis LAMB') 
    265 ylabel('emissivity threshold - AMSUA 23GHz') 
     276ylabel('emissivity threshold') 
    266277grid() 
    267278xticks(np.arange(0, M, 1), month, rotation = 25) 
    268279legend(loc = 2, prop = fontP) 
    269280xlim(-1, M) 
    270 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_spec_and_lamb_AMSUA23_2009.png') 
     281#plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_spec_and_lamb_AMSUA50_2009.png') 
    271282 
    272283 
    … …  
    276287# delimitation ice - no_ice with emissivity or ratio threshold # 
    277288################################################################ 
     289print 'definition of ice extent'  
    278290x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    279291x_coast = x_ind 
    280292y_coast = y_ind 
    281293z_coast = z_ind 
    282 ice_zone = np.zeros([M, 151, 139], float) 
     294#### SPEC emiss #### 
     295print 'spec' 
     296ice_zone_spec = np.zeros([M, 151, 139], float) 
    283297for imo in range (0, M): 
    284298    print 'month: ' + str(month[imo]) 
    … …  
    286300        for ilon in range (0, 139): 
    287301            if (isnan(emis_spec_month[imo, ilat, ilon]) == True): 
    288                 ice_zone[imo, ilat, ilon] = NaN  
     302                ice_zone_spec[imo, ilat, ilon] = NaN  
    289303            else: 
    290304                if (emis_spec_month[imo, ilat, ilon] <= emis_lim_spec[imo]): # use the monthly SPEC emissivity threshold  
    291                     ice_zone[imo, ilat, ilon] = 0. 
     305                    ice_zone_spec[imo, ilat, ilon] = 0. 
    292306                else: 
    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  
     307                    ice_zone_spec[imo, ilat, ilon] = 1.    
     308    '''map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_spec[imo, :, :], -1, 2., 1., colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_SPEC > ' + str(emis_lim_spec[imo])[0:6] + ')') 
     309    title(str(month[imo]) + ' 2009 - AMSUA 89GHz') 
     310    plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/AMSUA89_ice_emis_spec_thresh' + '_' + month[imo] + '2009.png')''' 
     311 
     312#### LAMB emiss #### 
     313print 'lamb' 
     314ice_zone_lamb = np.zeros([M, 151, 139], float) 
     315for imo in range (0, M): 
     316    print 'month: ' + str(month[imo]) 
     317    for ilat in range (0, 151): 
     318        for ilon in range (0, 139): 
     319            if (isnan(emis_lamb_month[imo, ilat, ilon]) == True): 
     320                ice_zone_lamb[imo, ilat, ilon] = NaN  
     321            else: 
     322                if (emis_lamb_month[imo, ilat, ilon] <= emis_lim_lamb[imo]): # use the monthly SPEC emissivity threshold  
     323                    ice_zone_lamb[imo, ilat, ilon] = 0. 
     324                else: 
     325                    ice_zone_lamb[imo, ilat, ilon] = 1.    
     326    '''map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_lamb[imo, :, :], -1, 2., 1., colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_LAMB > ' + str(emis_lim_lamb[imo])[0:6] + ')') 
     327    title(str(month[imo]) + ' 2009 - AMSUA 89GHz') 
     328    plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/AMSUA89_ice_emis_lamb_thresh' + '_' + month[imo] + '2009.png')''' 
    299329 
    300330 
    … …  
    304334########################### 
    305335# nb of pixels near pole = 926 
     336print 'calculation of ice area' 
    306337pix_area = 40. * 40. 
    307 nb_pix = np.zeros([M], float) 
    308 total_ice_area = np.zeros([M], float) 
    309 for 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  
    315 spec_ice_area = total_ice_area 
    316  
    317  
    318 figure() 
    319 plot(lamb_ice_area, 'g', label = 'lamb') 
    320 plot(spec_ice_area, 'b', label = 'spec') 
    321 plot(total_monthly_ice_area_osisaf, 'r', label = 'OSISAF') 
    322 ylabel('total ice area (square km)') 
    323 grid() 
    324 xticks(np.arange(0, M, 1), month, rotation = 25) 
    325 legend(loc = 3, prop = fontP) 
    326 title('AMSUA 23GHz') 
    327 xlim(-1, M) 
    328 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/total_ice_area_AMSUA23_OSISAF_2009.png') 
     338#### SPEC #### 
     339print 'spec' 
     340nb_pix_spec = np.zeros([M], float) 
     341total_ice_area_spec = np.zeros([M], float) 
     342for imo in range (0, M): 
     343    ice_spec = reshape(ice_zone_spec[imo, :, :], size(ice_zone_spec[imo, :, :])) 
     344    nb_pix_spec[imo] = len(np.where(ice_spec == 1.)[0]) 
     345    total_ice_area_spec[imo] = (pix_area * nb_pix_spec[imo]) + (926. * pix_area) 
     346 
     347 
     348#### LAMB #### 
     349print 'lamb' 
     350nb_pix_lamb = np.zeros([M], float) 
     351total_ice_area_lamb = np.zeros([M], float) 
     352for imo in range (0, M): 
     353    ice_lamb = reshape(ice_zone_lamb[imo, :, :], size(ice_zone_lamb[imo, :, :])) 
     354    nb_pix_lamb[imo] = len(np.where(ice_lamb == 1.)[0]) 
     355    total_ice_area_lamb[imo] = (pix_area * nb_pix_lamb[imo]) + (926. * pix_area) 
     356 
     357 
     358 
     359########################################## STACK DATA ################################################### 
     360 
    329361 
    330362########################## 
    331363# stack of ice area data # 
    332364########################## 
    333 # stack in .dat file OSISAF data 
    334 print 'start stacking in .dat file' 
    335 data_osisaf = open ('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_daily_ice_no-ice_2009.dat', 'a') 
    336 for 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  
    345 data_osisaf.close() 
    346  
    347  
    348365# stack in .dat file  
    349366print 'start stacking in .dat file' 
    350 data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/2009.dat', 'a') 
    351 for 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' % { 
     367data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/AMSUA89_data_classification_parameters_ice_no-ice_2009.dat', 'a') 
     368for imo in range (0, M): 
     369    for ii in range (0, 50): 
     370        data_classif.write(('%(months)10s    %(hist_vals_spec)10.5f    %(corresp_emis_spec)10.5f    %(hist_vals_lamb)10.5f    %(corresp_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' % { 
    354371                                            'months':month[imo], 
    355372                                            'hist_vals_spec':hist_vals_spec[imo, ii], 
    … …  
    358375                                            'corresp_emis_lamb':corresp_emis_lamb[imo, ii], 
    359376                                            'hist_vals_rate':hist_vals_rate[imo, ii], 
     377                                            'corresp_emis_rate':corresp_emis_rate[imo, ii], 
    360378                                            'emis_lim_spec':emis_lim_spec[imo], 
    361379                                            'emis_lim_lamb':emis_lim_lamb[imo], 
    362                                             'spec_ice_area':spec_ice_area[imo], 
    363                                             'lamb_ice_area':lamb_ice_area[imo], 
     380                                            'spec_ice_area':total_ice_area_spec[imo], 
     381                                            'lamb_ice_area':total_ice_area_lamb[imo], 
    364382                                            })) 
    365383 
    366384data_classif.close() 
    367  
     385''' 
    368386# stack in netcdf file 
    369387print 'start stacking' 
    370388for imo in range (0, M): 
    371389    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') 
     390    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + month[imo] + '2009_AMSUA50_spec_lamb_thresholds.nc', 'w', format='NETCDF3_CLASSIC') 
    373391    rootgrp.createDimension('longitude', len(xdist)) 
    374392    rootgrp.createDimension('latitude', len(ydist)) 
    375393    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
    376394    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
    377     nc_ice = rootgrp.createVariable('ice_area', 'f', ('latitude', 'longitude')) 
     395    nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude')) 
     396    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude')) 
    378397    nc_lon[:] = xdist 
    379398    nc_lat[:] = ydist 
    380     nc_ice[:] = ice_zone[imo, :, :] 
     399    nc_ice_spec[:] = ice_zone_spec[imo, :, :] 
     400    nc_ice_lamb[:] = ice_zone_lamb[imo, :, :] 
    381401    rootgrp.close() 
    382  
    383  
    384 ############################################################################################################# 
    385  
    386 ######################################### 
    387 # comparison with OSISAF sea ice extent # 
    388 ######################################### 
    389 osi_ice = np.zeros([M, 31, 151, 139]) 
    390 for 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)): 
    400             for ijr in range (0, month_day[imo]): 
    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)): 
    402                     osi_ice[imo, ijr, ilat, ilon] = 1. 
    403                 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): 
    408                             osi_ice[imo, ijr, ilat, ilon] = NaN 
    409                         else:  
    410                             osi_ice[imo, ijr, ilat, ilon] = 0. 
    411  
    412  
    413 # test: 
    414 colormap = colors.ListedColormap(['cyan','blue']) 
    415 map_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 ################################## 
    421 pix_area = 40. * 40. 
    422 nb_pix_osisaf = np.zeros([M, 31], float) 
    423 total_ice_area_osisaf = np.zeros([M, 31], float) 
    424 total_monthly_ice_area_osisaf = np.zeros([M], float) 
    425 for 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  
     402''' 
     403 
     404 
     405 
     406 
     407''' 
    436408################################### 
    437409# plot time evolution of ice area # 
    … …  
    448420grid() 
    449421plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_pix_area/total_ice_area_AMSUB_SPEC_LAMB_OSISAF_2009.png') 
    450  
    451  
    452  
    453  
    454  
    455  
    456  
    457  
     422''' 
     423 
     424 
     425 
     426 
     427 
     428 
     429 

  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_cartesian_gridded_AMSUA23.py

    r44 r45  
    3434x0 = -3000. # min limit of grid 
    3535x1 = 2500. # max limit of grid 
    36 dx = 40. 
     36dx = 40. # resolution for AMSUA data = 40km // resultion for AMSUB data = 20km 
    3737xvec = np.arange(x0, x1+dx, dx) 
    3838nx = len(xvec)  
    3939y0 = -3000. # min limit of grid 
    4040y1 = 3000. # max limit of grid 
    41 dy = 40. 
     41dy = 40.  # resolution for AMSUA data = 40km // resultion for AMSUB data = 20km 
    4242yvec = np.arange(y0, y1+dy, dy) 
    4343ny = len(yvec) 
    … …  
    5353for imo in range (0, M): 
    5454    print 'open file ' + str(month[imo]) 
    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') 
     55    file_amsua = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA50_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    5656    lon_amsua = file_amsua.variables['longitude'][:] 
    5757    lat_amsua = file_amsua.variables['latitude'][:] 
    … …  
    7878for imo in range (0, M): 
    7979    # emiss SPEC 
    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') 
     80    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, es_a_month[imo, :, :], 0.4, 1.02, 0.01, colormap, 'emissivity SPEC') 
     81    title(month[imo] + ' 2009 - AMSUA 50GHz') 
     82    savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec_AMSUA50_' + month[imo] + '2009.png') 
    8383    # emis SPEC-LAMB 
    84     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') 
    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') 
     84    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, esl_a_month[imo, :, :], 0., 0.12, 0.001, colormap, 'emissivity LAMB - SPEC') 
     85    title(month[imo] + ' 2009 - AMSUA 50GHz') 
     86    savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec-lamb_AMSUA50_' + month[imo] + '2009.png') 
    8787 
    8888 

  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_osisaf_ice_types.py

    r42 r45  
    3131x0 = -3000. # min limit of grid 
    3232x1 = 2500. # max limit of grid 
    33 dx = 40. 
     33dx = 10. 
    3434xvec = np.arange(x0, x1+dx, dx) 
    3535nx = len(xvec)  
    3636y0 = -3000. # min limit of grid 
    3737y1 = 3000. # max limit of grid 
    38 dy = 40. 
     38dy = 10. 
    3939yvec = np.arange(y0, y1+dy, dy) 
    4040ny = len(yvec) 
    … …  
    7676    ################################################ 
    7777    print 'open file month' + str(month[imo]) 
    78     rootgrp = 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', 'w', format='NETCDF3_CLASSIC') 
     78    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_cartesian_grid/OSISAF_ice_types_cartes_grid_res-10_' + str(month[imo]) + '_2009.nc', 'w', format='NETCDF3_CLASSIC') 
    7979    rootgrp.createDimension('longitude', len(xvec)) 
    8080    rootgrp.createDimension('latitude', len(yvec)) 

  • trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_spec_lamb_nadir.py

    r44 r45  
    3131x0 = -3000. # min limit of grid 
    3232x1 = 2500. # max limit of grid 
    33 dx = 40. 
     33dx = 20. 
    3434xvec = np.arange(x0, x1+dx, dx) 
    3535nx = len(xvec)  
    3636y0 = -3000. # min limit of grid 
    3737y1 = 3000. # max limit of grid 
    38 dy = 40. 
     38dy = 20. 
    3939yvec = np.arange(y0, y1+dy, dy) 
    4040ny = len(yvec) 
    … …  
    5757esl75 = np.zeros([ny, nx, 31, M], float) 
    5858esl100 = np.zeros([ny, nx, 31, M], float)''' 
    59 for imo in range (0, M): 
     59for imo in range (10, 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_AMSUA89.dat','r') 
     61    fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.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_AMSUA89.dat','r')   
     65    fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.dat','r')   
    6666    nbtotal = numlines-1 
    6767    iligne = 0 
    … …  
    189189    ############################################### 
    190190    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') 
     191    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_res-20_monthly_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 
    192192    rootgrp.createDimension('longitude', nx) 
    193193    rootgrp.createDimension('latitude', ny) 
Note: See TracChangeset for help on using the changeset viewer.