Changeset 54 for trunk/src/scripts_Laura


Ignore:
Timestamp:
08/13/14 19:00:17 (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/cartesian_grid_monthly.py

    r41 r54  
    88from mpl_toolkits.basemap import Basemap 
    99from mpl_toolkits.basemap import shiftgrid, cm 
    10 import draw_map 
     10 
    1111 
    1212#################################################### 

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

    r49 r54  
    4444################################################################################################################## 
    4545# We devide the loop in two steps :  
    46 # - first loop concerns all months except for AUGUST and SEPTEMBER - use of AMSUA23GHz SPEC emissivity to seperate ice from no-ice zones  
    47 # - second loop concerns AUGUST and SEPTEMBER - use of AMSUA30GHz SPEC emissivity to seperate ice from no_ice zones 
     46# - first loop concerns JANUARY, FEBRUARY, MARCH, APRIL, SEPTEMBER, OCTOBER, NOVEMBER, DECEMBER - use of AMSUA23GHz SPEC emissivity to seperate ice from no-ice zones  
     47# - second loop concerns MAY, JUNE, JULY, AUGUST - use of AMSUA89GHz SPEC emissivity to seperate ice from no_ice zones 
    4848################################################################################################################## 
    49 frequ = 23 # apply threshold on this frequency 
    50 # open .dat file to stack data (see end of loop) 
    51 #data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/AMSUA'+str(frequ)+'_data_classification_parameters_ice_no-ice_with_AMSUA23-and-30-spec_2009.dat', 'a') 
     49frequ = 89 # apply threshold on this frequency 
     50''' 
     51#open .dat file to stack data (see end of loop) 
     52data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/AMSUA'+str(frequ)+'_data_classification_parameters_ice_no-ice_with_AMSUA23-and-30-spec_2009.dat', 'a') 
    5253bin = 50 
    53  
    54  
    55 # monthly mean parameter (2D-array) on ARCTIC area 
    56 spec_month = np.zeros([M, ny, nx], float) 
    57 lamb_month = np.zeros([M, ny, nx], float) 
    58 diff_month = np.zeros([M, ny, nx], float) 
    59 ratio_month = np.zeros([M, ny, nx], float) 
    60 # monthly mean parameter (2D-array) on ARCTIC SEA ICE area 
    61 spec_ice = np.zeros([M, ny, nx], float) 
    62 lamb_ice = np.zeros([M, ny, nx], float) 
    63 diff_ice = np.zeros([M, ny, nx], float) 
    64 ratio_ice = np.zeros([M, ny, nx], float) 
     54''' 
     55 
     56 
     57# daily parameter (2D-array) on ARCTIC area 
     58emis_spec = np.zeros([M, ny, nx, 31], float) 
     59emis_lamb = np.zeros([M, ny, nx, 31], float) 
     60emis_diff = np.zeros([M, ny, nx, 31], float) 
     61emis_ratio = np.zeros([M, ny, nx, 31], float) 
     62 
     63# daily parameter (2D-array) on ARCTIC SEA ICE area 
     64daily_spec_ice = np.zeros([M, ny, nx, 31], float) 
     65daily_lamb_ice = np.zeros([M, ny, nx, 31], float) 
     66daily_diff_ice = np.zeros([M, ny, nx, 31], float) 
     67daily_ratio_ice = np.zeros([M, ny, nx, 31], float) 
     68 
     69''' 
    6570# monthly mean parameter (1D-array) on ARCTIC SEA ICE area transformed into vector 
    6671spec_vec = np.zeros([M, ny * nx], float) 
    … …  
    6873diff_vec = np.zeros([M, ny * nx], float) 
    6974ratio_vec = np.zeros([M, ny * nx], float) 
     75 
    7076# histogram distribution (intensity of occurence) of parameter in SEA ICE area (1D-array, bins = 200) 
    7177hist_vals_spec = np.zeros([M, bin], float) 
    … …  
    7379hist_vals_diff = np.zeros([M, bin], float) 
    7480hist_vals_ratio = np.zeros([M, bin], float) 
     81 
    7582# histogram distribution (emissivity value) of parameter in SEA ICE area (1D-array, bins = 200) 
    7683corresp_emis_spec = np.zeros([M, bin], float) 
    … …  
    7885corresp_emis_diff = np.zeros([M, bin], float) 
    7986corresp_emis_ratio = np.zeros([M, bin], float) 
    80 months1 = np.array([0, 1, 2, 3, 4, 5, 6, 9, 10, 11]) # use AMSUA 23GHz to delimit ice/no_ice for all months except for AUGUST and SEPTEMBER 
     87''' 
     88months1 = np.array([0, 1, 2, 3, 8, 9, 10, 11]) # use AMSUA 23GHz to delimit ice/no_ice for JANUARY, FEBRUARY, MARCH, APRIL, SEPTEMBER, OCTOBER, NOVEMBER, DECEMBER 
    8189for imo in months1: 
    8290    print 'month ' + month[imo] 
    … …  
    94102    print 'open emissivity to sub_classify file' 
    95103    ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB) 
    96     fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
     104    fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    97105    day = fichier_e.variables['days'][:] 
    98     emis_spec = fichier_e.variables['e_spec'][:] 
    99     emis_lamb = fichier_e.variables['e_lamb'][:] 
    100     emis_diff = fichier_e.variables['e_spec_lamb'][:] 
     106    emis_spec[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec'][:] 
     107    emis_lamb[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_lamb'][:] 
     108    emis_diff[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec_lamb'][:] 
    101109    fichier_e.close() 
    102     for ilon in range (0, nx): 
    103         for ilat in range (0, ny): 
    104             spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, :][nonzero(isnan(emis_spec[ilat, ilon, :]) == False)]) 
    105             lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, :][nonzero(isnan(emis_lamb[ilat, ilon, :]) == False)]) 
    106             diff_month[imo, ilat, ilon] = mean(emis_diff[ilat, ilon, :][nonzero(isnan(emis_diff[ilat, ilon, :]) == False)]) 
    107     ## data of emis RATIO 
    108     fichier_r = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    109     ratio_month[imo, :, :] = fichier_r.variables['emis_ratio'][:] 
    110     fichier_r.close() 
    111     print 'compute matrix of parameter on SEA ICE area' 
    112     for ilon in range (0, nx): 
    113         for ilat in range (0, ny): 
    114             if (isnan(spec_lim[ilat, ilon]) == True): 
    115                 spec_ice[imo, ilat, ilon] = NaN 
    116                 lamb_ice[imo, ilat, ilon] = NaN 
    117                 diff_ice[imo, ilat, ilon] = NaN 
    118                 ratio_ice[imo, ilat, ilon] = NaN 
    119             else: 
    120                 spec_ice[imo, ilat, ilon] = spec_month[imo, ilat, ilon] 
    121                 lamb_ice[imo, ilat, ilon] = lamb_month[imo, ilat, ilon] 
    122                 diff_ice[imo, ilat, ilon] = diff_month[imo, ilat, ilon] 
    123                 ratio_ice[imo, ilat, ilon] = ratio_month[imo, ilat, ilon] 
     110    # calculate emis ratio daily 
     111    for ijr in range (0, month_day[imo]): 
     112        for ilon in range (0, nx): 
     113            for ilat in range (0, ny): 
     114                emis_ratio[imo, ilat, ilon, ijr] = ((emis_lamb[imo, ilat, ilon, ijr] - emis_spec[imo, ilat, ilon, ijr]) / emis_spec[imo, ilat, ilon, ijr]) * 100. 
     115                if (isnan(spec_lim[ilat, ilon]) == True): 
     116                    daily_spec_ice[imo, ilat, ilon, ijr] = NaN 
     117                    daily_lamb_ice[imo, ilat, ilon, ijr] = NaN 
     118                    daily_diff_ice[imo, ilat, ilon, ijr] = NaN 
     119                    daily_ratio_ice[imo, ilat, ilon, ijr] = NaN 
     120                else: 
     121                    daily_spec_ice[imo, ilat, ilon, ijr] = emis_spec[imo, ilat, ilon, ijr] 
     122                    daily_lamb_ice[imo, ilat, ilon, ijr] = emis_lamb[imo, ilat, ilon, ijr] 
     123                    daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr] 
     124                    daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr] 
     125    ''' 
     126    # ATTENTION : previous part of script has been modified ==> cannot be applied directly to this following part of script (modification of arrays and names.... 
    124127    print 'compute SPEC distribution' 
    125128    ######## 
    … …  
    163166    print 'start stacking in .dat file' 
    164167    #data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/AMSUB'+str(frequ)+'_data_classification_parameters_ice_no-ice_with_AMSUA23-spec_2009.dat', 'a') 
    165     '''for ii in range (0, bin): 
     168    for ii in range (0, bin): 
    166169        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_diff)10.5f    %(corresp_emis_diff)10.5f    %(hist_vals_rate)10.5f    %(corresp_emis_rate)10.5f    \n' % { 
    167170                                            'months':month[imo], 
    … …  
    174177                                            'hist_vals_rate':hist_vals_ratio[imo, ii], 
    175178                                            'corresp_emis_rate':corresp_emis_ratio[imo, ii], 
    176                                             }))''' 
     179                                            })) 
     180    ''' 
    177181    ######################## 
    178182    # stack in netcdf file # 
    179183    ########################  
    180184    print 'stack in file month ' + str(month[imo]) 
    181     rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-30_' + month[imo] + '2009_AMSUA' + str(frequ) + '_spec_lamb_thresholds.nc', 'w', format='NETCDF3_CLASSIC') 
    182     rootgrp.createDimension('longitude', len(xvec)) 
    183     rootgrp.createDimension('latitude', len(yvec)) 
     185    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUB' + str(frequ) + '_spec_thresholds.nc', 'w', format='NETCDF3_CLASSIC') 
     186    rootgrp.createDimension('longitude', nx) 
     187    rootgrp.createDimension('latitude', ny) 
     188    rootgrp.createDimension('days', month_day[imo]) 
    184189    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
    185190    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
    186     nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude')) 
    187     nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude')) 
    188     nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude')) 
    189     nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude')) 
     191    nc_days = rootgrp.createVariable('days', 'f', ('days',)) 
     192    nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     193    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     194    nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     195    nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude', 'days')) 
    190196    nc_lon[:] = xvec 
    191197    nc_lat[:] = yvec 
    192     nc_ice_spec[:] = spec_ice[imo, :, :] 
    193     nc_ice_lamb[:] = lamb_ice[imo, :, :] 
    194     nc_ice_diff[:] = diff_ice[imo, :, :] 
    195     nc_ice_ratio[:] = ratio_ice[imo, :, :] 
     198    nc_days[:] = np.arange(0, month_day[imo]) 
     199    nc_ice_spec[:] = daily_spec_ice[imo, :, :, 0 : month_day[imo]] 
     200    nc_ice_lamb[:] = daily_lamb_ice[imo, :, :, 0 : month_day[imo]] 
     201    nc_ice_diff[:] = daily_diff_ice[imo, :, :, 0 : month_day[imo]] 
     202    nc_ice_ratio[:] = daily_ratio_ice[imo, :, :, 0 : month_day[imo]] 
    196203    rootgrp.close() 
    197204 
    198205 
    199 months2 = np.array([7, 8])# use AMSUA 30GHz to delimit ice/no_ice for AUGUST and SEPTEMBER 
     206 
     207 
     208months2 = np.array([4, 5, 6, 7])# use AMSUA 89GHz to delimit ice/no_ice for MAY, JUNE, JULY, AUGUST 
    200209for imo in months2: 
    201210    print 'month ' + month[imo] 
    … …  
    204213    ################################################################################## 
    205214    print 'open threshold file' 
    206     fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA23_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 
     215    fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA89_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 
    207216    spec_lim = fichier_emis.variables['spec_ice_area'][:] 
    208217    #lamb_lim = fichier_emis.variables['lamb_ice_area'][:] 
    … …  
    213222    print 'open emissivity to sub_classify file' 
    214223    ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB) 
    215     fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
     224    fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    216225    day = fichier_e.variables['days'][:] 
    217     emis_spec = fichier_e.variables['e_spec'][:] 
    218     emis_lamb = fichier_e.variables['e_lamb'][:] 
    219     emis_diff = fichier_e.variables['e_spec_lamb'][:] 
     226    emis_spec[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec'][:] 
     227    emis_lamb[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_lamb'][:] 
     228    emis_diff[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec_lamb'][:] 
    220229    fichier_e.close() 
    221     for ilon in range (0, nx): 
    222         for ilat in range (0, ny): 
    223             spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, :][nonzero(isnan(emis_spec[ilat, ilon, :]) == False)]) 
    224             lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, :][nonzero(isnan(emis_lamb[ilat, ilon, :]) == False)]) 
    225             diff_month[imo, ilat, ilon] = mean(emis_diff[ilat, ilon, :][nonzero(isnan(emis_diff[ilat, ilon, :]) == False)]) 
    226     ## data of emis RATIO 
    227     fichier_r = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    228     ratio_month[imo, :, :] = fichier_r.variables['emis_ratio'][:] 
    229     fichier_r.close() 
    230     print 'compute matrix of parameter on SEA ICE area' 
    231     for ilon in range (0, nx): 
    232         for ilat in range (0, ny): 
    233             if (isnan(spec_lim[ilat, ilon]) == True): 
    234                 spec_ice[imo, ilat, ilon] = NaN 
    235                 lamb_ice[imo, ilat, ilon] = NaN 
    236                 diff_ice[imo, ilat, ilon] = NaN 
    237                 ratio_ice[imo, ilat, ilon] = NaN 
    238             else: 
    239                 spec_ice[imo, ilat, ilon] = spec_month[imo, ilat, ilon] 
    240                 lamb_ice[imo, ilat, ilon] = lamb_month[imo, ilat, ilon] 
    241                 diff_ice[imo, ilat, ilon] = diff_month[imo, ilat, ilon] 
    242                 ratio_ice[imo, ilat, ilon] = ratio_month[imo, ilat, ilon] 
     230    # calculate emis ratio daily 
     231    for ijr in range (0, month_day[imo]): 
     232        for ilon in range (0, nx): 
     233            for ilat in range (0, ny): 
     234                emis_ratio[imo, ilat, ilon, ijr] = ((emis_lamb[imo, ilat, ilon, ijr] - emis_spec[imo, ilat, ilon, ijr]) / emis_spec[imo, ilat, ilon, ijr]) * 100. 
     235                if (isnan(spec_lim[ilat, ilon]) == True): 
     236                    daily_spec_ice[imo, ilat, ilon, ijr] = NaN 
     237                    daily_lamb_ice[imo, ilat, ilon, ijr] = NaN 
     238                    daily_diff_ice[imo, ilat, ilon, ijr] = NaN 
     239                    daily_ratio_ice[imo, ilat, ilon, ijr] = NaN 
     240                else: 
     241                    daily_spec_ice[imo, ilat, ilon, ijr] = emis_spec[imo, ilat, ilon, ijr] 
     242                    daily_lamb_ice[imo, ilat, ilon, ijr] = emis_lamb[imo, ilat, ilon, ijr] 
     243                    daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr] 
     244                    daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr] 
     245    ''' 
     246    # ATTENTION : previous part of script has been modified ==> cannot be applied directly to this following part of script (modification of arrays and names.... 
    243247    print 'compute SPEC distribution' 
    244248    ######## 
    … …  
    282286    print 'start stacking in .dat file' 
    283287    #data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/AMSUB'+str(frequ)+'_data_classification_parameters_ice_no-ice_with_AMSUA23-spec_2009.dat', 'a') 
    284     '''for ii in range (0, bin): 
     288    for ii in range (0, bin): 
    285289        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_diff)10.5f    %(corresp_emis_diff)10.5f    %(hist_vals_rate)10.5f    %(corresp_emis_rate)10.5f    \n' % { 
    286290                                            'months':month[imo], 
    … …  
    293297                                            'hist_vals_rate':hist_vals_ratio[imo, ii], 
    294298                                            'corresp_emis_rate':corresp_emis_ratio[imo, ii], 
    295                                             }))''' 
     299                                            })) 
     300    ''' 
    296301    ######################## 
    297302    # stack in netcdf file # 
    298303    ########################  
    299304    print 'stack in file month ' + str(month[imo]) 
    300     rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-30_' + month[imo] + '2009_AMSUA' + str(frequ) + '_spec_lamb_thresholds.nc', 'w', format='NETCDF3_CLASSIC') 
    301     rootgrp.createDimension('longitude', len(xvec)) 
    302     rootgrp.createDimension('latitude', len(yvec)) 
     305    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUB' + str(frequ) + '_spec_thresholds.nc', 'w', format='NETCDF3_CLASSIC') 
     306    rootgrp.createDimension('longitude', nx) 
     307    rootgrp.createDimension('latitude', ny) 
     308    rootgrp.createDimension('days', month_day[imo]) 
    303309    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
    304310    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
    305     nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude')) 
    306     nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude')) 
    307     nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude')) 
    308     nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude')) 
     311    nc_days = rootgrp.createVariable('days', 'f', ('days',)) 
     312    nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     313    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     314    nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     315    nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude', 'days')) 
    309316    nc_lon[:] = xvec 
    310317    nc_lat[:] = yvec 
    311     nc_ice_spec[:] = spec_ice[imo, :, :] 
    312     nc_ice_lamb[:] = lamb_ice[imo, :, :] 
    313     nc_ice_diff[:] = diff_ice[imo, :, :] 
    314     nc_ice_ratio[:] = ratio_ice[imo, :, :] 
     318    nc_days[:] = np.arange(0, month_day[imo]) 
     319    nc_ice_spec[:] = daily_spec_ice[imo, :, :, 0 : month_day[imo]] 
     320    nc_ice_lamb[:] = daily_lamb_ice[imo, :, :, 0 : month_day[imo]] 
     321    nc_ice_diff[:] = daily_diff_ice[imo, :, :, 0 : month_day[imo]] 
     322    nc_ice_ratio[:] = daily_ratio_ice[imo, :, :, 0 : month_day[imo]] 
    315323    rootgrp.close() 
    316324 
    … …  
    320328 
    321329 
    322  
    323 ''' 
     330''' 
     331fichier = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUA' + str(frequ) + '_spec_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 
     332ice_spec = fichier.variables['spec_ice_area'][:] 
     333ice_lamb = fichier.variables['lamb_ice_area'][:] 
     334ice_ratio = fichier.variables['ratio_ice_area'][:] 
     335fichier.close() 
     336mean_ratio = np.zeros([ny, nx], float) 
     337for ilon in range (0, nx): 
     338    for ilat in range (0, ny): 
     339        mean_ratio[ilat, ilon] = mean(ice_ratio[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_ratio[ilat, ilon, 0 : month_day[imo]]) == False)]) 
     340 
     341 
     342ion() 
     343x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
     344x_coast = x_ind 
     345y_coast = y_ind 
     346z_coast = z_ind 
     347map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, mean_ratio[:, :], -3., 5., 0.1, cm.s3pcpn_l_r, 'test') 
     348 
     349 
     350 
     351 
    324352# test: 
    325353ion() 

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

    r48 r54  
    127127# plot time evolution of ice area # 
    128128################################### 
     129ion() 
    129130figure() 
    130131plot(AS[0, :], 'c', label = 'AMSUA spec_23GHz') 
    … …  
    138139plot(area_s_B[:], 'b', label = 'AMSUB spec_89GHz') 
    139140plot(area_l_B[:], 'b--', label = 'AMSUB lamb_89GHz') 
    140 plot(area_osi, 'k', label = 'OSISAF') 
     141plot(area_osi + 1500000, 'k', label = 'OSISAF + 1.5e7 (correction)') 
    141142fontP = FontProperties() 
    142143fontP.set_size('small') 
    … …  
    146147xlim(-1, M) 
    147148grid() 
    148 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/total_ice_area/compar_total_ice_area_AMSUA_AMSUB_SPEC_LAMB_OSISAF_2009.png') 
     149plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/total_ice_area/compar_total_ice_area_AMSUA_AMSUB_SPEC_LAMB_corrected_OSISAF_2009.png') 
    149150 
    150151 

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

    r47 r54  
    326326                else: 
    327327                    ice_zone_lamb[imo, ilat, ilon] = emis_lamb_month[imo, ilat, ilon]      
    328     map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_lamb[8, :, :], 0.5, 1.02, 0.01,  cm.s3pcpn_l_r, 'Sea ice emissivity lamb (threshold : emis_LAMB > ' + str(emis_lim_lamb[imo])[0:6] + ')') 
     328    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_lamb[imo, :, :], 0.5, 1.02, 0.01,  cm.s3pcpn_l_r, 'Sea ice emissivity lamb (threshold : emis_LAMB > ' + str(emis_lim_lamb[imo])[0:6] + ')') 
    329329    title(str(month[imo]) + ' 2009 - AMSUA ' + str(frequ) + 'GHz') 
    330330    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/AMSUA' + str(frequ) + '_ice_emis_lamb_thresh' + '_' + month[imo] + '2009.png') 

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

    r46 r54  
    5252    osi_type = fichier_osisaf.variables['osi_ice_type'][:] 
    5353    fichier_osisaf.close() 
    54     for ilon in range (0, len(xdist)): 
    55         for ilat in range (0, len(ydist)): 
     54    for ilon in range (0, nx): 
     55        for ilat in range (0, ny): 
    5656            for ijr in range (0, month_day[imo]): 
    5757                if ((xdist[ilon] >= -312.5) & (xdist[ilon] <= 328.5) & (ydist[ilat] >=-334.3) & (ydist[ilat] <= 338.5)): 
    5858                    osi_ice[imo, ijr, ilat, ilon] = 1. 
    5959                else: 
    60                     if (osi_type[ijr, ilat, ilon] >= 2.): 
     60                    if (osi_type[ijr, ilat, ilon] > 1.): 
    6161                        osi_ice[imo, ijr, ilat, ilon] = 1. 
    6262                    else: 
    … …  
    6464                            osi_ice[imo, ijr, ilat, ilon] = NaN 
    6565                        else:  
    66                             osi_ice[imo, ijr, ilat, ilon] = 0. 
     66                            osi_ice[imo, ijr, ilat, ilon] = NaN 
    6767 
    6868 
    6969# test: 
     70ion() 
    7071x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    7172x_coast = x_ind 
    … …  
    7374z_coast = z_ind 
    7475colormap = colors.ListedColormap(['0.5', 'cyan','blue', 'green']) 
    75 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_type[10, :, :], 0, 5, 1, colormap, 'ice_type') 
    76 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_ice[3, 7, :, :], -1., 2, 1, cm.s3pcpn_l_r, 'test') 
    77  
    78  
    79  
     76map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_type[1, :, :], 0, 5, 1, cm.s3pcpn_l_r, 'ice_type') 
     77map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_ice[0, 0, :, :], -1., 2, 1, cm.s3pcpn_l_r, 'test') 
    8078 
    8179 

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

    r49 r54  
    4747 
    4848frequ = np.array(['23', '30', '50', '89']) 
    49 ice_spec = np.zeros([4, M, ny, nx], float) 
    50 ice_lamb = np.zeros([4, M, ny, nx], float) 
    51 ice_diff = np.zeros([4, M, ny, nx], float) 
    52 ice_ratio = np.zeros([4, M, ny, nx], float) 
     49 
     50mean_ice_spec = np.zeros([4, M, ny, nx], float) 
     51mean_ice_lamb = np.zeros([4, M, ny, nx], float) 
     52mean_ice_diff = np.zeros([4, M, ny, nx], float) 
     53mean_ice_ratio = np.zeros([4, M, ny, nx], float) 
     54 
    5355ice_spec_anom = np.zeros([4, M, ny, nx], float) 
    5456ice_lamb_anom = np.zeros([4, M, ny, nx], float) 
    … …  
    5961    for imo in range (0, M): 
    6062        print 'month ' + month[imo] 
    61         fichier_ice = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-30_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 
    62         ice_spec[ifr, imo, :, :] = fichier_ice.variables['spec_ice_area'][:] 
    63         ice_lamb[ifr, imo, :, :] = fichier_ice.variables['lamb_ice_area'][:] 
    64         ice_diff[ifr, imo, :, :] = fichier_ice.variables['diff_ice_area'][:] 
    65         ice_ratio[ifr, imo, :, :] = fichier_ice.variables['ratio_ice_area'][:] 
     63        fichier_ice = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '_spec_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 
     64        ice_spec = fichier_ice.variables['spec_ice_area'][:] 
     65        ice_lamb = fichier_ice.variables['lamb_ice_area'][:] 
     66        ice_diff = fichier_ice.variables['diff_ice_area'][:] 
     67        ice_ratio = fichier_ice.variables['ratio_ice_area'][:] 
    6668        fichier_ice.close() 
    67         a = reshape(ice_spec[ifr, imo, :, :], size(ice_spec[ifr, imo, :, :])) 
    68         b = reshape(ice_lamb[ifr, imo, :, :], size(ice_lamb[ifr, imo, :, :])) 
    69         c = reshape(ice_diff[ifr, imo, :, :], size(ice_diff[ifr, imo, :, :])) 
    70         d = reshape(ice_ratio[ifr, imo, :, :], size(ice_ratio[ifr, imo, :, :])) 
     69        for ilat in range (0, ny): 
     70            for ilon in range (0, nx): 
     71                mean_ice_spec[ifr, imo, ilat, ilon] = mean(ice_spec[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_spec[ilat, ilon, 0 : month_day[imo]]) == False)]) 
     72                mean_ice_lamb[ifr, imo, ilat, ilon] = mean(ice_lamb[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_lamb[ilat, ilon, 0 : month_day[imo]]) == False)]) 
     73                mean_ice_diff[ifr, imo, ilat, ilon] = mean(ice_diff[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_diff[ilat, ilon, 0 : month_day[imo]]) == False)]) 
     74                mean_ice_ratio[ifr, imo, ilat, ilon] = mean(ice_ratio[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_ratio[ilat, ilon, 0 : month_day[imo]]) == False)]) 
     75        a = reshape(mean_ice_spec[ifr, imo, :, :], size(mean_ice_spec[ifr, imo, :, :])) 
     76        b = reshape(mean_ice_lamb[ifr, imo, :, :], size(mean_ice_lamb[ifr, imo, :, :])) 
     77        c = reshape(mean_ice_diff[ifr, imo, :, :], size(mean_ice_diff[ifr, imo, :, :])) 
     78        d = reshape(mean_ice_ratio[ifr, imo, :, :], size(mean_ice_ratio[ifr, imo, :, :])) 
    7179        print 'calculate anomaly' 
    7280        for ilat in range (0, ny): 
    7381            for ilon in range (0, nx): 
    74                 ice_spec_anom[ifr, imo, ilat, ilon] = ice_spec[ifr, imo, ilat, ilon] - mean(a[nonzero(isnan(a) == False)]) 
    75                 ice_lamb_anom[ifr, imo, ilat, ilon] = ice_lamb[ifr, imo, ilat, ilon] - mean(b[nonzero(isnan(b) == False)]) 
    76                 ice_diff_anom[ifr, imo, ilat, ilon] = ice_diff[ifr, imo, ilat, ilon] - mean(c[nonzero(isnan(c) == False)]) 
    77                 ice_ratio_anom[ifr, imo, ilat, ilon] = ice_ratio[ifr, imo, ilat, ilon] - mean(d[nonzero(isnan(d) == False)]) 
     82                ice_spec_anom[ifr, imo, ilat, ilon] = mean_ice_spec[ifr, imo, ilat, ilon] - mean(a[nonzero(isnan(a) == False)]) 
     83                ice_lamb_anom[ifr, imo, ilat, ilon] = mean_ice_lamb[ifr, imo, ilat, ilon] - mean(b[nonzero(isnan(b) == False)]) 
     84                ice_diff_anom[ifr, imo, ilat, ilon] = mean_ice_diff[ifr, imo, ilat, ilon] - mean(c[nonzero(isnan(c) == False)]) 
     85                ice_ratio_anom[ifr, imo, ilat, ilon] = mean_ice_ratio[ifr, imo, ilat, ilon] - mean(d[nonzero(isnan(d) == False)]) 
    7886        ######################## 
    7987        # stack in netcdf file # 
    8088        ########################  
    8189        print 'stack in file month ' + str(month[imo]) 
    82         rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-30_param_anomaly_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '.nc', 'w', format='NETCDF3_CLASSIC') 
    83         rootgrp.createDimension('longitude', len(xvec)) 
    84         rootgrp.createDimension('latitude', len(yvec)) 
     90        rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_param_anomaly_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '.nc', 'w', format='NETCDF3_CLASSIC') 
     91        rootgrp.createDimension('longitude', nx) 
     92        rootgrp.createDimension('latitude', ny) 
    8593        nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
    8694        nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
    … …  
    109117for ifr in range (0, 2): 
    110118    for imo in range (0, M): 
    111         map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ice_diff_anom[ifr, imo, :, :], ice_diff_anom[nonzero(isnan(ice_diff_anom) == False)].min(), ice_diff_anom[nonzero(isnan(ice_diff_anom) == False)].max(), 0.0001, cm.s3pcpn_l_r, 'Emissivity diff anomaly') 
     119        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ice_ratio_anom[ifr, imo, :, :], ice_ratio_anom[nonzero(isnan(ice_ratio_anom) == False)].min(), ice_ratio_anom[nonzero(isnan(ice_ratio_anom) == False)].max(), 0.0001, cm.s3pcpn_l_r, 'Emissivity diff anomaly') 
    112120        title('AMSUA ' + str(frequ) + ' - ' + str(month[imo]) + ' 2009') 
    113121        plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/emiss_diff_anomaly_map_AMSUA'+str(frequ[ifr])+'_with_AMSUA23-and-30-spec_'+str(month[imo])+'_2009.png') 

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

    r42 r54  
    99from netCDF4 import Dataset 
    1010import scipy.special 
    11 import ffgrid2 
    12 import map_ffgrid 
    1311from matplotlib import colors 
    1412import cartesian_grid_monthly 
    … …  
    3028# grid parameters # 
    3129################### 
    32 #x0 = -3000. # min limit of grid 
    33 #x1 = 2500. # max limit of grid 
    34 #dx = 40. 
    35 #xvec = np.arange(x0, x1+dx, dx) 
    36 #nx = len(xvec)  
    37 #y0 = -3000. # min limit of grid 
    38 #y1 = 3000. # max limit of grid 
    39 #dy = 40. 
    40 #yvec = np.arange(y0, y1+dy, dy) 
    41 #ny = len(yvec) 
     30x0 = -3000. # min limit of grid 
     31x1 = 2500. # max limit of grid 
     32dx = 12.5 
     33xvec = np.arange(x0, x1+dx, dx) 
     34nx = len(xvec)  
     35y0 = -3000. # min limit of grid 
     36y1 = 3000. # max limit of grid 
     37dy = 12.5 
     38yvec = np.arange(y0, y1+dy, dy) 
     39ny = len(yvec) 
    4240 
    4341 
    … …  
    5351        print 'date : ' + str(DAY[ijr]) + '_' + str(month[imo]) + '_2009' 
    5452# Concentration  
    55 ifremer_conc = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/CERSAT-IFREMER/sea_ice_concentration/2009' + MONTH[imo] + DAY[ijr] + '.nc', 'r', format='NETCDF3_CLASSIC') 
    56 conc = reshape(ifremer_conc.variables['concentration'][:], 608*896) 
     53ifremer_conc = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/CERSAT-IFREMER/sea_ice_concentration/20090101.nc', 'r', format='NETCDF3_CLASSIC') 
     54conc = ifremer_conc.variables['concentration'][:] 
    5755ifremer_conc.close() 
    5856# Backscatter 
    59 ifremer_back = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/CERSAT-IFREMER/sea_ice_backscatter/ASCAT_2009' + MONTH[imo] + DAY[ijr] + '.nc', 'r', format='NETCDF3_CLASSIC') 
    60 lon = reshape(ifremer_back.variables['longitude'][:], 608*896) 
    61 lat = reshape(ifremer_back.variables['latitude'][:], 608*896) 
    62 back = reshape(ifremer_back.variables['sigma_40'][:], 608*896) 
     57ifremer_back = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/CERSAT-IFREMER/sea_ice_backscatter/ASCAT_20090101.nc', 'r', format='NETCDF3_CLASSIC') 
     58lon = ifremer_back.variables['longitude'][:] 
     59lat = ifremer_back.variables['latitude'][:] 
     60back = ifremer_back.variables['sigma_40'][:] 
    6361ifremer_back.close() 
    6462bb = nonzero((lon >= -180.) & (lon <= 180.)) 
    … …  
    6664llat = lat[bb] 
    6765bback = back[bb] 
     66 
     67 
     68ion() 
     69figure() 
     70m = Basemap(width=15.e6,height=15.e6,\ 
     71            projection='gnom',lat_0=60.,lon_0=-30.) 
     72m.drawcoastlines(linewidth=1) 
     73cs=m.contourf(lon[0, :, :], lat[0, :, :], back[0, :, :], clevs, cmap = cm.s3pcpn_l_r) 
     74 
     75m = Basemap(llcrnrlon=-180,urcrnrlon=180,llcrnrlat=30,urcrnrlat=90,projection='cyl',resolution='c',fix_aspect=True) 
     76xii, yii = m(*np.meshgrid(x,y)) 
     77clevs=arange(-0.06,0.06,0.001) 
     78 
     79emissivity 
     80cbar =colorbar(cs) 
     81cbar.set_label('test', color='k') 
     82 
     83         
     84 
     85 
    6886 
    6987contourf(llon, llat, bback, cm.s3pcpn_l_r) 
Note: See TracChangeset for help on using the changeset viewer.