Changeset 56


Ignore:
Timestamp:
08/22/14 15:06:02 (10 years ago)
Author:
lahlod
Message:

other_scripts

Location:
trunk/src/scripts_Laura/ARCTIC/Travail_CEN
Files:
5 added
4 deleted
20 edited

Legend:

Unmodified
Added
Removed

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

    r41 r56  
    1111 
    1212 
     13####################################################### 
     14# Defines every point from bathymetry NETCDF file     # 
     15# that corresponds to altitude = 0 (coast) and stacks # 
     16# these points in cartesian coordinate vectors x_ind  # 
     17# and y_ind                                           # 
     18####################################################### 
    1319 
    1420def arctic_map_lat50(): 
    … …  
    4046    z_ind = Z_LIM[indices] 
    4147 
    42     volume = z_ind / 50. 
     48    volume = z_ind / 50. # size of the points defining the coast (with z_ind = 1.) 
    4349    #plt.ion() 
    4450    #fig = plt.figure() 

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

    r46 r56  
    3232    # definition of the new cartesian grid (x, y) # 
    3333    ###############################################  
    34     x0 = -3000. # min limit of grid 
    35     x1 = 2500. # max limit of grid 
     34    x0 = -3000. # min x limit of grid (in km) 
     35    x1 = 2500. # max x limit of grid (in km) 
    3636    xvec = np.arange(x0, x1+dx, dx) 
    3737    nx = len(xvec)  
    38     y0 = -3000. # min limit of grid 
    39     y1 = 3000. # max limit of grid 
     38    y0 = -3000. # min y limit of grid (in km) 
     39    y1 = 3000. # max y limit of grid (in km) 
    4040    yvec = np.arange(y0, y1+dy, dy) 
    4141    ny = len(yvec) 
    42     xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole) 
     42    xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole [x, y] = [0, 0]) 
    4343 
    4444 
    … …  
    4646    # calculation of the new gridding # 
    4747    ###################################  
    48     zgrid_output = np.zeros([ny, nx, nb_days], float) 
    49     ngrid_output = np.zeros([ny, nx, nb_days], float) 
    50     z2grid_output = np.zeros([ny, nx, nb_days], float) 
    51     sigmagrid_output = np.zeros([ny, nx, nb_days], float) 
    52     bblat = nonzero(lati >= 65.) # only select high latitude values of z 
     48    zgrid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of data in new cartesian grid 
     49    ngrid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of number of data in each grid cell of new cartesian grid 
     50    z2grid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of data*data in new cartesian grid 
     51    sigmagrid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of std of data in each grid cell of new cartesian grid 
     52    bblat = nonzero(lati >= 65.) # only select high latitude observations (above 65° lat) 
    5353    for ijr in range (0, nb_days): # loop on time - here time = days 
    5454        bbjour = nonzero(jour[bblat] == ijr + 1.)[0] 
    … …  
    6161        # origin of the (x, y) grid : (x=0, y=0) <=> (lon=0, lat=0)  
    6262        L = len(z) 
    63         Rt = 6371. 
    64         alpha = (pi*Rt)/180. 
     63        Rt = 6371. # radius of earth 
     64        alpha = (pi*Rt)/180. # convertion from deg to rad angles 
    6565        beta = pi/180. 
    66         x = np.zeros([L], float) 
    67         y = np.zeros([L], float) 
     66        x = np.zeros([L], float) # x coordinates of new cartesian grid 
     67        y = np.zeros([L], float) # y coordinates of new cartesian grid 
    6868        for k in range (0, L): 
    6969                if ((longitude[k] >= 0.) & (longitude[k] <= 90.)): # 4eme quadrant 
    … …  
    9090                        ix[i] = 0 
    9191                else: 
    92                         ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x vakue to a cell number 
     92                        ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x value to a cell number 
    9393        i = 0 
    9494        iy = np.zeros([L], int)  
    … …  
    101101        # calculation of distances between point M(x,y) and 4 grid points of its belonging cell # 
    102102        ######################################################################################### 
    103         close_point = np.zeros([L, 2], int) 
     103        close_point = np.zeros([L, 2], int) # for x- and y-axis, 2D-array of closest grid cell of new cartesian grid to observation  
    104104        k = 0 
    105105        for k in range (0, L): # boucle sur les points M (x et y) 
    … …  
    108108                d3 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2)) 
    109109                d4 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2)) 
    110                 d_vec = np.array([d1, d2, d3, d4]) 
     110                d_vec = np.array([d1, d2, d3, d4]) # vector of distances between observation and each closer grid_cell of new cartesian grid 
    111111                ind = np.where(d_vec == min(d_vec)) # finds the smallest distance between the 4 points of the grid 
    112112                point_vec = np.array([(ix[k], iy[k]), (ix[k] + 1, iy[k]), (ix[k], iy[k] + 1), (ix[k] + 1, iy[k] + 1)])  

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

    r54 r56  
    4747# - second loop concerns MAY, JUNE, JULY, AUGUST - use of AMSUA89GHz SPEC emissivity to seperate ice from no_ice zones 
    4848################################################################################################################## 
     49 
     50 
    4951frequ = 89 # apply threshold on this frequency 
    50 ''' 
    51 #open .dat file to stack data (see end of loop) 
    52 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') 
    53 bin = 50 
    54 ''' 
    5552 
    5653 
    … …  
    6764daily_ratio_ice = np.zeros([M, ny, nx, 31], float) 
    6865 
    69 ''' 
    70 # monthly mean parameter (1D-array) on ARCTIC SEA ICE area transformed into vector 
    71 spec_vec = np.zeros([M, ny * nx], float) 
    72 lamb_vec = np.zeros([M, ny * nx], float) 
    73 diff_vec = np.zeros([M, ny * nx], float) 
    74 ratio_vec = np.zeros([M, ny * nx], float) 
    75  
    76 # histogram distribution (intensity of occurence) of parameter in SEA ICE area (1D-array, bins = 200) 
    77 hist_vals_spec = np.zeros([M, bin], float) 
    78 hist_vals_lamb = np.zeros([M, bin], float) 
    79 hist_vals_diff = np.zeros([M, bin], float) 
    80 hist_vals_ratio = np.zeros([M, bin], float) 
    81  
    82 # histogram distribution (emissivity value) of parameter in SEA ICE area (1D-array, bins = 200) 
    83 corresp_emis_spec = np.zeros([M, bin], float) 
    84 corresp_emis_lamb = np.zeros([M, bin], float) 
    85 corresp_emis_diff = np.zeros([M, bin], float) 
    86 corresp_emis_ratio = np.zeros([M, bin], float) 
    87 ''' 
     66 
    8867months1 = 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 
    8968for imo in months1: 
    … …  
    9473    print 'open threshold file' 
    9574    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') 
     75    spec_lim = fichier_emis.variables['spec_ice_area'][:] # sea ice pixels defined with spec emis at 23GHz  
     76    #lamb_lim = fichier_emis.variables['lamb_ice_area'][:] 
     77    fichier_emis.close() 
     78    ######################################################### 
     79    # application of AMSUA 23GHz delimitation to other data # 
     80    ######################################################### 
     81    print 'open emissivity to sub_classify file' 
     82    ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB) 
     83    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') 
     84    day = fichier_e.variables['days'][:] 
     85    emis_spec[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec'][:] 
     86    emis_lamb[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_lamb'][:] 
     87    emis_diff[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec_lamb'][:] 
     88    fichier_e.close() 
     89    # calculate emis ratio daily 
     90    for ijr in range (0, month_day[imo]): 
     91        for ilon in range (0, nx): 
     92            for ilat in range (0, ny): 
     93                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. 
     94                # create 2D-array of emissivity spec, lamb, diff and ratio on sea ice extent only, defined by AMSUA 23GHz spec emiss threshold    
     95                if (isnan(spec_lim[ilat, ilon]) == True): # if pixel of sea ice extent defined with spec_emiss_23_threshold corresponds to 'no_ice', then compute NaN in pixel 
     96                    daily_spec_ice[imo, ilat, ilon, ijr] = NaN 
     97                    daily_lamb_ice[imo, ilat, ilon, ijr] = NaN 
     98                    daily_diff_ice[imo, ilat, ilon, ijr] = NaN 
     99                    daily_ratio_ice[imo, ilat, ilon, ijr] = NaN 
     100                else: # if pixel of sea ice extent defined with spec_emiss_23_threshold corresponds to 'ice', then compute value of emis spec, emis lamb or emis diff in pixel 
     101                    daily_spec_ice[imo, ilat, ilon, ijr] = emis_spec[imo, ilat, ilon, ijr] 
     102                    daily_lamb_ice[imo, ilat, ilon, ijr] = emis_lamb[imo, ilat, ilon, ijr] 
     103                    daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr] 
     104                    daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr] 
     105    ######################## 
     106    # stack in netcdf file # 
     107    ########################  
     108    print 'stack in file month ' + str(month[imo]) 
     109    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') 
     110    rootgrp.createDimension('longitude', nx) 
     111    rootgrp.createDimension('latitude', ny) 
     112    rootgrp.createDimension('days', month_day[imo]) 
     113    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
     114    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
     115    nc_days = rootgrp.createVariable('days', 'f', ('days',)) 
     116    nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     117    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     118    nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     119    nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude', 'days')) 
     120    nc_lon[:] = xvec 
     121    nc_lat[:] = yvec 
     122    nc_days[:] = np.arange(0, month_day[imo]) 
     123    nc_ice_spec[:] = daily_spec_ice[imo, :, :, 0 : month_day[imo]] 
     124    nc_ice_lamb[:] = daily_lamb_ice[imo, :, :, 0 : month_day[imo]] 
     125    nc_ice_diff[:] = daily_diff_ice[imo, :, :, 0 : month_day[imo]] 
     126    nc_ice_ratio[:] = daily_ratio_ice[imo, :, :, 0 : month_day[imo]] 
     127    rootgrp.close() 
     128 
     129 
     130 
     131 
     132months2 = np.array([4, 5, 6, 7])# use AMSUA 89GHz to delimit ice/no_ice for MAY, JUNE, JULY, AUGUST 
     133for imo in months2: 
     134    print 'month ' + month[imo] 
     135    ################################################################################## 
     136    # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study # 
     137    ################################################################################## 
     138    print 'open threshold file' 
     139    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') 
    96140    spec_lim = fichier_emis.variables['spec_ice_area'][:] 
    97141    #lamb_lim = fichier_emis.variables['lamb_ice_area'][:] 
    … …  
    123167                    daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr] 
    124168                    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.... 
    127     print 'compute SPEC distribution' 
    128     ######## 
    129     # SPEC # 
    130     ######## 
    131     cs = reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))[nonzero(isnan(reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))) == False)] 
    132     spec_vec[imo, 0 : len(cs)] = cs 
    133     hist_vals_spec[imo, :] = hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[0] 
    134     for ibin in range (0, bin): 
    135         corresp_emis_spec[imo, ibin] = mean(hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    136     print 'compute LAMB distribution' 
    137     ######## 
    138     # LAMB # 
    139     ######## 
    140     cl = reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))[nonzero(isnan(reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))) == False)] 
    141     lamb_vec[imo, 0 : len(cl)] = cl 
    142     hist_vals_lamb[imo, :] = hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[0] 
    143     for ibin in range (0, bin): 
    144         corresp_emis_lamb[imo, ibin] = mean(hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    145     print 'compute DIFF distribution' 
    146     ######## 
    147     # DIFF # 
    148     ######## 
    149     cd = reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))[nonzero(isnan(reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))) == False)] 
    150     diff_vec[imo, 0 : len(cd)] = cd 
    151     hist_vals_diff[imo, :] = hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[0] 
    152     for ibin in range (0, bin): 
    153         corresp_emis_diff[imo, ibin] = mean(hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    154     print 'compute RATIO distribution' 
    155     ######### 
    156     # RATIO # 
    157     ######### 
    158     cr = reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))[nonzero(isnan(reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))) == False)] 
    159     ratio_vec[imo, 0 : len(cr)] = cr 
    160     hist_vals_ratio[imo, :] = hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[0] 
    161     for ibin in range (0, bin): 
    162         corresp_emis_ratio[imo, ibin] = mean(hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    163     ###################### 
    164     # stack in .dat file # 
    165     ###################### 
    166     print 'start stacking in .dat file' 
    167     #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') 
    168     for ii in range (0, bin): 
    169         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' % { 
    170                                             'months':month[imo], 
    171                                             'hist_vals_spec':hist_vals_spec[imo, ii], 
    172                                             'corresp_emis_spec':corresp_emis_spec[imo, ii], 
    173                                             'hist_vals_lamb':hist_vals_lamb[imo, ii], 
    174                                             'corresp_emis_lamb':corresp_emis_lamb[imo, ii], 
    175                                             'hist_vals_diff':hist_vals_diff[imo, ii], 
    176                                             'corresp_emis_diff':corresp_emis_diff[imo, ii], 
    177                                             'hist_vals_rate':hist_vals_ratio[imo, ii], 
    178                                             'corresp_emis_rate':corresp_emis_ratio[imo, ii], 
    179                                             })) 
    180     ''' 
    181169    ######################## 
    182170    # stack in netcdf file # 
    … …  
    206194 
    207195 
    208 months2 = np.array([4, 5, 6, 7])# use AMSUA 89GHz to delimit ice/no_ice for MAY, JUNE, JULY, AUGUST 
    209 for imo in months2: 
    210     print 'month ' + month[imo] 
    211     ################################################################################## 
    212     # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study # 
    213     ################################################################################## 
    214     print 'open threshold file' 
    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') 
    216     spec_lim = fichier_emis.variables['spec_ice_area'][:] 
    217     #lamb_lim = fichier_emis.variables['lamb_ice_area'][:] 
    218     fichier_emis.close() 
    219     ######################################################### 
    220     # application of AMSUA 23GHz delimitation to other data # 
    221     ######################################################### 
    222     print 'open emissivity to sub_classify file' 
    223     ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB) 
    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') 
    225     day = fichier_e.variables['days'][:] 
    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'][:] 
    229     fichier_e.close() 
    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.... 
    247     print 'compute SPEC distribution' 
    248     ######## 
    249     # SPEC # 
    250     ######## 
    251     cs = reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))[nonzero(isnan(reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))) == False)] 
    252     spec_vec[imo, 0 : len(cs)] = cs 
    253     hist_vals_spec[imo, :] = hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[0] 
    254     for ibin in range (0, bin): 
    255         corresp_emis_spec[imo, ibin] = mean(hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    256     print 'compute LAMB distribution' 
    257     ######## 
    258     # LAMB # 
    259     ######## 
    260     cl = reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))[nonzero(isnan(reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))) == False)] 
    261     lamb_vec[imo, 0 : len(cl)] = cl 
    262     hist_vals_lamb[imo, :] = hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[0] 
    263     for ibin in range (0, bin): 
    264         corresp_emis_lamb[imo, ibin] = mean(hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    265     print 'compute DIFF distribution' 
    266     ######## 
    267     # DIFF # 
    268     ######## 
    269     cd = reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))[nonzero(isnan(reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))) == False)] 
    270     diff_vec[imo, 0 : len(cd)] = cd 
    271     hist_vals_diff[imo, :] = hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[0] 
    272     for ibin in range (0, bin): 
    273         corresp_emis_diff[imo, ibin] = mean(hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    274     print 'compute RATIO distribution' 
    275     ######### 
    276     # RATIO # 
    277     ######### 
    278     cr = reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))[nonzero(isnan(reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))) == False)] 
    279     ratio_vec[imo, 0 : len(cr)] = cr 
    280     hist_vals_ratio[imo, :] = hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[0] 
    281     for ibin in range (0, bin): 
    282         corresp_emis_ratio[imo, ibin] = mean(hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    283     ###################### 
    284     # stack in .dat file # 
    285     ###################### 
    286     print 'start stacking in .dat file' 
    287     #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') 
    288     for ii in range (0, bin): 
    289         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' % { 
    290                                             'months':month[imo], 
    291                                             'hist_vals_spec':hist_vals_spec[imo, ii], 
    292                                             'corresp_emis_spec':corresp_emis_spec[imo, ii], 
    293                                             'hist_vals_lamb':hist_vals_lamb[imo, ii], 
    294                                             'corresp_emis_lamb':corresp_emis_lamb[imo, ii], 
    295                                             'hist_vals_diff':hist_vals_diff[imo, ii], 
    296                                             'corresp_emis_diff':corresp_emis_diff[imo, ii], 
    297                                             'hist_vals_rate':hist_vals_ratio[imo, ii], 
    298                                             'corresp_emis_rate':corresp_emis_ratio[imo, ii], 
    299                                             })) 
    300     ''' 
    301     ######################## 
    302     # stack in netcdf file # 
    303     ########################  
    304     print 'stack in file month ' + str(month[imo]) 
    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]) 
    309     nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 
    310     nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 
    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')) 
    316     nc_lon[:] = xvec 
    317     nc_lat[:] = yvec 
    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]] 
    323     rootgrp.close() 
    324196 
    325197''' 
    326 data_classif.close() 
    327 ''' 
    328  
    329  
    330 ''' 
    331 fichier = 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') 
    332 ice_spec = fichier.variables['spec_ice_area'][:] 
    333 ice_lamb = fichier.variables['lamb_ice_area'][:] 
    334 ice_ratio = fichier.variables['ratio_ice_area'][:] 
    335 fichier.close() 
    336 mean_ratio = np.zeros([ny, nx], float) 
    337 for 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  
    342 ion() 
    343 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    344 x_coast = x_ind 
    345 y_coast = y_ind 
    346 z_coast = z_ind 
    347 map_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  
    352198# test: 
    353199ion() 

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

    r54 r56  
    1919 
    2020 
    21  
     21######################## 
     22# Time characteristics # 
     23######################## 
    2224MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) 
    2325month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER']) 
    … …  
    2527M = len(month) 
    2628 
    27 frequ = np.array([23, 30, 50, 89]) 
     29 
     30frequ = np.array([23, 30, 50, 89]) # frequencies of AMSUA used for the study 
    2831 
    2932 
    … …  
    4851         liste = line.split() 
    4952         # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es) 
    50          tot_area_spec[iligne] = float(liste[9]) 
    51          tot_area_lamb[iligne] = float(liste[10]) 
     53         tot_area_spec[iligne] = float(liste[9]) # total sea ice area defined with emis_spec threshold 
     54         tot_area_lamb[iligne] = float(liste[10]) # total sea ice area defined with emis_lamb threshold 
    5255         iligne=iligne+1 
    5356    fichier.close() 
    5457    vec = np.arange(0, nbtotal + 51, 50) 
    55     area_s = np.zeros([M], float) 
    56     area_l = np.zeros([M], float) 
     58    area_s = np.zeros([M], float) # vector of monthly mean total ice area defined with emis_spec threshold  
     59    area_l = np.zeros([M], float) # vector of monthly mean total ice area defined with emis_spec threshold  
    5760    for imo in range (0, M): 
    5861        area_s[imo] = tot_area_spec[imo + vec[imo]] 
    5962        area_l[imo] = tot_area_lamb[imo + vec[imo]] 
    60     AS[ifr, :] = area_s[:] 
    61     AL[ifr, :] = area_l[:] 
     63    AS[ifr, :] = area_s[:] # create 2D-array of monthly mean total sea ice area defined with emis_spec threshold, for each frequency 
     64    AL[ifr, :] = area_l[:] # create 2D-array of monthly mean total sea ice area defined with emis_lamb threshold, for each frequency 
    6265 
    6366 
    … …  
    127130# plot time evolution of ice area # 
    128131################################### 
     132# plot various total sea ice areas for each data, monthly in 2009 
    129133ion() 
    130134figure() 
    … …  
    139143plot(area_s_B[:], 'b', label = 'AMSUB spec_89GHz') 
    140144plot(area_l_B[:], 'b--', label = 'AMSUB lamb_89GHz') 
    141 plot(area_osi + 1500000, 'k', label = 'OSISAF + 1.5e7 (correction)') 
     145plot(area_osi + 1500000, 'k', label = 'OSISAF + 1.5e7 (correction)') # '+1500000' corresponds to translation of OSISAF curve due to constant bias between OSISAF and AMSU total sea ice areas - quantity of translation is arbitrary 
    142146fontP = FontProperties() 
    143147fontP.set_size('small') 
    … …  
    150154 
    151155 
    152  
    153 ############################### 
    154 # calculation of bias and std # 
    155 ############################### 
     156''' 
     157################################################################################## 
     158# calculation of bias between each AMSU total ice area and OSISAF total ice area # 
     159################################################################################## 
    156160a = np.zeros([M], float) 
    157161b = np.zeros([M], float) 
    … …  
    197201grid() 
    198202plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/total_ice_area/bias_total_ice_area_AMSU_OSI_2009.png') 
    199  
    200  
    201  
    202  
    203  
     203''' 
     204 
     205 
     206 
     207 

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

    r55 r56  
    3030# grid characteristics # 
    3131######################## 
     32### X AXIS 
    3233x0 = -3000. # min limit of grid 
    3334x1 = 2500. # max limit of grid 
    … …  
    4142nx_b = len(xvec_b)  
    4243 
     44### Y AXIS 
    4345y0 = -3000. # min limit of grid 
    4446y1 = 3000. # max limit of grid 
    … …  
    8587 
    8688 
    87  
     89''' 
     90##################### 
     91# map area of study # 
     92##################### 
    8893ion() 
    8994x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    … …  
    9297z_coast = z_ind 
    9398map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec_b[xxi_b:xxf_b+1], yvec_b[yyi_b:yyf_b+1], spec_b[yyi_b:yyf_b+1, xxi_b:xxf_b+1, 0], 0.6, 1.02, 0.01, cm.s3pcpn_l_r, 'test') 
    94  
     99''' 
    95100 
    96101 
    … …  
    171176 
    172177 
     178 
     179 
     180######################################################### 
     181# Join daily evolution of each month in a yearly vector # 
     182######################################################### 
    173183mean_year_spec_a = np.append(mean_spec_a_zone[0, 0 : month_day[0]], mean_spec_a_zone[1, 0 : month_day[1]]) 
    174184mean_year_spec_a23 = np.append(mean_spec_a_zone23[0, 0 : month_day[0]], mean_spec_a_zone23[1, 0 : month_day[1]]) 
    … …  
    204214# plot daily parameters # 
    205215######################### 
     216### plot difference between the two zonal std (of amsua and amsub at 89GHz) 
    206217vec_months = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) 
    207218ion() 
    … …  
    222233 
    223234 
    224  
     235### plot zonal mean of emissivity spec of AMSUA and B at 89GHz and the difference between them (subplot) 
    225236vec_months = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) 
    226237ion() 
    … …  
    264275#title('Chukchi Sea') 
    265276plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/study_by_zones/monthly_evolution_emis_AMSUA89_AMSUB89_AMSUA23_bias_std_mean_zone_Chukchi_Sea.png') 
    266 ''' 
    267  
    268  
    269  
    270  
    271  
    272  
    273  
    274 ######################################## 
    275 # monthly mean and std in whole arctic #  
    276 ######################################## 
    277 std_spec_a89 = np.zeros([M, ny_a, nx_a], float) 
    278 std_spec_a23 = np.zeros([M, ny_a, nx_a], float) 
    279 std_spec_b = np.zeros([M, ny_a, nx_a], float) 
    280 std_lamb_a89 = np.zeros([M, ny_a, nx_a], float) 
    281 std_lamb_a23 = np.zeros([M, ny_a, nx_a], float) 
    282 std_lamb_b = np.zeros([M, ny_a, nx_a], float) 
    283 '''bias_spec = np.zeros([M, ny_a, nx_a], float) 
    284 bias_anom = np.zeros([M, ny_a, nx_a], float) 
    285 bias_mean = np.zeros([M], float)''' 
    286 for imo in range (0, M): 
    287     print 'month ' + month[imo] 
    288     print 'read daily emis spec lamb and diff AMSUA and AMSUB' 
    289     # AMSUA 89 
    290     fichier_a89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') 
    291     spec_a89 = fichier_a89.variables['mean_spec'][:] 
    292     lamb_a89 = fichier_a89.variables['mean_lamb'][:] 
    293     std_spec_a89[imo, :, :] = fichier_a89.variables['std_spec'][:] 
    294     std_lamb_a89[imo, :, :] = fichier_a89.variables['std_lamb'][:] 
    295     fichier_a89.close() 
    296     # AMSUA 23 
    297     fichier_b89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') 
    298     spec_a23 = fichier_a23.variables['mean_spec'][:] 
    299     lamb_a23 = fichier_a23.variables['mean_lamb'][:] 
    300     std_spec_a23[imo, :, :] = fichier_a23.variables['std_spec'][:] 
    301     std_lamb_a23[imo, :, :] = fichier_a23.variables['std_lamb'][:] 
    302     fichier_a23.close() 
    303     # AMSUB 
    304     fichier_b = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') 
    305     spec_b = fichier_b.variables['mean_spec'][:] 
    306     lamb_b = fichier_b.variables['mean_lamb'][:] 
    307     std_spec_b[imo, :, :] = fichier_b.variables['std_spec'][:] 
    308     std_lamb_b[imo, :, :] = fichier_b.variables['std_lamb'][:] 
    309     fichier_b.close() 
    310     '''bias_spec[imo, :, :] = spec_a - spec_b 
    311     bias_mean[imo] = mean(bias_spec[imo, :, :][nonzero(isnan(bias_spec[imo, :, :]) == False)]) 
    312     for ilon in range (0, nx_a): 
    313         for ilat in range (0, ny_a): 
    314             bias_anom[imo, ilat, ilon] = bias_spec[imo, ilat, ilon] - bias_mean[imo]''' 
    315     
    316 # calculate difference of std between spec and lamb for amsua and amsub 
    317 c = np.zeros([M], float) 
    318 f = np.zeros([M], float) 
    319 for imo in range (0, M): 
    320    a = mean(std_spec_a89[imo, :, :][nonzero(isnan(std_spec_a89[imo, :, :])==False)])  
    321    b = mean(std_lamb_a89[imo, :, :][nonzero(isnan(std_lamb_a89[imo, :, :])==False)])  
    322    c[imo] = a - b 
    323    d = mean(std_spec_b[imo, :, :][nonzero(isnan(std_spec_b[imo, :, :])==False)])  
    324    e = mean(std_lamb_b[imo, :, :][nonzero(isnan(std_lamb_b[imo, :, :])==False)]) 
    325    f[imo] = d - e 
    326  
    327 # calculate difference of std between spec amsua and spec amsub 
    328 k = np.zeros([M], float) 
    329 for imo in range (0, M): 
    330    g = mean(std_spec_a89[imo, :, :][nonzero(isnan(std_spec_a89[imo, :, :])==False)])  
    331    h = mean(std_spec_b[imo, :, :][nonzero(isnan(std_spec_b[imo, :, :])==False)])  
    332    k[imo] = abs(g - h) 
    333  
    334  
    335  
    336 ################################## 
    337 # map bias spec_AMSUA-spec_AMSUB # 
    338 ################################## 
    339 #ion() 
    340 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    341 x_coast = x_ind 
    342 y_coast = y_ind 
    343 z_coast = z_ind 
    344 for imo in range (0, M): 
    345     print 'map for month ' + month[imo] 
    346     map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec_a, yvec_a, std_lamb_a89[imo, :, :], 0., 0.12, 0.001, cm.s3pcpn_l_r, 'std emis lamb AMSUA89') 
    347     title(month[imo] + ' 2009') 
    348     plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/monthly_std/std_emis_lamb_AMSUA89_' + month[imo] + '2009.png') 
    349  

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

    r52 r56  
    4343 
    4444 
    45 ''' 
     45 
    4646############################################ 
    4747# time evolution (monthly) in a given zone # 
    … …  
    6666len(yvec[yyi:yyf+1]) 
    6767 
    68 mean_grad_ratio_zone = np.zeros([M, 31], float) 
    69 std_grad_ratio_zone = np.zeros([M, 31], float) 
    70  
    71 mean_spec_anom_23_zone = np.zeros([M, 31], float) 
    72 std_spec_anom_23_zone = np.zeros([M, 31], float) 
    73  
    74 mean_spec_anom_89_zone = np.zeros([M, 31], float) 
    75 std_spec_anom_89_zone = np.zeros([M, 31], float) 
    76  
    77 mean_ratio_anom_89_zone = np.zeros([M, 31], float) 
    78 std_ratio_anom_89_zone = np.zeros([M, 31], float) 
    79  
    80 S = np.zeros([M, 31], float) 
     68 
     69mean_grad_ratio_zone = np.zeros([M, 31], float) # 2D-array of zonal mean gradient ratio for each day in each month  
     70std_grad_ratio_zone = np.zeros([M, 31], float) # 2D-array of zonal std of gradient ratio for each day in each month  
     71 
     72mean_spec_anom_23_zone = np.zeros([M, 31], float) # 2D-array of zonal mean emis_spec spatial anomaly for each day in each month (at 23GHz) 
     73std_spec_anom_23_zone = np.zeros([M, 31], float) # 2D-array of zonal std of emis_spec spatial anomaly for each day in each month (at 23GHz) 
     74 
     75mean_spec_anom_89_zone = np.zeros([M, 31], float ) # 2D-array of zonal mean emis_spec spatial anomaly for each day in each month (at 89GHz) 
     76std_spec_anom_89_zone = np.zeros([M, 31], float) # 2D-array of zonal std of emis_spec spatial anomaly for each day in each month (at 89GHz) 
     77 
     78mean_ratio_anom_89_zone = np.zeros([M, 31], float) # 2D-array of zonal mean emis_ratio spatial anomaly for each day in each month (at 89GHz) 
     79std_ratio_anom_89_zone = np.zeros([M, 31], float) # 2D-array of zonal std of emis_ratio spatial anomaly for each day in each month (at 89GHz) 
     80 
     81S = np.zeros([M, 31], float) # number of data pixels in selected zone, per day and per month 
    8182 
    8283for imo in range (0, M):  
    … …  
    165166 
    166167 
    167 ############################################################################ 
    168 # calculate standard deviation of emissivity parameters over the year 2009 # 
    169 ############################################################################ 
     168###################################################################################### 
     169# calculate standard deviation of emissivity parameters over the year 2009 (364 days)# 
     170###################################################################################### 
    170171print 'calculate standard deviation of emissivity parameters over the year 2009' 
    171172time_std1 = sqrt((1./364.)*sum((mean_year_grad_ratio[nonzero(isnan(mean_year_grad_ratio) == False)] - mean(mean_year_grad_ratio[nonzero(isnan(mean_year_grad_ratio) == False)]))**2)) 
    … …  
    244245title('area of study - Chukchi Sea') 
    245246plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/study_by_zones/map_grad_ratio_zone_Chukchi_Sea.png') 
    246 ''' 
    247  
    248  
    249 ################################ 
    250 # map parameters in all Arctic # 
    251 ################################ 
     247 
     248 
     249 
     250############################################ 
     251# read emissivity parameters in all Arctic # 
     252############################################ 
    252253cumul_params = np.zeros([M, ny, nx], float) 
    253254grad_ratio = np.zeros([M, ny, nx], float) 
    … …  
    275276    for ilon in range (0, nx): 
    276277        for ilat in range (0, ny): 
     278            # calculate monthly mean of emissivity parameters 
    277279            grad_ratio[imo, ilat, ilon] = mean(gr[:, ilat, ilon][nonzero(isnan(gr[:, ilat, ilon]) == False)]) 
    278280            spec_anom_23[imo, ilat, ilon] = mean(sa_23[:, ilat, ilon][nonzero(isnan(sa_23[:, ilat, ilon]) == False)]) 
    279281            spec_anom_89[imo, ilat, ilon] = mean(sa_89[:, ilat, ilon][nonzero(isnan(sa_89[:, ilat, ilon]) == False)]) 
    280282            ratio_anom_89[imo, ilat, ilon] = mean(ra_89[:, ilat, ilon][nonzero(isnan(ra_89[:, ilat, ilon]) == False)]) 
     283            # calculate monthly cumulation index = (sum(abs(all emissivity parameters))/maximum of this sum)*100 (to have a percentage) 
    281284            cumul_params[imo, ilat, ilon] = abs(grad_ratio[imo, ilat, ilon]) + abs(spec_anom_23[imo, ilat, ilon]) + abs(spec_anom_89[imo, ilat, ilon]) + abs(ratio_anom_89[imo, ilat, ilon]) 
    282285    CP[imo, :, :] = (cumul_params[imo, :, :]/max(cumul_params[imo, :, :][nonzero(isnan(cumul_params[imo, :, :]) == False)])) * 100. 
    283286 
    284287 
     288####################### 
     289# map cuulation index # 
     290####################### 
    285291#ion() 
    286292x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 

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

    r51 r56  
    4949    # 23GHz 
    5050    print 'read daily emis spec 23GHz' 
    51     fichier_23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') 
     51    fichier_23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')  
    5252    emis_spec_23 = fichier_23.variables['e_spec'][:] 
    5353    fichier_23.close() 

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

    r40 r56  
    1515from matplotlib import colors 
    1616 
    17  
     17############################################# 
     18# This function gives the length of the arc # 
     19# which links two points (1 and 2) defined  # 
     20# by their polar lon/lat coordinates        # 
     21############################################# 
    1822 
    1923def distance_on_unit_sphere(lat1, lon1, lat2, lon2): 

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

    r40 r56  
    1717 
    1818##################### 
    19 # read fichier .DAT # 
     19# read file .DAT # 
    2020##################### 
    2121fichier=open('/mma/hermozol/Documents/Data_tests/GLACE/GLACE_AMSUB89_EMIS_SEPTEMBER2009.DAT','r') 
    … …  
    9595####################################### 
    9696r = pi / 180. 
    97 opac = - cos(zen_angle * r) * np.log(transm)  
    98 E3 = scipy.special.expn(3, opac) 
    99 theta_eff = np.zeros([L], float) 
    100 tdn_lamb = np.zeros([L], float) 
    101 tb_spec = np.zeros([L], float) 
    102 tb_lamb = np.zeros([L], float) 
    103 e_spec = np.zeros([L], float) 
    104 e_lamb = np.zeros([L], float) 
    105 e_spec_lamb = np.zeros([L], float) 
     97opac = - cos(zen_angle * r) * np.log(transm) # opacity 
     98E3 = scipy.special.expn(3, opac) # exponentielle d ordre 3 
     99theta_eff = np.zeros([L], float) # effective incidence angle from Matzler method 
     100tdn_lamb = np.zeros([L], float) # downwelling radiation 
     101tb_spec = np.zeros([L], float) # specular brightness temperature 
     102tb_lamb = np.zeros([L], float) # lambertian brightness temperature 
     103e_spec = np.zeros([L], float) # specular emissivity 
     104e_lamb = np.zeros([L], float) # lambertian emissivity 
     105e_spec_lamb = np.zeros([L], float) # difference between lambertian and specular emissivity 
    106106for ii in range (0, L): 
    107107    theta_eff[ii] = math.acos(- opac[ii] / (np.log(2. * E3[ii]))) * (1. / r) 

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

    r40 r56  
    1 # griddata.py - 2010-07-11 ccampo 
    21import numpy as np 
    32from numpy import * 

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

    r40 r56  
    1515 
    1616 
     17 
     18 
     19####################################################################### 
     20# Open .dat file to read emis spec and emis lamb in September and fov # 
     21####################################################################### 
    1722fichier=open('/mma/hermozol/Documents/Data_tests/monthly_GLACE/lamb_spec_param_all_angles_SEPTEMBER2009.dat','r') 
    1823numlines = 0 
    … …  
    2126 
    2227fichier.close 
    23  
    2428fichier = open('/mma/hermozol/Documents/Data_tests/monthly_GLACE/lamb_spec_param_all_angles_SEPTEMBER2009.dat','r')    
    2529nbtotal = numlines-1 
    … …  
    4044     iligne=iligne+1 
    4145 
    42  
    43  
    4446fichier.close() 
    4547 
    46  
    47  
     48# compute fov as x axis in scatter plot 
    4849x0 = fov.min() 
    4950x1 = fov.max() 
     51# compute diff emis_lamb-emis_spec as y axis in scatter plot 
    5052y0 = e_spec_lamb.min() 
    5153y1 = e_spec_lamb.max() 
     54# to each observation of coord x and y, add z = value of emis_spec 
    5255z0 = e_spec.min() 
    5356z1 = e_spec.max() 
     57# step of x and y axis 
    5458dx = 1. 
    5559dy = 0.001 
     60# use ffgrid to create a new grid of fov/diff(lamb-spec) 
    5661zgrid, xvec, yvec = ffgrid2.ffgrid(fov, e_spec_lamb, e_spec, dx, dy, x0, x1, y0, y1, z0, z1) 
    5762xii, yii = np.meshgrid(xvec, yvec) 
     63 
     64############# 
     65# draw grid # 
     66############# 
    5867plt.ion() 
    5968figure() 

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

    r54 r56  
    4141 
    4242 
    43 frequ = 89 
     43frequ = 89 # define which frequency to use for the study 
    4444 
    4545##################### 
    4646# read NETCDF files # 
    4747##################### 
    48 emis_spec_month = np.zeros([M, ny, nx], float) 
    49 emis_lamb_month = np.zeros([M, ny, nx], float) 
    50 ratio_emis = np.zeros([M, ny, nx], float) 
     48 
     49emis_spec_month = np.zeros([M, ny, nx], float) # monthly mean of emis_spec 
     50emis_lamb_month = np.zeros([M, ny, nx], float) # monthly mean of emis_lamb 
     51ratio_emis = np.zeros([M, ny, nx], float) # monthly mean of emis_ratio 
     52 
     53# conversion of previous data into 1D arrays 
    5154ratio_emis_vec = np.zeros([M, ny * nx], float) 
    5255spec_emis_vec = np.zeros([M, ny * nx], float) 
    5356lamb_emis_vec = np.zeros([M, ny * nx], float) 
     57 
    5458print 'read data from files' 
    5559for imo in range (0, M): 
    … …  
    9397''' 
    9498 
    95  
    96 c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) 
     99############################################################## 
     100# definition of parameters of the histograms of distribution # 
     101############################################################## 
     102c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) # definition of one color per month 
    97103#limit_coef_spec = np.array([0.6, 0.6, 0.7, 0.8, 0.75]) # i1 = AMSUA23GHz / i2 = AMSUA30GHz / i3 = AMSUA50GHz / i4 = AMSUA89GHz / i5 = AMSUB89GHz 
    98104#limit_coef_lamb = np.array([0.6, 0.6, 0.8, 0.8, 0.77]) # i1 = AMSUA23GHz / i2 = AMSUA30GHz / i3 = AMSUA50GHz / i4 = AMSUA89GHz / i5 = AMSUB89GHz 
    99 #idata = 0 
    100 fontP = FontProperties() 
     105#idata = 0 # defines which frequency we use (0=A23, 1=A30, 2=A50, 3=A89, 4=B89) 
     106fontP = FontProperties() # concerns legend font in plots 
    101107fontP.set_size('small') 
     108 
     109 
    102110print 'distribution of emissivity values' 
    103111################################### 
    … …  
    108116corresp_emis_spec = np.zeros([M, 50], float) 
    109117for imo in range (0, M): 
    110     hist_vals_spec[imo, :] = hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 
     118    hist_vals_spec[imo, :] = hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] # calculate histogram distribution of emissivity 
    111119    for ibin in range (0, 50): 
    112         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]) 
     120        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]) # calculate corresponding emissivity value to histogram distribution (in frequency of occurence) 
    113121 
    114122## plot first six months of spec emissivity histograms ## 

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

    r40 r56  
    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 
    1210import scipy.special 
    1311import ffgrid2 
    … …  
    2220 
    2321 
    24 dx=1. 
    25 dy=0.5 
     22###################################### 
     23# polar lon/lat grid characteristics # 
     24###################################### 
     25dx=1. # longitude every 1 deg 
     26dy=0.5 # latitude every 0.5 deg 
    2627x0, x1 = -180., 180. 
    27 y0, y1 = 65., 90. 
     28y0, y1 = 65., 90. # only high latitudes 
    2829xvec = np.arange(x0, x1 + dx, dx) 
    2930yvec = np.arange(y0, y1 + dy, dy) 
    … …  
    3233 
    3334 
    34 e_spec_lamb_month = np.zeros([ny, nx, M], float) 
     35#################################################################################### 
     36# read emis_spec, emis_lamb and emis_diff data in monthly NETCDF files for mapping # 
     37#################################################################################### 
     38 
     39e_spec_lamb_month = np.zeros([ny, nx, M], float) # monthly mean of diff(emis_lamb-emis_spec) 
    3540e_spec_month = np.zeros([ny, nx, M], float) 
     41 
    3642for imo in range (0, M): 
    3743    fichier = Dataset('/mma/hermozol/Documents/Data_tests/monthly_GLACE/gridded_data/grid_monthly_data_lamb_spec_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 
    … …  
    4551    e_lamb = fichier.variables['e_lamb'][:] 
    4652    e_spec_lamb = fichier.variables['e_spec_lamb'][:] 
    47     e_sl_00 = fichier.variables['e_mixed_s00'][:] 
    48     e_sl_25 = fichier.variables['e_mixed_s25'][:] 
    49     e_sl_50 = fichier.variables['e_mixed_s50'][:] 
    50     e_sl_75 = fichier.variables['e_mixed_s75'][:] 
    51     e_sl_100 = fichier.variables['e_mixed_s100'][:] 
     53    e_sl_00 = fichier.variables['e_mixed_s00'][:] # emissivity specular and lambertian with specular coefficient of 0% (only lambertian emissivity) 
     54    e_sl_25 = fichier.variables['e_mixed_s25'][:] # emissivity specular and lambertian with specular coefficient of 25%  
     55    e_sl_50 = fichier.variables['e_mixed_s50'][:] # emissivity specular and lambertian with specular coefficient of 50%  
     56    e_sl_75 = fichier.variables['e_mixed_s75'][:] # emissivity specular and lambertian with specular coefficient of 75%  
     57    e_sl_100 = fichier.variables['e_mixed_s100'][:] # emissivity specular and lambertian with specular coefficient of 100% (only specular emissivity)  
    5258    fichier.close() 
    53     e_spec_month[:, :, imo] = e_spec[:, :] 
    54     e_spec_lamb_month[:, :, imo] = e_spec_lamb[:, :] 
     59    e_spec_month[:, :, imo] = e_spec[:, :] # monthly mean of emis_spec 
     60    e_spec_lamb_month[:, :, imo] = e_spec_lamb[:, :] # monthly mean of diff(emis_lamb-emis_spec) 
    5561 
    5662 
    57  
     63####### 
     64# map # 
     65####### 
     66print 'start mapping' 
    5867cm = plt.cm.get_cmap('RdBu') 
    5968#cmap = cm.s3pcpn_l_r 
    6069plt.ion() 
    6170for imo in range (0, M): 
     71    print 'month ' + month[imo] 
    6272    map_ffgrid.draw_map_l(lon, lat, 0.6, 1.02, 0.01, e_spec_month[:,:,imo], cm.s3pcpn_l_r, 'emissivity SPEC') 
    6373    title(month[imo] + ' 2009 - AMSUB') 

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

    r55 r56  
    4141 
    4242 
    43 frequ = np.array(['23', '30', '50', '89']) 
     43frequ = np.array(['23', '30', '50', '89']) # choose studied frequency  
    4444 
     45 
     46####################################################################### 
     47# Reads monthly mean and std of emis spec, lamb and diff NETCDF files # 
     48####################################################################### 
    4549mean_spec = np.zeros([4, M, ny, nx], float) 
    4650mean_lamb = np.zeros([4, M, ny, nx], float) 
    4751mean_lamb_spec = np.zeros([4, M, ny, nx], float) 
     52 
    4853std_spec = np.zeros([4, M, ny, nx], float) 
    4954std_lamb = np.zeros([4, M, ny, nx], float) 
     55 
    5056for ifr in range (0, 4): 
    5157    print 'frequency ' + frequ[ifr] 
    … …  
    6167 
    6268 
    63 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, mean_spec[0, 10, :, :], 0.45, 1.02, 0.001, colormap, 'emissivity LAMB') 
    6469 
    6570 
    66  
     71################################################# 
     72# map of spec, lamb or diff of emis std or mean # 
     73################################################# 
    6774ion() 
    6875x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 

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

    r40 r56  
    1010 
    1111 
    12  
     12###################################################### 
     13# Function that maps 2D-array data in a lon/lat grid # 
     14# computed in ffgrid function                        # 
     15###################################################### 
    1316 
    1417 
    1518def draw_map_l (xcoord, ycoord, c0, c1, dc, zvar, colormap, text): 
    16  
     19    
     20    ############################################################## 
     21    # xcoord = vector of longitudes                              # 
     22    # ycoord = vector of latitudes                               # 
     23    # c0, c1, dc = boundaries and step of scale of the colorbar  #  
     24    # zvar = 2D-array of data to map                             # 
     25    # colormap = colors to use for mapping od data               # 
     26    # text = title of the colorbar                               # 
     27    ############################################################## 
     28    
    1729    figure() 
    1830    m = Basemap(projection = 'nplaea', boundinglat = 65., lon_0 = 0., resolution = 'l', fix_aspect = True) 

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

    r51 r56  
    6363        print 'stack in file month ' + str(month[imo]) 
    6464        print 'open anomaly file' 
    65         fichier_anom = 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', 'r', format='NETCDF3_CLASSIC') 
     65        fichier_anom = 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', 'r', format='NETCDF3_CLASSIC') 
    6666        anom_spec[ifr, imo, :, :] = fichier_anom.variables['spec_anomaly'][:] 
    6767        anom_lamb[ifr, imo, :, :] = fichier_anom.variables['lamb_anomay'][:] 
    … …  
    7070        fichier_anom.close() 
    7171        print 'open ice file' 
    72         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') 
     72        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_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 
    7373        ice_spec[ifr, imo, :, :] = fichier_ice.variables['spec_ice_area'][:] 
    7474        ice_lamb[ifr, imo, :, :] = fichier_ice.variables['lamb_ice_area'][:] 
    … …  
    9999 
    100100 
    101 ################################################ 
    102 # calculate gradient ratio for two frequencies # 
    103 ################################################ 
     101#################################################### 
     102# calculate gradient ratio between two frequencies # 
     103#################################################### 
    104104print 'calculate gradient ratio' 
    105105grad_ratio = np.zeros([M, ny, nx], float) 
    … …  
    111111        for ilon in range (0, nx): 
    112112            grad_ratio[imo, ilat, ilon] = (ice_spec[0, imo, ilat, ilon] - ice_spec[3, imo, ilat, ilon]) / (ice_spec[0, imo, ilat, ilon] + ice_spec[3, imo, ilat, ilon]) 
    113     ''' 
    114     # compute histogram distribution of gradient ratio to find characteristic threshold 
    115     grad_ratio_vec = reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))[nonzero(isnan(reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))) == False)] 
    116     hist_ratio[imo, :] = hist(grad_ratio_vec, bins = 50, normed = True, histtype='step')[0] 
    117     for ibin in range (0, 50): 
    118         corresp_ratio[imo, ibin] = mean(hist(grad_ratio_vec, bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 
    119  
    120 # plot histogram 
    121 c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) 
    122 figure() 
    123 for imo in range (0, 6): 
    124     plot(corresp_ratio[imo], hist_ratio[imo], c = str(c[imo]), label = str(month[imo])) 
    125  
    126 grid() 
    127 xlim(corresp_ratio.min(), corresp_ratio.max()) 
    128 xlabel('emissivity ratio') 
    129 ylabel('frequency of occurence') 
    130 fontP = FontProperties() 
    131 fontP.set_size('small') 
    132 legend(loc = 2, prop = fontP) 
    133 #plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_spec_AMSUA' + str(frequ) + '_JANUARY-JUNE_2009.png') 
    134 ## plot six following months of spec emissivity histograms ## 
    135 figure() 
    136 for imo in range (6, M): 
    137     plot(corresp_ratio[imo], hist_ratio[imo], c = str(c[imo - 6]), label = str(month[imo])) 
    138  
    139 grid() 
    140 xlim(corresp_ratio.min(), corresp_ratio.max()) 
    141 xlabel('emissivity ratio') 
    142 ylabel('frequency of occurence') 
    143 legend(loc = 1, prop = fontP) 
    144 #plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_spec_AMSUA' + str(frequ) + '_JULY-DECEMBER_2009.png') 
    145 ''' 
     113 
    146114 
    147115################################ 
    … …  
    184152 
    185153 
    186  
    187 ''' 
    188 figure() 
    189 scatter(x_coast, y_coast, c = z_coast, s = volume) 
    190 #cmap = cm.s3pcpn_l_r 
    191 levels = np.arange(-1, 11, 1) 
    192 cs = contour(xvec, yvec, cumul[imo, :, :])  
    193 cs = contour(xvec, yvec, new_anom_spec0[imo, :, :], levels)  
    194 cs = contour(xvec, yvec, new_anom_ratio[imo, :, :], levels)  
    195 cs = contour(xvec, yvec, new_grad_ratio[imo, :, :], levels)  
    196 cbar = colorbar(cs, ticks = [-1, 0, 1, 2, 3, 4, 5]) 
    197 cbar.ax.set_yticklabels(['emis_spec_89', 'emis_spec_23', 'emis_ratio_89', 'grad_ratio_23-89', '-']) 
    198 #cbar.set_label(text) 
    199 xlim(-3500, 2700.) 
    200 ylim(-4000, 2800.) 
    201 ''' 
    202  
    203  
    204  
    205  
    206  
    207  
    208  
    209 ''' 
    210 ############################################ 
    211 # time evolution (monthly) in a given zone # 
    212 ############################################ 
    213 #zone 1 : yi = 400. // yf = 720. // xi = -1080. // xf = -920. (Canadian Archipelago) 
    214 #zone 2 : yi = 920. // yf = 1280. // xi = -720. // xf = -560. (Chukchi Sea, Bering Detroit) 
    215 # select borders of zone 
    216 yi = 920.  
    217 yf = 1280. 
    218 xi = -720. 
    219 xf = -560. 
    220 #find corresponding index in xvec and yvec 
    221 xxi = np.where(xvec == xi)[0][0] 
    222 xxf = np.where(xvec == xf)[0][0] 
    223 yyi = np.where(yvec == yi)[0][0] 
    224 yyf = np.where(yvec == yf)[0][0] 
    225  
    226 # define vectors 
    227 anom_spec0_zone = np.zeros([M], float) 
    228 anom_spec3_zone = np.zeros([M], float) 
    229 anom_ratio_zone = np.zeros([M], float) 
    230 grad_ratio_zone = np.zeros([M], float) 
    231 std_as0 = np.zeros([M], float) 
    232 std_as3 = np.zeros([M], float) 
    233 std_ar = np.zeros([M], float) 
    234 std_gr = np.zeros([M], float) 
    235 for imo in range (0, M): 
    236     # anom spec emis from 23GHz 
    237     anom_spec0_vec = reshape(anom_spec[0, imo, yyi : yyf + 1, xxi : xxf + 1], size(anom_spec[0, imo, yyi : yyf + 1, xxi : xxf + 1])) 
    238     if (len(anom_spec0_vec[nonzero(isnan(anom_spec0_vec) == True)]) != 0): 
    239         anom_spec0_zone[imo] = NaN 
    240         std_as0[imo] = NaN 
    241         anom_spec3_zone[imo] = NaN 
    242         std_as3[imo] = NaN 
    243         anom_ratio_zone[imo] = NaN 
    244         std_ar[imo] = NaN 
    245         grad_ratio_zone[imo] = NaN 
    246         std_gr[imo] = NaN 
    247         continue 
    248     else: 
    249         anom_spec0_zone[imo] = anom_spec0_vec.mean() 
    250         std_as0[imo] = sqrt((1. / (size(anom_spec[0, imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((anom_spec0_vec[:] - anom_spec0_zone[imo])**2)) 
    251         # anom spec emis from 89GHz 
    252         anom_spec3_vec = reshape(anom_spec[3, imo, yyi : yyf + 1, xxi : xxf + 1], size(anom_spec[3, imo, yyi : yyf + 1, xxi : xxf + 1])) 
    253         anom_spec3_zone[imo] = anom_spec3_vec.mean() 
    254         std_as3[imo] = sqrt((1. / (size(anom_spec[3, imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((anom_spec3_vec[:] - anom_spec3_zone[imo])**2)) 
    255         # anom emis ratio from 89GHz 
    256         anom_ratio_vec = reshape(anom_ratio[3, imo, yyi : yyf + 1, xxi : xxf + 1], size(anom_ratio[3, imo, yyi : yyf + 1, xxi : xxf + 1])) 
    257         anom_ratio_zone[imo] = anom_ratio_vec.mean() 
    258         std_ar[imo] = sqrt((1. / (size(anom_ratio[3, imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((anom_ratio_vec[:] - anom_ratio_zone[imo])**2)) 
    259         # gradient ratio from 23 and 89GHz 
    260         grad_ratio_vec = reshape(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1], size(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1])) 
    261         grad_ratio_zone[imo] = grad_ratio_vec.mean() 
    262         std_gr[imo] = sqrt((1. / (size(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((grad_ratio_vec[:] - grad_ratio_zone[imo])**2)) 
    263  
    264  
    265 figure() 
    266 fontP = FontProperties() 
    267 fontP.set_size('small') 
    268 subplot(2,1,1) 
    269 errorbar(arange(0, M, 1), anom_spec0_zone, yerr = std_as0, label = 'emis spec anomaly 23GHz') 
    270 errorbar(arange(0, M, 1), anom_spec3_zone, yerr = std_as3, label = 'emis spec anomaly 89GHz') 
    271 errorbar(arange(0, M, 1), grad_ratio_zone, yerr = std_gr, label = 'gradient ratio 23-89GHz') 
    272 plot(arange(0, M, 1), np.zeros([M], float), '--k') 
    273 legend(loc = 3, prop = fontP) 
    274 grid() 
    275 xticks(arange(0, M, 1), month, rotation = 15) 
    276 xlim(-1, M) 
    277 subplot(2,1,2) 
    278 errorbar(arange(0, M, 1), anom_ratio_zone, yerr = std_ar, label = 'emis ratio anomaly 89GHz') 
    279 plot(arange(0, M, 1), np.zeros([M], float), '--k') 
    280 legend(loc = 2, prop = fontP) 
    281 grid() 
    282 xticks(arange(0, M, 1), month, rotation = 15) 
    283 xlim(-1, M) 
    284 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/study_by_zones/monthly_evolution_emis_params_zone_North_Beaufort_Sea.png') 
    285 ### map the selected zone ### 
    286 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec[xxi : xxf], yvec[yyi : yyf], grad_ratio[imo, yyi : yyf, xxi : xxf], grad_ratio[imo, yyi : yyf, xxi : xxf][nonzero(isnan(grad_ratio[imo, yyi : yyf, xxi : xxf]) == False)].min(), grad_ratio[imo, yyi : yyf, xxi : xxf][nonzero(isnan(grad_ratio[imo, yyi : yyf, xxi : xxf]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'selected zone') 
    287 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/study_by_zones/map_zone_North_Beaufort_Sea.png') 
    288 ''' 
    289  
    290  
    291  
    292 ''' 
    293154############################################################################ 
    294155# calculate correlation between anom ratio and gradient ratio (spec emiss) # 
    … …  
    304165        for jj in range (0, 4): 
    305166            corr_mat[imo, ii, jj] = corrcoef(params[ii], params[jj])[0][1] 
     167 
    306168 
    307169figure() 
    … …  
    322184grid() 
    323185plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/correlation_params_2009.png') 
    324 ''' 
     186 
    325187 
    326188''' 

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

    r55 r56  
    5757ice_diff_anom = np.zeros([4, M, ny, nx], float) 
    5858ice_ratio_anom = np.zeros([4, M, ny, nx], float) 
     59 
    5960for ifr in range (3, 4): 
    6061    print 'frequency ' + frequ[ifr] 
    … …  
    6970        for ilat in range (0, ny): 
    7071            for ilon in range (0, nx): 
     72                # calculate monthly mean of emis spec, lamb, diff and ratio 
    7173                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)]) 
    7274                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)]) 
    7375                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)]) 
    7476                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)]) 
     77        # transforms 2D-arrays of monthly mean of emissivity into vectors, in order to calculate the average emissivity over total sea ice extent, for each month  
    7578        a = reshape(mean_ice_spec[ifr, imo, :, :], size(mean_ice_spec[ifr, imo, :, :])) 
    7679        b = reshape(mean_ice_lamb[ifr, imo, :, :], size(mean_ice_lamb[ifr, imo, :, :])) 
    … …  
    8083        for ilat in range (0, ny): 
    8184            for ilon in range (0, nx): 
     85                # calculate monthly spatial anomaly in each pixel 
    8286                ice_spec_anom[ifr, imo, ilat, ilon] = mean_ice_spec[ifr, imo, ilat, ilon] - mean(a[nonzero(isnan(a) == False)]) 
    8387                ice_lamb_anom[ifr, imo, ilat, ilon] = mean_ice_lamb[ifr, imo, ilat, ilon] - mean(b[nonzero(isnan(b) == False)]) 
    … …  
    107111 
    108112 
    109 ''' 
    110  
    111 # test: 
    112 ion() 
    113 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    114 x_coast = x_ind 
    115 y_coast = y_ind 
    116 z_coast = z_ind 
    117 for ifr in range (0, 2): 
    118     for imo in range (0, M): 
    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') 
    120         title('AMSUA ' + str(frequ) + ' - ' + str(month[imo]) + ' 2009') 
    121         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') 
    122 #cm.s3pcpn_l_r 
    123  
    124 ''' 
    125113 
    126114 
    127  
    128  
    129 ''' 
    130 frequ = 89 
    131 spec_ice = np.zeros([M, ny, nx], float) 
    132 lamb_ice = np.zeros([M, ny, nx], float) 
    133 diff_ice = np.zeros([M, ny, nx], float) 
    134 spec_month = np.zeros([M, ny, nx], float) 
    135 lamb_month = np.zeros([M, ny, nx], float) 
    136 diff_month = np.zeros([M, ny, nx], float) 
    137 std_spec = np.zeros([M, ny, nx], float) 
    138 std_lamb = np.zeros([M, ny, nx], float) 
    139 std_diff = np.zeros([M, ny, nx], float) 
    140 for imo in range (0, 1): 
    141     print 'month ' + month[imo] 
    142     ################################################################################## 
    143     # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study # 
    144     ################################################################################## 
    145     print 'open threshold file' 
    146     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') 
    147     spec_lim = fichier_emis.variables['spec_ice_area'][:] 
    148     fichier_emis.close() 
    149     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') 
    150     day = fichier_e.variables['days'][:] 
    151     emis_spec = fichier_e.variables['e_spec'][:] 
    152     emis_lamb = fichier_e.variables['e_lamb'][:] 
    153     emis_diff = fichier_e.variables['e_spec_lamb'][:] 
    154     fichier_e.close() 
    155     for ilon in range (0, nx): 
    156         for ilat in range (0, ny): 
    157            if (isnan(spec_lim[ilat, ilon]) == True): 
    158                 spec_ice[imo, ilat, ilon] = NaN 
    159                 lamb_ice[imo, ilat, ilon] = NaN 
    160                 diff_ice[imo, ilat, ilon] = NaN 
    161             else: 
    162                 spec_ice[imo, ilat, ilon] = spec_month[imo, ilat, ilon] 
    163                 lamb_ice[imo, ilat, ilon] = lamb_month[imo, ilat, ilon] 
    164                 diff_ice[imo, ilat, ilon] = diff_month[imo, ilat, ilon] 
    165             spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, :][nonzero(isnan(emis_spec[ilat, ilon, :]) == False)]) 
    166             lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, :][nonzero(isnan(emis_lamb[ilat, ilon, :]) == False)]) 
    167             diff_month[imo, ilat, ilon] = mean(emis_diff[ilat, ilon, :][nonzero(isnan(emis_diff[ilat, ilon, :]) == False)]) 
    168             Ss = sum((emis_spec[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_spec[ilat, ilon, 0 : month_day[imo]]) == False)] - spec_month[imo, ilat, ilon])**2) 
    169             std_spec[imo, ilat, ilon] = sqrt((1. / month_day[imo] - 1.) * Ss) 
    170             Sl = sum((emis_lamb[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_lamb[ilat, ilon, 0 : month_day[imo]]) == False)] - lamb_month[imo, ilat, ilon])**2) 
    171             std_lamb[imo, ilat, ilon] = sqrt((1. / month_day[imo] - 1.) * Sl) 
    172      
    173  
    174  
    175  
    176 # test: 
     115####################### 
     116# map of emis anomaly # 
     117####################### 
    177118ion() 
    178119x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 
    … …  
    186127#cm.s3pcpn_l_r 
    187128 
    188 imo = 1 
    189 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec[imo, :, :], std_spec[imo,:,:][nonzero(isnan(std_spec[imo,:,:]) == False)].min(), std_spec[imo, :, :][nonzero(isnan(std_spec[imo , :, :]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'std spec') 
    190 ''' 
     129 

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

    r40 r56  
    3838 
    3939 
    40 osi_type = np.zeros([M, 31, ny, nx], float) 
    41 spec = np.zeros([M, 31, ny, nx], float) 
    42 lamb = np.zeros([M, 31, ny, nx], float) 
     40osi_type = np.zeros([M, 31, ny, nx], float) # daily ice type in polar lon/lat grid 
     41spec = np.zeros([M, 31, ny, nx], float) # daily emissivity spec in polar lon/lat grid 
     42lamb = np.zeros([M, 31, ny, nx], float) # daily emissivity lamb in polar lon/lat grid 
     43 
    4344for imo in range (3, 4): 
    4445    print 'month ' + month[imo] 
    … …  
    7677         iligne=iligne+1 
    7778    fichier_amsub.close() 
    78     for ijr in range (0, month_day[imo]): 
     79    for ijr in range (0, month_day[imo]): # read daily data 
     80        # stop computing for this month if day > maximum day in month 
    7981        if ((ijr + 1) > month_day[imo]): continue 
    8082        else: 
    … …  
    8486               continue 
    8587           else: 
     88               # re-grid in polar lon/lat daily emis spec and emis_lamb  
    8689               print 'date ', jjr[bbjr][0]  
    8790               lat_jr = lat[bbjr] 
    … …  
    129132 
    130133 
    131  
     134''' 
    132135######### 
    133136# tests # 
    … …  
    140143cmap2 = colors.ListedColormap(['1.', '0.25', '0.50', '0.75']) 
    141144map_ffgrid.draw_map_l (xvec, yvec, 0, 5, 1, osi_type[11, 19, :, :] , cmap2, 'OSISAF ice type') 
     145''' 
    142146 
    143147 
    … …  
    148152 
    149153 
    150  

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

    r40 r56  
    1717month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) 
    1818M = len(month) 
     19 
     20 
    1921 
    2022 
    … …  
    9597 
    9698 
    97 bbjr = nonzero(jjr == 1.) 
    98 bbzen = nonzero(zen[bbjr] > 20.) 
    99 z0 = e[nonzero(isnan(e) == False)].min() 
    100 z1 = e[nonzero(isnan(e) == False)].max() 
    101 zgrid, xvec, yvec = ffgrid2.ffgrid(llon, llat, e, 1., 1., -180., 180., 65., 90., z0, z1) 
    102 emis_amsub_89 = zgrid   
    103 plt.ion() 
    104  
    105 map_ffgrid.draw_map_l(, lat, emis.min(), emis.max(), 0.01, emis, 'emissivity') 
    106          
    107  
    108  
    109  
    110  
    111  
    11299                 
    113100 

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

    r53 r56  
    3131x0 = -3000. # min limit of grid 
    3232x1 = 2500. # max limit of grid 
    33 dx = 100. 
     33dx = 40. # x resolution of grid 
    3434xvec = np.arange(x0, x1+dx, dx) 
    3535nx = len(xvec)  
    3636y0 = -3000. # min limit of grid 
    3737y1 = 3000. # max limit of grid 
    38 dy = 100. 
     38dy = 40. # y resolution of grid 
    3939yvec = np.arange(y0, y1+dy, dy) 
    4040ny = len(yvec) 
Note: See TracChangeset for help on using the changeset viewer.