Changeset 52

Timestamp:

08/11/14 18:51:39 (10 years ago)

Author:

lahlod

Message:

modifs

Location:

trunk/src/scripts_Laura/ARCTIC/Travail_CEN

Files:

• daily_emis_AMSUA_AMSUB_89.py (modified) (12 diffs)
• daily_evolution_anomaly_grad_ratio.py (modified) (4 diffs)

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

@@ -32 +32 @@
 x0 = -3000. # min limit of grid
 x1 = 2500. # max limit of grid
-dx = 40.
-xvec = np.arange(x0, x1+dx, dx)
-nx = len(xvec) 
+# AMSUA grid
+dx_a = 40.
+xvec_a = np.arange(x0, x1+dx_a, dx_a)
+nx_a = len(xvec_a) 
+# AMSUB grid
+dx_b = 20.
+xvec_b = np.arange(x0, x1+dx_b, dx_b)
+nx_b = len(xvec_b) 
+
 y0 = -3000. # min limit of grid
 y1 = 3000. # max limit of grid
-dy = 40.
-yvec = np.arange(y0, y1+dy, dy)
-ny = len(yvec)
-
-
-
-'''
+# AMSUA grid
+dy_a = 40.
+yvec_a = np.arange(y0, y1+dy_a, dy_a)
+ny_a = len(yvec_a)
+# AMSUB grid
+dy_b = 20.
+yvec_b = np.arange(y0, y1+dy_b, dy_b)
+ny_b = len(yvec_b)
+
+
 
 ############################################
@@ -53 +62 @@
 
 # select borders of zone
-yi = 1880. 
-yf = 2280.
-xi = -480.
-xf = -120.
-
-#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]
-
-len(xvec[xxi:xxf+1])
-len(yvec[yyi:yyf+1])
+yi = 960. 
+yf = 1360.
+xi = -680.
+xf = -320.
+
+#find corresponding index in xvec and yvec...
+# ... in AMSUA grid
+xxi_a = np.where(xvec_a == xi)[0][0]
+xxf_a = np.where(xvec_a == xf)[0][0]
+yyi_a = np.where(yvec_a == yi)[0][0]
+yyf_a = np.where(yvec_a == yf)[0][0]
+len(xvec_a[xxi_a:xxf_a+1])
+len(yvec_a[yyi_a:yyf_a+1])
+
+# ... in AMSUB grid
+xxi_b = np.where(xvec_b == xi)[0][0]
+xxf_b = np.where(xvec_b == xf)[0][0]
+yyi_b = np.where(yvec_b == yi)[0][0]
+yyf_b = np.where(yvec_b == yf)[0][0]
+len(xvec_b[xxi_b:xxf_b+1])
+len(yvec_b[yyi_b:yyf_b+1])
+
+
+
+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
+map_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')
 
 
@@ -71 +97 @@
 
 mean_spec_a_zone = np.zeros([M, 31], float)
+mean_spec_a_zone23 = np.zeros([M, 31], float)
 std_spec_a_zone = np.zeros([M, 31], float)
 
@@ -85 +112 @@
     print 'month ' + month[imo]
     print 'read daily emis spec lamb and diff AMSUA and AMSUB'
-    # AMSUA
+    # AMSUA 89
     fichier_a = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
     spec_a = fichier_a.variables['e_spec'][:]
@@ -91 +118 @@
     diff_a = fichier_a.variables['e_spec_lamb'][:]
     fichier_a.close()
+    # AMSUA 23
+    fichier_a23 = 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')
+    spec_a23 = fichier_a23.variables['e_spec'][:]
+    lamb_a23 = fichier_a23.variables['e_lamb'][:]
+    diff_a23 = fichier_a23.variables['e_spec_lamb'][:]
+    fichier_a23.close()
     # AMSUB
-    fichier_b = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
+    fichier_b = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_20/cartesian_grid_res-20_monthly_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
     spec_b = fichier_b.variables['e_spec'][:]
     lamb_b = fichier_b.variables['e_lamb'][:]
     diff_b = fichier_b.variables['e_spec_lamb'][:]
     fichier_b.close()
-    spec_a_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float)
-    lamb_a_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float)
-    spec_b_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float)
-    lamb_b_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float)
+    spec_a_vec = np.zeros([month_day[imo], len(xvec_a[xxi_a : xxf_a+1]) * len(yvec_a[yyi_a : yyf_a+1])], float)
+    spec_a_vec23 = np.zeros([month_day[imo], len(xvec_a[xxi_a : xxf_a+1]) * len(yvec_a[yyi_a : yyf_a+1])], float)
+    lamb_a_vec = np.zeros([month_day[imo], len(xvec_a[xxi_a : xxf_a+1]) * len(yvec_a[yyi_a : yyf_a+1])], float)
+    spec_b_vec = np.zeros([month_day[imo], len(xvec_b[xxi_b : xxf_b+1]) * len(yvec_b[yyi_b : yyf_b+1])], float)
+    lamb_b_vec = np.zeros([month_day[imo], len(xvec_b[xxi_b : xxf_b+1]) * len(yvec_b[yyi_b : yyf_b+1])], float)
     # calculate daily mean and std zonal spec and lamb AMSUA and AMSUB
     print 'calculate daily mean and std'
@@ -106 +140 @@
         print 'date' + str(ijr + 1) + ' ' + month[imo] + ' 2009'
         # spec AMSUA
-        spec_a_vec[ijr, :] = reshape(spec_a[yyi : yyf + 1, xxi : xxf + 1, ijr], size(spec_a[yyi : yyf + 1, xxi : xxf + 1, ijr]))
+        print 'spec amsua'
+        spec_a_vec[ijr, :] = reshape(spec_a[yyi_a : yyf_a + 1, xxi_a : xxf_a + 1, ijr], size(spec_a[yyi_a : yyf_a + 1, xxi_a : xxf_a + 1, ijr]))
+        spec_a_vec23[ijr, :] = reshape(spec_a23[yyi_a : yyf_a + 1, xxi_a : xxf_a + 1, ijr], size(spec_a23[yyi_a : yyf_a + 1, xxi_a : xxf_a + 1, ijr]))
         mean_spec_a_zone[imo, ijr] = mean(spec_a_vec[ijr, :][nonzero(isnan(spec_a_vec[ijr, :]) == False)])
-        std_spec_a_zone[imo, ijr] = sqrt((1. / (size(spec_a[yyi : yyf+1, xxi : xxf+1, ijr]) - 1.)) * sum((spec_a_vec[ijr, :][nonzero(isnan(spec_a_vec[ijr, :]) == False)] - mean_spec_a_zone[imo, ijr])**2))
+        mean_spec_a_zone23[imo, ijr] = mean(spec_a_vec23[ijr, :][nonzero(isnan(spec_a_vec23[ijr, :]) == False)])
+        std_spec_a_zone[imo, ijr] = sqrt((1. / (size(spec_a[yyi_a : yyf_a+1, xxi_a : xxf_a+1, ijr]) - 1.)) * sum((spec_a_vec[ijr, :][nonzero(isnan(spec_a_vec[ijr, :]) == False)] - mean_spec_a_zone[imo, ijr])**2))
         # lamb AMSUA
-        lamb_a_vec[ijr, :] = reshape(lamb_a[yyi : yyf + 1, xxi : xxf + 1, ijr], size(lamb_a[yyi : yyf + 1, xxi : xxf + 1, ijr])) 
+        print 'lamb amsua'
+        lamb_a_vec[ijr, :] = reshape(lamb_a[yyi_a : yyf_a + 1, xxi_a : xxf_a + 1, ijr], size(lamb_a[yyi_a : yyf_a + 1, xxi_a : xxf_a + 1, ijr])) 
         mean_lamb_a_zone[imo, ijr] = mean(lamb_a_vec[ijr, :][nonzero(isnan(lamb_a_vec[ijr, :]) == False)])
-        std_lamb_a_zone[imo, ijr] = sqrt((1. / (size(lamb_a[yyi : yyf+1, xxi : xxf+1, ijr]) - 1.)) * sum((lamb_a_vec[ijr, :][nonzero(isnan(lamb_a_vec[ijr, :]) == False)] - mean_lamb_a_zone[imo, ijr])**2))
+        std_lamb_a_zone[imo, ijr] = sqrt((1. / (size(lamb_a[yyi_a : yyf_a+1, xxi_a : xxf_a+1, ijr]) - 1.)) * sum((lamb_a_vec[ijr, :][nonzero(isnan(lamb_a_vec[ijr, :]) == False)] - mean_lamb_a_zone[imo, ijr])**2))
         # spec AMSUB
-        spec_b_vec[ijr, :] = reshape(spec_b[yyi : yyf + 1, xxi : xxf + 1, ijr], size(spec_b[yyi : yyf + 1, xxi : xxf + 1, ijr]))
+        print 'spec amsub'
+        spec_b_vec[ijr, :] = reshape(spec_b[yyi_b : yyf_b + 1, xxi_b : xxf_b + 1, ijr], size(spec_b[yyi_b : yyf_b + 1, xxi_b : xxf_b + 1, ijr]))
         mean_spec_b_zone[imo, ijr] = mean(spec_b_vec[ijr, :][nonzero(isnan(spec_b_vec[ijr, :]) == False)])
-        std_spec_b_zone[imo, ijr] = sqrt((1. / (size(spec_b[yyi : yyf+1, xxi : xxf+1, ijr]) - 1.)) * sum((spec_b_vec[ijr, :][nonzero(isnan(spec_b_vec[ijr, :]) == False)] - mean_spec_b_zone[imo, ijr])**2))
+        std_spec_b_zone[imo, ijr] = sqrt((1. / (size(spec_b[yyi_b : yyf_b+1, xxi_b : xxf_b+1, ijr]) - 1.)) * sum((spec_b_vec[ijr, :][nonzero(isnan(spec_b_vec[ijr, :]) == False)] - mean_spec_b_zone[imo, ijr])**2))
         # lamb AMSUB
-        lamb_b_vec[ijr, :] = reshape(lamb_b[yyi : yyf + 1, xxi : xxf + 1, ijr], size(lamb_b[yyi : yyf + 1, xxi : xxf + 1, ijr])) 
+        print 'lamb amsub'
+        lamb_b_vec[ijr, :] = reshape(lamb_b[yyi_b : yyf_b + 1, xxi_b : xxf_b + 1, ijr], size(lamb_b[yyi_b : yyf_b + 1, xxi_b : xxf_b + 1, ijr])) 
         mean_lamb_b_zone[imo, ijr] = mean(lamb_b_vec[ijr, :][nonzero(isnan(lamb_b_vec[ijr, :]) == False)])
-        std_lamb_b_zone[imo, ijr] = sqrt((1. / (size(lamb_b[yyi : yyf+1, xxi : xxf+1, ijr]) - 1.)) * sum((lamb_b_vec[ijr, :][nonzero(isnan(lamb_b_vec[ijr, :]) == False)] - mean_lamb_b_zone[imo, ijr])**2))
+        std_lamb_b_zone[imo, ijr] = sqrt((1. / (size(lamb_b[yyi_b : yyf_b+1, xxi_b : xxf_b+1, ijr]) - 1.)) * sum((lamb_b_vec[ijr, :][nonzero(isnan(lamb_b_vec[ijr, :]) == False)] - mean_lamb_b_zone[imo, ijr])**2))
         if (isnan(mean_spec_a_zone[imo, ijr]) == True):
             std_spec_a_zone[imo, ijr] = NaN
@@ -132 +172 @@
 
 mean_year_spec_a = np.append(mean_spec_a_zone[0, 0 : month_day[0]], mean_spec_a_zone[1, 0 : month_day[1]])
+mean_year_spec_a23 = np.append(mean_spec_a_zone23[0, 0 : month_day[0]], mean_spec_a_zone23[1, 0 : month_day[1]])
 std_year_spec_a = np.append(std_spec_a_zone[0, 0 : month_day[0]], std_spec_a_zone[1, 0 : month_day[1]])
 
@@ -146 +187 @@
     # spec_a
     mean_year_spec_a = np.append(mean_year_spec_a, mean_spec_a_zone[imo, 0 : month_day[imo]])
+    mean_year_spec_a23 = np.append(mean_year_spec_a23, mean_spec_a_zone23[imo, 0 : month_day[imo]])
     std_year_spec_a = np.append(std_year_spec_a, std_spec_a_zone[imo, 0 : month_day[imo]])
     # lamb_a
@@ -164 +206 @@
 vec_months = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334])
 ion()
+figure()
+plot(mean_year_spec_a, '+-r', label = 'AMSUA89')
+plot(mean_year_spec_a23, '+-m', label = 'AMSUA23')
+plot(mean_year_spec_b, '+-g', label = 'AMSUB89')
+xlim(0, 365)
+ylim(0.5, 1.)
+xticks(vec_months, month, rotation = 25)
+yticks(np.arange(0.5, 1., 0.05))
+fontP = FontProperties()
+fontP.set_size('small')
+legend(loc = 3, prop = fontP)
+grid()
+ylabel('emissivity spec')
+title('Beaufort Sea')
+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_AMSUA89_AMSUB89_zone_Beaufort_Sea.png')
+################################
+# plot mean difference and std #
+################################
 figure()
 plot(mean_year_spec_a - mean_year_spec_b, 'r', label = 'mean spec AMSUA - mean spec AMSUB]')
@@ -194 +254 @@
 # monthly mean and std in whole arctic # 
 ########################################
-bias_spec = np.zeros([M, ny, nx], float)
+bias_spec = np.zeros([M, ny_a, nx_a], float)
+bias_anom = np.zeros([M, ny_a, nx_a], float)
+bias_mean = np.zeros([M], float)
 for imo in range (0, M):
     print 'month ' + month[imo]
     print 'read daily emis spec lamb and diff AMSUA and AMSUB'
-    # AMSUA
-    fichier_a = 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')
-    spec_a = fichier_a.variables['mean_spec'][:]
-    lamb_a = fichier_a.variables['mean_lamb'][:]
-    std_spec_a = fichier_a.variables['std_spec'][:]
-    std_lamb_a = fichier_a.variables['std_spec'][:]
-    fichier_a.close()
+    # AMSUA 89
+    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')
+    spec_a89 = fichier_a89.variables['mean_spec'][:]
+    lamb_a89 = fichier_a89.variables['mean_lamb'][:]
+    std_spec_a89 = fichier_a89.variables['std_spec'][:]
+    std_lamb_a89 = fichier_a89.variables['std_spec'][:]
+    fichier_a89.close()
+    # AMSUA 23
+    fichier_a23 = 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')
+    spec_a23 = fichier_a23.variables['mean_spec'][:]
+    lamb_a23 = fichier_a23.variables['mean_lamb'][:]
+    std_spec_a23 = fichier_a23.variables['std_spec'][:]
+    std_lamb_a23 = fichier_a23.variables['std_spec'][:]
+    fichier_a23.close()
     # AMSUB
-    fichier_b = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_20/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
+    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')
     spec_b = fichier_b.variables['mean_spec'][:]
     lamb_b = fichier_b.variables['mean_lamb'][:]
@@ -213 +282 @@
     fichier_b.close()
     bias_spec[imo, :, :] = spec_a - spec_b
+    bias_mean[imo] = mean(bias_spec[imo, :, :][nonzero(isnan(bias_spec[imo, :, :]) == False)])
+    for ilon in range (0, nx_a):
+        for ilat in range (0, ny_a):
+            bias_anom[imo, ilat, ilon] = bias_spec[imo, ilat, ilon] - bias_mean[imo]
+   
 
 
@@ -225 +299 @@
 z_coast = z_ind
 for imo in range (0, M):
-    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, bias_spec[imo, :, :], -0.02, 0.02, 0.001, cm.s3pcpn_l_r, 'Bias of emis spec AMSUA89-AMSUB89')
+    print 'map for month ' + month[imo]
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec_a, yvec_a, bias_anom[imo, :, :], -0.018, 0.012, 0.001, cm.s3pcpn_l_r, 'Bias anomaly of emis spec AMSUA89-AMSUB89')
     title(month[imo] + ' 2009')
-    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/map_bias_AMSUA89-AMSUB89_arctic_' + month[imo] + '2009.png')
-
+    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/map_bias_anomaly_AMSUA89-AMSUB89_arctic_' + month[imo] + '2009.png')
+

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

@@ -43 +43 @@
 
 
-
+'''
 ############################################
 # time evolution (monthly) in a given zone #
@@ -221 +221 @@
 
 
-'''
+
 subplot(2, 1, 2)
 plot(std_year_grad_ratio, '--b', label = 'std grad ratio')
@@ -230 +230 @@
 legend(loc = 2, prop = fontP)
 xlim(0, 365)
-'''
+
 
 
@@ -244 +244 @@
 title('area of study - Chukchi Sea')
 plt.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')
-
-
-
-map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, gr[15, :, :], gr[15, :, :][nonzero(isnan(gr[15, :, :]) == False)].min(), gr[15, :, :][nonzero(isnan(gr[15, :, :]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'daily gradient ratio (01-01-2009)')
-contour(xvec[xxi:xxf+1], yvec[yyi:yyf+1], )
+'''
+
+
+################################
+# map parameters in all Arctic #
+################################
+cumul_params = np.zeros([M, ny, nx], float)
+grad_ratio = np.zeros([M, ny, nx], float)
+spec_anom_23 = np.zeros([M, ny, nx], float)
+spec_anom_89 = np.zeros([M, ny, nx], float)
+ratio_anom_89 = np.zeros([M, ny, nx], float)
+CP = np.zeros([M, ny, nx], float)
+for imo in range (0, M): 
+    # daily read gradient ratio file
+    print 'read daily gradient ratio for month ' + month[imo]
+    fichier_grad_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_grad_ratio_spec_23-89_near_nadir_AMSUA_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
+    gr = fichier_grad_ratio.variables['grad_ratio'][:]
+    fichier_grad_ratio.close()
+    # read daily emis anomaly file for 23GHz
+    print 'read daily emis anomaly 23GHz for month ' + month[imo]
+    fichier_anom23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
+    sa_23 = fichier_anom23.variables['spec_anomaly'][:]
+    fichier_anom23.close()
+    # read daily emis anomaly file for 89GHz
+    print 'read daily emis anomaly 89GHz for month ' + month[imo]
+    fichier_anom89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
+    sa_89 = fichier_anom89.variables['spec_anomaly'][:]
+    ra_89 = fichier_anom89.variables['ratio_anomaly'][:]
+    fichier_anom89.close()
+    for ilon in range (0, nx):
+        for ilat in range (0, ny):
+            grad_ratio[imo, ilat, ilon] = mean(gr[:, ilat, ilon][nonzero(isnan(gr[:, ilat, ilon]) == False)])
+            spec_anom_23[imo, ilat, ilon] = mean(sa_23[:, ilat, ilon][nonzero(isnan(sa_23[:, ilat, ilon]) == False)])
+            spec_anom_89[imo, ilat, ilon] = mean(sa_89[:, ilat, ilon][nonzero(isnan(sa_89[:, ilat, ilon]) == False)])
+            ratio_anom_89[imo, ilat, ilon] = mean(ra_89[:, ilat, ilon][nonzero(isnan(ra_89[:, ilat, ilon]) == False)])
+            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])
+    CP[imo, :, :] = (cumul_params[imo, :, :]/max(cumul_params[imo, :, :][nonzero(isnan(cumul_params[imo, :, :]) == False)])) * 100.
+
+
+#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):
+    print 'map month ' + month[imo]
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, CP[imo, :, :], 0., 100., 1., cm.s3pcpn_l_r, 'Cumulation index (%)')
+    title('AMSUA - ' + month[imo] + ' 2009')
+    savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/cumul_params/cumul_all_parameters/map_cumulation_index_'+ str(MONTH[imo]) + month[imo] + '2009.png')