Changeset 50

Timestamp:

08/06/14 18:06:12 (10 years ago)

Author:

lahlod

Message:

new scripts

Location:

trunk/src/scripts_Laura/ARCTIC/Travail_CEN

Files:

• map_emis_ratio_arctic.py (added)
• map_statistic_parameters.py (modified) (3 diffs)
• read_cartesian_gridded_AMSUA23.py (modified) (2 diffs)

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

@@ -78 +78 @@
 
 
-################################################
-# calculate gradient ratio for two frequencies #
-################################################
-grad_ratio = np.zeros([M, ny, nx], float)
-for imo in range (0, M):
-    for ilat in range (0, ny):
-        for ilon in range (0, nx):
-            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])
-     
-
-# test:
+'''
+###############
+# map anomaly #
+###############
+print 'map monthly anomaly'
 ion()
 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
@@ -94 +88 @@
 y_coast = y_ind
 z_coast = z_ind
-ifr = 0
-for imo in range (0, M):
-    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, anom_ratio[ifr, imo, :, :], anom_ratio[ifr, imo,:,:][nonzero(isnan(anom_ratio[ifr, imo,:,:]) == False)].min(), anom_ratio[ifr, imo, :, :][nonzero(isnan(anom_ratio[ifr, imo , :, :]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'Anomaly of emissivity ratio')
+ifr = 0 # anomaly of emissivity ratio from AMSUA89GHz database
+for imo in range (0, M):
+    print 'month ' + str(month[imo])
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, anom_spec[ifr, imo, :, :], anom_spec[ifr, imo,:,:][nonzero(isnan(anom_spec[ifr, imo,:,:]) == False)].min(), anom_spec[ifr, imo, :, :][nonzero(isnan(anom_spec[ifr, imo , :, :]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'Anomaly of emissivity spec')
     title('AMSUA ' + str(frequ[ifr]) + 'GHz - ' + str(month[imo]) + ' 2009')
-    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/anomaly_ratio_map_AMSUA' + str(frequ[ifr]) + '_with_AMSUA23-and-30-spec_'+str(month[imo])+'_2009.png')
+    print 'save figure in .png file'
+    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_anomaly/spec/anomaly_spec_map_AMSUA' + str(frequ[ifr]) + '_with_AMSUA23-and-30-spec_'+str(month[imo])+'_2009.png')
 #cm.s3pcpn_l_r
-
-
-
-
+'''
+
+
+################################################
+# calculate gradient ratio for two frequencies #
+################################################
+print 'calculate gradient ratio'
+grad_ratio = np.zeros([M, ny, nx], float)
+hist_ratio = np.zeros([M, 50], float)
+corresp_ratio = np.zeros([M, 50], float)
+for imo in range (0, M):
+    print 'month ' + month[imo]
+    for ilat in range (0, ny):
+        for ilon in range (0, nx):
+            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])
+    '''
+    # compute histogram distribution of gradient ratio to find characteristic threshold
+    grad_ratio_vec = reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))[nonzero(isnan(reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))) == False)]
+    hist_ratio[imo, :] = hist(grad_ratio_vec, bins = 50, normed = True, histtype='step')[0]
+    for ibin in range (0, 50):
+        corresp_ratio[imo, ibin] = mean(hist(grad_ratio_vec, bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
+
+# plot histogram
+c = np.array(['r', 'b', 'c', 'm', 'y', 'g'])
+figure()
+for imo in range (0, 6):
+    plot(corresp_ratio[imo], hist_ratio[imo], c = str(c[imo]), label = str(month[imo]))
+
+grid()
+xlim(corresp_ratio.min(), corresp_ratio.max())
+xlabel('emissivity ratio')
+ylabel('frequency of occurence')
+fontP = FontProperties()
+fontP.set_size('small')
+legend(loc = 2, prop = fontP)
+#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')
+## plot six following months of spec emissivity histograms ##
+figure()
+for imo in range (6, M):
+    plot(corresp_ratio[imo], hist_ratio[imo], c = str(c[imo - 6]), label = str(month[imo]))
+
+grid()
+xlim(corresp_ratio.min(), corresp_ratio.max())
+xlabel('emissivity ratio')
+ylabel('frequency of occurence')
+legend(loc = 1, prop = fontP)
+#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')
+'''
+
+################################
+# map cumulation of parameters #
+################################
+cumul1 = anom_spec[0, :, :, :] - grad_ratio # correlation
+cumul2 = anom_spec[3, :, :, :] + grad_ratio # anti correlation
+cumul3 = anom_spec[3, :, :, :] + anom_ratio[3, :, :, :] # anti correlation
+cumul4 = abs(anom_spec[0, :, :, :]) + abs(anom_spec[3, :, :, :]) + abs(grad_ratio) + abs(anom_ratio[3, :, :, :])
+
+
+
+
+ion()
+x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
+x_coast = x_ind
+y_coast = y_ind
+z_coast = z_ind
+for imo in range (0, M):
+    '''
+    # map cumul1
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, cumul1[imo, :, :], -0.25, 0.14, 0.01, cm.s3pcpn_l_r, 'Cumulation [emis spec anomaly 23 - grad ratio 23-89]')
+    title('AMSUA - ' + month[imo] + ' 2009')
+    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/cumul_params/map_cumulation_spec-anom-23_grad-ratio_' + month[imo] +'2009.png')
+    # map cumul2
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, cumul2[imo, :, :], -0.17, 0.17, 0.01, cm.s3pcpn_l_r, 'Cumulation [emis spec anomaly 89 + grad ratio 23-89]')
+    title('AMSUA - ' + month[imo] + ' 2009')
+    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/cumul_params/map_cumulation_spec-anom-89_grad-ratio_' + month[imo] +'2009.png')
+    # map cumul3
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, cumul3[imo, :, :], -5.46, 2.14, 0.01, cm.s3pcpn_l_r, 'Cumulation [emis spec anomaly 89 + emis ratio 89]')
+    title('AMSUA - ' + month[imo] + ' 2009')
+    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/cumul_params/map_cumulation_spec-anom-89_ratio-anom-89_' + month[imo] +'2009.png')
+    '''
+    # map cumul4 
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, cumul4[imo, :, :], cumul4[imo, :, :][nonzero(isnan(cumul4[imo, :, :]) == False)].min(), cumul4[imo, :, :][nonzero(isnan(cumul4[imo, :, :]) == False)].max(), 0.1, cm.s3pcpn_l_r, 'Cumulation all parameters')
+    title('AMSUA - ' + month[imo] + ' 2009')
+    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/cumul_params/map_cumulation_all_parameters_' + month[imo] +'2009.png')
+
+
+
+
+
+'''
+figure()
+scatter(x_coast, y_coast, c = z_coast, s = volume)
+#cmap = cm.s3pcpn_l_r
+levels = np.arange(-1, 11, 1)
+cs = contour(xvec, yvec, cumul[imo, :, :]) 
+cs = contour(xvec, yvec, new_anom_spec0[imo, :, :], levels) 
+cs = contour(xvec, yvec, new_anom_ratio[imo, :, :], levels) 
+cs = contour(xvec, yvec, new_grad_ratio[imo, :, :], levels) 
+cbar = colorbar(cs, ticks = [-1, 0, 1, 2, 3, 4, 5])
+cbar.ax.set_yticklabels(['emis_spec_89', 'emis_spec_23', 'emis_ratio_89', 'grad_ratio_23-89', '-'])
+#cbar.set_label(text)
+xlim(-3500, 2700.)
+ylim(-4000, 2800.)
+'''
+
+
+
+
+
+
+
+'''
+############################################
+# time evolution (monthly) in a given zone #
+############################################
+# read monthly std of emis 
+ifr = 3
+fichier_temporal_std = 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_AMSUB' + str(frequ[ifr]) + '_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC')
+
+
+
+#zone 1 : yi = 400. // yf = 720. // xi = -1080. // xf = -920. (Canadian Archipelago)
+#zone 2 : yi = 920. // yf = 1280. // xi = -720. // xf = -560. (Chukchi Sea, Bering Detroit)
+# select borders of zone
+yi = 920. 
+yf = 1280.
+xi = -720.
+xf = -560.
+#find corresponding index in xvec and yvec
+xxi = np.where(xvec == xi)[0][0]
+xxf = np.where(xvec == xf)[0][0]
+yyi = np.where(yvec == yi)[0][0]
+yyf = np.where(yvec == yf)[0][0]
+
+# define vectors
+anom_spec0_zone = np.zeros([M], float)
+anom_spec3_zone = np.zeros([M], float)
+anom_ratio_zone = np.zeros([M], float)
+grad_ratio_zone = np.zeros([M], float)
+std_as0 = np.zeros([M], float)
+std_as3 = np.zeros([M], float)
+std_ar = np.zeros([M], float)
+std_gr = np.zeros([M], float)
+for imo in range (0, M):
+    # anom spec emis from 23GHz
+    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]))
+    if (len(anom_spec0_vec[nonzero(isnan(anom_spec0_vec) == True)]) != 0):
+        anom_spec0_zone[imo] = NaN
+        std_as0[imo] = NaN
+        anom_spec3_zone[imo] = NaN
+        std_as3[imo] = NaN
+        anom_ratio_zone[imo] = NaN
+        std_ar[imo] = NaN
+        grad_ratio_zone[imo] = NaN
+        std_gr[imo] = NaN
+        continue
+    else:
+        anom_spec0_zone[imo] = anom_spec0_vec.mean()
+        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))
+        # anom spec emis from 89GHz
+        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]))
+        anom_spec3_zone[imo] = anom_spec3_vec.mean()
+        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))
+        # anom emis ratio from 89GHz
+        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]))
+        anom_ratio_zone[imo] = anom_ratio_vec.mean()
+        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))
+        # gradient ratio from 23 and 89GHz
+        grad_ratio_vec = reshape(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1], size(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1]))
+        grad_ratio_zone[imo] = grad_ratio_vec.mean()
+        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))
+
+
+figure()
+fontP = FontProperties()
+fontP.set_size('small')
+subplot(2,1,1)
+errorbar(arange(0, M, 1), anom_spec0_zone, yerr = std_as0, label = 'emis spec anomaly 23GHz')
+errorbar(arange(0, M, 1), anom_spec3_zone, yerr = std_as3, label = 'emis spec anomaly 89GHz')
+errorbar(arange(0, M, 1), grad_ratio_zone, yerr = std_gr, label = 'gradient ratio 23-89GHz')
+plot(arange(0, M, 1), np.zeros([M], float), '--k')
+legend(loc = 3, prop = fontP)
+grid()
+xticks(arange(0, M, 1), month, rotation = 15)
+xlim(-1, M)
+subplot(2,1,2)
+errorbar(arange(0, M, 1), anom_ratio_zone, yerr = std_ar, label = 'emis ratio anomaly 89GHz')
+plot(arange(0, M, 1), np.zeros([M], float), '--k')
+legend(loc = 2, prop = fontP)
+grid()
+xticks(arange(0, M, 1), month, rotation = 15)
+xlim(-1, M)
+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')
+### map the selected zone ###
+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')
+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')
+'''
+
+
+
+'''
 ############################################################################
 # calculate correlation between anom ratio and gradient ratio (spec emiss) #
 ############################################################################
+corr_mat = np.zeros([M, 4, 4], float)
+for imo in range (0, M):
+    a = reshape(anom_spec[0, imo, :, :], size(anom_spec[0, imo, :, :]))[nonzero(isnan(reshape(anom_spec[0, imo, :, :], size(anom_spec[0, imo, :, :]))) == False)]
+    b = reshape(anom_spec[3, imo, :, :], size(anom_spec[3, imo, :, :]))[nonzero(isnan(reshape(anom_spec[3, imo, :, :], size(anom_spec[3, imo, :, :]))) == False)]
+    c = reshape(anom_ratio[3, imo, :, :], size(anom_ratio[3, imo, :, :]))[nonzero(isnan(reshape(anom_ratio[3, imo, :, :], size(anom_ratio[3, imo, :, :]))) == False)]
+    d = reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))[nonzero(isnan(reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))) == False)]
+    params = np.array([a, b, c, d])
+    for ii in range (0, 4):
+        for jj in range (0, 4):
+            corr_mat[imo, ii, jj] = corrcoef(params[ii], params[jj])[0][1]
+
+figure()
+fontP = FontProperties()
+fontP.set_size('small')
+plot(arange(0, M, 1), corr_mat[:, 0, 1], 'r', label = 'spec anomaly 23GHz - spec anomaly 89GHz')
+plot(arange(0, M, 1), corr_mat[:, 0, 2], 'g', label = 'spec anomaly 23GHz - ratio anomaly 89GHz')
+plot(arange(0, M, 1), corr_mat[:, 0, 3], 'b', label = 'spec anomaly 23GHz - gradient ratio')
+plot(arange(0, M, 1), corr_mat[:, 1, 2], 'c', label = 'spec anomaly 89GHz - ratio anomaly 89GHz')
+plot(arange(0, M, 1), corr_mat[:, 1, 3], 'm', label = 'spec anomaly 89GHz - gradient ratio')
+plot(arange(0, M, 1), corr_mat[:, 2, 3], 'y', label = 'ratio anomaly 89GHz - gradient ratio')
+plot(arange(0, M, 1), zeros([M], float), '--k')
+xlim(-1, M)
+ylim(-2, 1.)
+legend(loc = 8, prop = fontP)
+xticks(arange(0, M, 1), month, rotation = 25)
+yticks(np.arange(-2., 1., 0.1))
+grid()
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/correlation_params_2009.png')
+'''
+
+'''
 corr_mat = np.zeros([4, ny, nx], float)
 for ifr in range (0, 4):
@@ -117 +341 @@
     title('AMSUA ' + str(frequ[ifr]) + 'GHz - year 2009')
     plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/correlation_map_grad-ratio_emis-spec-anom_AMSUA' + str(frequ[ifr]) + '_with_AMSUA23-and-30-spec_year2009.png')
-
+'''
+
+
+
+
+
+
+
+

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

@@ -9 +9 @@
 from netCDF4 import Dataset
 import scipy.special
-import ffgrid2
-import map_ffgrid
 from matplotlib import colors
-import cartesian_grid_monthly
 import arctic_map
 import map_cartesian_grid
@@ -44 +41 @@
 
 
-
-###################
-# Read AMSUA data #
-###################
-es_a_month = np.zeros([M, ny, nx], float)
-el_a_month = np.zeros([M, ny, nx], float)
-esl_a_month = np.zeros([M, ny, nx], float)
-for imo in range (0, M):
-    print 'open file ' + str(month[imo])
-    file_amsua = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA30_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC')
-    lon_amsua = file_amsua.variables['longitude'][:]
-    lat_amsua = file_amsua.variables['latitude'][:]
-    days = file_amsua.variables['days'][:]
-    es_a = file_amsua.variables['e_spec'][:]
-    el_a = file_amsua.variables['e_lamb'][:]
-    esl_a = file_amsua.variables['e_spec_lamb'][:]
-    file_amsua.close()
-    print 'compute monthly mean of data for map'
-    for ilon in range (0, len(lon_amsua)):
-        for ilat in range (0, len(lat_amsua)):
-            es_a_month[imo, ilat, ilon] = mean(es_a[ilat, ilon, :][nonzero(isnan(es_a[ilat, ilon, :]) == False)])
-            el_a_month[imo, ilat, ilon] = mean(el_a[ilat, ilon, :][nonzero(isnan(el_a[ilat, ilon, :]) == False)]) 
-            esl_a_month[imo, ilat, ilon] = mean(esl_a[ilat, ilon, :][nonzero(isnan(esl_a[ilat, ilon, :]) == False)])
-
-
-print 'map of data'
-ion()
-x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
-x_coast = x_ind
-y_coast = y_ind
-z_coast = z_ind
-colormap = cm.s3pcpn_l_r
-for imo in range (0, M):
-    # emiss SPEC
-    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, es_a_month[imo, :, :], 0.45, 1.02, 0.01, colormap, 'emissivity SPEC')
-    title(month[imo] + ' 2009 - AMSUA 30GHz')
-    savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec_AMSUA30_' + month[imo] + '2009.png')
-    # emis SPEC-LAMB
-    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, esl_a_month[imo, :, :], 0., 0.03, 0.001, colormap, 'emissivity LAMB - SPEC')
-    title(month[imo] + ' 2009 - AMSUA 30GHz')
-    savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec-lamb_AMSUA30_' + month[imo] + '2009.png')
+frequ = np.array(['23', '30', '50', '89'])
+for ifr in range (3, 4):
+    print 'frequency ' + str(frequ[ifr]) + 'GHz'
+    ###################
+    # Read AMSUA data #
+    ###################
+    es_month = np.zeros([M, ny, nx], float)
+    el_month = np.zeros([M, ny, nx], float)
+    esl_month = np.zeros([M, ny, nx], float)
+    std_es_month = np.zeros([M, ny, nx], float)
+    std_el_month = np.zeros([M, ny, nx], float)
+    std_esl_month = np.zeros([M, ny, nx], float)
+    for imo in range (0, M):
+        print 'open file ' + str(month[imo])
+        fichier = 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[ifr]) + '_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC')
+        lon = fichier.variables['longitude'][:]
+        lat = fichier.variables['latitude'][:]
+        days = fichier.variables['days'][:]
+        es = fichier.variables['e_spec'][:]
+        el = fichier.variables['e_lamb'][:]
+        esl = fichier.variables['e_spec_lamb'][:]
+        fichier.close()
+        print 'compute monthly mean and std of data for mapping'
+        for ilon in range (0, len(lon)):
+            for ilat in range (0, len(lat)):
+                es_month[imo, ilat, ilon] = mean(es[ilat, ilon, :][nonzero(isnan(es[ilat, ilon, :]) == False)])
+                el_month[imo, ilat, ilon] = mean(el[ilat, ilon, :][nonzero(isnan(el[ilat, ilon, :]) == False)]) 
+                esl_month[imo, ilat, ilon] = mean(esl[ilat, ilon, :][nonzero(isnan(esl[ilat, ilon, :]) == False)])
+                std_es_month[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1)) * sum((es[ilat, ilon, :][nonzero(isnan(es[ilat, ilon, :]) == False)] - es_month[imo, ilat, ilon])**2))
+                std_el_month[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1)) * sum((el[ilat, ilon, :][nonzero(isnan(el[ilat, ilon, :]) == False)] - el_month[imo, ilat, ilon])**2))
+                std_esl_month[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1)) * sum((esl[ilat, ilon, :][nonzero(isnan(esl[ilat, ilon, :]) == False)] - esl_month[imo, ilat, ilon])**2))
+    ##################################
+    # stack std data in netcdf files #
+    ##################################
+    print 'start stacking'
+    for imo in range (0, M):
+        print 'stack in file month ' + str(month[imo])
+        rootgrp = 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_AMSUB' + str(frequ[ifr]) + '_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC')
+        rootgrp.createDimension('longitude', nx)
+        rootgrp.createDimension('latitude', ny)
+        nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
+        nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
+        nc_mean_spec = rootgrp.createVariable('mean_spec', 'f', ('latitude', 'longitude'))
+        nc_mean_lamb = rootgrp.createVariable('mean_lamb', 'f', ('latitude', 'longitude'))
+        nc_mean_lamb_spec = rootgrp.createVariable('mean_lamb_spec', 'f', ('latitude', 'longitude'))
+        nc_std_spec = rootgrp.createVariable('std_spec', 'f', ('latitude', 'longitude'))
+        nc_std_lamb = rootgrp.createVariable('std_lamb', 'f', ('latitude', 'longitude'))
+        nc_std_lamb_spec = rootgrp.createVariable('std_lamb_spec', 'f', ('latitude', 'longitude'))
+        nc_lon[:] = xvec
+        nc_lat[:] = yvec
+        nc_mean_spec[:] = es_month[imo, :, :]
+        nc_mean_lamb[:] = el_month[imo, :, :]
+        nc_mean_lamb_spec[:] = esl_month[imo, :, :]
+        nc_std_spec[:] = std_es_month[imo, :, :]
+        nc_std_lamb[:] = std_el_month[imo, :, :]
+        nc_std_lamb_spec[:] = std_esl_month[imo, :, :]
+        rootgrp.close()
+    '''
+    print 'map of data'
+    ion()
+    x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
+    x_coast = x_ind
+    y_coast = y_ind
+    z_coast = z_ind
+    colormap = cm.s3pcpn_l_r
+    for imo in range (0, M):
+        print 'map month ' + month[imo]
+        # emiss SPEC
+        print 'emis spec'
+        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, es_a_month[imo, :, :], 0.45, 1.02, 0.01, colormap, 'emissivity SPEC')
+        title(month[imo] + ' 2009 - AMSUA 30GHz')
+        savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec_AMSUA' + str(frequ[ifr]) + '_' + month[imo] + '2009.png')
+        # emis SPEC-LAMB
+        print 'emis lamb-spec'
+        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, esl_a_month[imo, :, :], 0., 0.03, 0.001, colormap, 'emissivity LAMB - SPEC')
+        title(month[imo] + ' 2009 - AMSUA 30GHz')
+        savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec-lamb_AMSUA' + str(frequ[ifr]) + '_' + month[imo] + '2009.png')
+        # std SPEC
+        print 'std emis spec'
+        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_es_month[imo, :, :], 0., std_es_month[nonzero(isnan(std_es_month) == False)].max(), 0.001, colormap, 'std of emissivity SPEC')
+        title(month[imo] + ' 2009 - AMSUA ' + str(frequ[ifr]) + 'GHz')
+        savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/std_emis_spec_AMSUA' + str(frequ[ifr]) + '_' + month[imo] + '2009.png')
+        # std LAMB
+        print 'std emis lamb'
+        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_el_month[imo, :, :], 0., std_el_month[nonzero(isnan(std_el_month) == False)].max(), 0.001, colormap, 'std of emissivity LAMB')
+        title(month[imo] + ' 2009 - AMSUA ' + str(frequ[ifr]) + 'GHz')
+        savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/std_emis_lamb_AMSUA' + str(frequ[ifr]) + '_' + month[imo] + '2009.png')
+        # std LAMB-SPEC
+        print 'std emis lamb-spec'
+        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_esl_month[imo, :, :], 0., std_esl_month[nonzero(isnan(std_esl_month) == False)].max(), 0.0001, colormap, 'std emissivity LAMB - SPEC')
+        title(month[imo] + ' 2009 - AMSUA ' + str(frequ[ifr]) + 'GHz')
+        savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/std_emis_lamb-spec_AMSUA' + str(frequ[ifr]) + '_' + month[imo] + '2009.png')
+        '''
 
 
 
-