Changeset 47

Timestamp:

07/30/14 18:21:31 (10 years ago)

Author:

lahlod

Message:

new_scripts

Location:

trunk/src/scripts_Laura/ARCTIC/Travail_CEN

Files:

• choose_new_classified_points.py (added)
• compare_ice_area_different_data.py (modified) (2 diffs)
• ice_class_delimit_AMSU_data.py (modified) (5 diffs)
• ice_delimit_AMSUB_data.py (modified) (21 diffs)
• import_ice_class_delimit_AMSU_data.py (modified) (5 diffs)
• test.py (added)

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

@@ -146 +146 @@
 xlim(-1, M)
 grid()
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/compar_total_ice_area_AMSUA_AMSUB_SPEC_LAMB_OSISAF_2009.png')
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/total_ice_area/compar_total_ice_area_AMSUA_AMSUB_SPEC_LAMB_OSISAF_2009.png')
 
 
@@ -153 +153 @@
 # calculation of bias and std #
 ###############################
+a = np.zeros([M], float)
+b = np.zeros([M], float)
+c = np.zeros([M], float)
+d = np.zeros([M], float)
+e = np.zeros([M], float)
+f = np.zeros([M], float)
+g = np.zeros([M], float)
+h = np.zeros([M], float)
+i = np.zeros([M], float)
+j = np.zeros([M], float)
+for imo in range (0, M):
+    a[imo] = AS[0, imo] - area_osi[imo]
+    b[imo] = AL[0, imo] - area_osi[imo]
+    c[imo] = AS[1, imo] - area_osi[imo]
+    d[imo] = AL[1, imo] - area_osi[imo]
+    e[imo] = AS[2, imo] - area_osi[imo]
+    f[imo] = AL[2, imo] - area_osi[imo]
+    g[imo] = AS[3, imo] - area_osi[imo]
+    h[imo] = AL[3, imo] - area_osi[imo]
+    i[imo] = area_s_B[imo] - area_osi[imo]
+    j[imo] = area_l_B[imo] - area_osi[imo]
+
+
 figure()
-plot((AS[0, :] - AS[1, :]) * 100. / ((mean(AS[0, :]) + mean(AS[0, :])) / 2.), 'r', label = 'AMSUA spec 23GHz - 50GHz')
-plot((AS[0, :] - AS[2, :]) * 100. / ((mean(AS[0, :]) + mean(AS[2, :])) / 2.), 'b', label = 'AMSUA spec 23GHz - 89GHz')
-plot((AS[1, :] - AS[2, :]) * 100. / ((mean(AS[1, :]) + mean(AS[2, :])) / 2.), 'g', label = 'AMSUA spec 50GHz - 89GHz')
+plot(a, 'c', label = 'AMSUA spec 23GHz')
+plot(b, 'c--', label = 'AMSUA lamb 23GHz')
+plot(c, 'r', label = 'AMSUA spec 30GHz')
+plot(d, 'r--', label = 'AMSUA lamb 30GHz')
+plot(e, 'm', label = 'AMSUA spec 50GHz')
+plot(f, 'm--', label = 'AMSUA lamb 50GHz')
+plot(g, 'g', label = 'AMSUA spec 89GHz')
+plot(h, 'g--', label = 'AMSUA lamb 89GHz')
+plot(i, 'b', label = 'AMSUB spec 89GHz')
+plot(j, 'b--', label = 'AMSUB lamb 89GHz')
+plot(np.arange(0, M, 1), np.zeros([M], float), 'k')
 fontP = FontProperties()
 fontP.set_size('small')
 legend(loc = 3, prop = fontP)
-ylabel('bias of total ice area (%)')
-yticks(np.arange(0, M, 1), month, rotation = 25)
+ylabel('bias of total ice area')
 xticks(np.arange(0, M, 1), month, rotation = 25)
 xlim(-1, M)
 grid()
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/total_ice_area/compar_bias_total_ice_area_AMSUA_SPEC_2009.png')
-
-
-
-
-
+plt.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')
+
+
+
+
+

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

@@ -85 +85 @@
     emis_ratio = np.zeros([ny, nx, M], float)
     for imo in range (0, M):
-        #print 'month ' + month[imo]
+        print 'month ' + month[imo]
         ##############
         # emissivity #
         ##############
-        #print 'read file for emiss'
+        print 'read file for emiss'
         fichier_emis = 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) + '_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC')
         xdist = fichier_emis.variables['longitude'][:]
@@ -97 +97 @@
         emis_lamb = fichier_emis.variables['e_lamb'][:]
         fichier_emis.close()
-        #print 'calculation of monthly mean'
+        print 'calculation of monthly mean'
         for ilon in range (0, nx):
             for ilat in range (0, ny):
                 emis_spec_moy[ilat, ilon, imo] = mean(emis_spec[ilat, ilon, :][nonzero(isnan(emis_spec[ilat, ilon, :]) == False)])
                 emis_lamb_moy[ilat, ilon, imo] = mean(emis_lamb[ilat, ilon, :][nonzero(isnan(emis_lamb[ilat, ilon, :]) == False)])
-        #print 'calculation of monthly mean difference lamb-spec'
+        print 'calculation of monthly mean difference lamb-spec'
         emis_diff[:, :, imo] = emis_lamb_moy[:, :, imo] - emis_spec_moy[:, :, imo]
         #########
         # ratio #
         #########
-        #print 'read file for ratio'
+        print 'read file for ratio'
         fichier_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA' + str(frequ) + '_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC')
         xdist = fichier_ratio.variables['longitude'][:]
@@ -119 +119 @@
     # emissivity distribution after filtering #
     ###########################################
-    '''
     hist_val_spec = np.zeros([50, M], float)
     hist_val_lamb = np.zeros([50, M], float)
@@ -127 +126 @@
     corresp_val_lamb = np.zeros([50, M], float)
     corresp_val_ratio = np.zeros([50, M], float)
-    corresp_val_diff = np.zeros([50, M], float)
-    '''
+    corresp_val_diff = np.zeros([50, M], float)   
     emis_spec_f = np.zeros([7000, M], float)
     emis_lamb_f = np.zeros([7000, M], float)
     emis_ratio_f = np.zeros([7000, M], float)
     emis_diff_f = np.zeros([7000, M], float)
+    X = np.zeros([7000, M], float)
+    Y = np.zeros([7000, M], float)
     L_spec = np.zeros([M], int)
     for imo in range (0, M):
-        #print 'month ' + month[imo]
+        print 'month ' + month[imo]
         # choice of spec emissivity as the threshold for the study // definition of x and y index corresponding to the points which emissivity value is over threshold
         y_index_spec = np.where(emis_spec_moy[:, :, imo] >= lim_spec[imo])[0] 
+        Y[0 : len(y_index_spec), imo] = y_index_spec[:]
         x_index_spec = np.where(emis_spec_moy[:, :, imo] >= lim_spec[imo])[1]
+        X[0 : len(x_index_spec), imo] = x_index_spec[:]
         L_spec[imo] = len(x_index_spec) # = len(y_index)
         #print 'length of x and y index ', L_spec[imo]
         # definition of filtered values (vectors) with the previous threshold // values of filtered emissivity SPEC, LAMB, rate and difference LAMB-SPEC 
+        print 'calculation of vectors of filtered data'
         for ii in range (0, L_spec[imo]):
             emis_spec_f[ii, imo] = emis_spec_moy[y_index_spec[ii], x_index_spec[ii], imo]
@@ -147 +150 @@
             emis_ratio_f[ii, imo] = emis_ratio[y_index_spec[ii], x_index_spec[ii], imo]
             emis_diff_f[ii, imo] = emis_diff[y_index_spec[ii], x_index_spec[ii], imo]
-        '''
+        print 'calculation of histogram distribution'
         # definition of their distribution within the new filtered values
-        hist_val_spec[:, imo] = hist(emis_spec_f, bins = 50, normed = True, histtype='step')[0]
-        hist_val_lamb[:, imo] = hist(emis_lamb_f, bins = 50, normed = True, histtype='step')[0]
-        hist_val_ratio[:, imo] = hist(emis_ratio_f, bins = 50, normed = True, histtype='step')[0]
-        hist_val_diff[:, imo] = hist(emis_diff_f, bins = 50, normed = True, histtype='step')[0]
+        hist_val_spec[:, imo] = hist(emis_spec_f[0 : L_spec[imo], imo], bins = 50, normed = True, histtype='step')[0]
+        hist_val_lamb[:, imo] = hist(emis_lamb_f[0 : L_spec[imo], imo], bins = 50, normed = True, histtype='step')[0]
+        hist_val_ratio[:, imo] = hist(emis_ratio_f[0 : L_spec[imo], imo], bins = 50, normed = True, histtype='step')[0]
+        hist_val_diff[:, imo] = hist(emis_diff_f[0 : L_spec[imo]
+
+, imo], bins = 50, normed = True, histtype='step')[0]
         for ibin in range (0, 50):
-            corresp_val_spec[ibin, imo] = mean(hist(emis_spec_f, bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
-            corresp_val_lamb[ibin, imo] = mean(hist(emis_lamb_f, bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
-            corresp_val_ratio[ibin, imo] = mean(hist(emis_ratio_f, bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
-            corresp_val_diff[ibin, imo] = mean(hist(emis_diff_f, bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
-        '''
+            corresp_val_spec[ibin, imo] = mean(hist(emis_spec_f[0 : L_spec[imo], imo], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
+            corresp_val_lamb[ibin, imo] = mean(hist(emis_lamb_f[0 : L_spec[imo], imo], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
+            corresp_val_ratio[ibin, imo] = mean(hist(emis_ratio_f[0 : L_spec[imo], imo], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
+            corresp_val_diff[ibin, imo] = mean(hist(emis_diff_f[0 : L_spec[imo], imo], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2])
 
-    return(emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec)      
+
+    return(emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec, X, Y, hist_val_spec, hist_val_lamb, hist_val_ratio, hist_val_diff, corresp_val_spec, corresp_val_lamb, corresp_val_ratio, corresp_val_diff)      
 
 

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

@@ -41 +41 @@
 
 
+frequ = 89
 
 #####################
@@ -56 +57 @@
     ### emissivity ###
     print 'open emissivity file'
-    fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA30_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
+    fichier_emis = 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')
     xdist = fichier_emis.variables['longitude'][:]
     ydist = fichier_emis.variables['latitude'][:]
@@ -70 +71 @@
     ### monthly emissivity ratio ###
     print 'open emissivity ratio file'
-    fichier_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA30_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
+    fichier_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUB' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
     xdist = fichier_ratio.variables['longitude'][:]
     ydist = fichier_ratio.variables['latitude'][:]
@@ -94 +95 @@
 
 c = np.array(['r', 'b', 'c', 'm', 'y', 'g'])
-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
-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
-idata = 1
+#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
+#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
+#idata = 0
+fontP = FontProperties()
+fontP.set_size('small')
 print 'distribution of emissivity values'
 ###################################
@@ -119 +122 @@
 xlabel('emissivity SPEC')
 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_AMSUA30_JANUARY-JUNE_2009.png')
+#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()
@@ -133 +134 @@
 ylabel('frequency of occurence')
 legend(loc = 9, prop = fontP)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_spec_AMSUA30_JULY-DECEMBER_2009.png')
+#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')
 
 # find the emissivity value corresponding to the threshold of ice/no_ice limit
@@ -139 +140 @@
 for imo in range (0, M):
     print 'month ' + str(month[imo])
-    bb = np.where((corresp_emis_spec[imo, :] < limit_coef_spec[idata]))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY
-    aa = np.where(hist_vals_spec[imo, bb] < max(hist_vals_spec[imo, bb]) / 2.)[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity
-    cc = np.where(hist_vals_spec[imo, bb] == max(hist_vals_spec[imo, bb]))[0] # which emissivity index corresponds to peak emissivity
-    dd = np.where(aa > cc[0])[0]
-    emis_lim_spec[imo] = corresp_emis_spec[imo, bb][aa][dd[0]]
+    bb = np.where((corresp_emis_spec[imo, :] > 0.7) & (corresp_emis_spec[imo, :] < 0.77))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY
+    aa = np.where(hist_vals_spec[imo, bb] == min(hist_vals_spec[imo, bb]))[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity
+    #cc = np.where(hist_vals_spec[imo, bb] == max(hist_vals_spec[imo, bb]))[0] # which emissivity index corresponds to peak emissivity
+    #dd = np.where(aa > cc[0])[0]
+    emis_lim_spec[imo] = corresp_emis_spec[imo, bb][aa][0]
 
 # plot monthly evolution of emissivity spec threshold used for delimitation
@@ -153 +154 @@
 legend(loc = 2)
 xlim(-1, M)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emis_lim_frequ_SPEC_AMSUA30_2009.png')
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUB/emis_lim_frequ_SPEC_AMSUB' + str(frequ) + '_2009.png')
 
 
@@ -177 +178 @@
 ylabel('frequency of occurence')
 legend(loc = 2, prop = fontP)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_lamb_AMSUA30_JANUARY-JUNE_2009.png')
+#plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_lamb_AMSUA' + str(frequ) + '_JANUARY-JUNE_2009.png')
 ## plot six following months of spec emissivity histograms ##
 figure()
@@ -188 +189 @@
 ylabel('frequency of occurence')
 legend(loc = 9, prop = fontP)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_lamb_AMSUA30_JULY-DECEMBER_2009.png')
+#plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_lamb_AMSUA' + str(frequ) + '_JULY-DECEMBER_2009.png')
 
 # find the emissivity value corresponding to the threshold of ice/no_ice limit
@@ -194 +195 @@
 for imo in range (0, M):
     print 'month ' + str(month[imo])
-    bb = np.where((corresp_emis_lamb[imo, :] < limit_coef_lamb[idata]))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY
-    aa = np.where(hist_vals_lamb[imo, bb] < max(hist_vals_lamb[imo, bb]) / 2.)[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity
-    cc = np.where(hist_vals_lamb[imo, bb] == max(hist_vals_lamb[imo, bb]))[0] # which emissivity index corresponds to peak emissivity
-    dd = np.where(aa > cc[0])[0]
-    emis_lim_lamb[imo] = corresp_emis_lamb[imo, bb][aa][dd[0]]
+    bb = np.where((corresp_emis_lamb[imo, :] > 0.72) & (corresp_emis_lamb[imo, :] < 0.8))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY
+    aa = np.where(hist_vals_lamb[imo, bb] == min(hist_vals_lamb[imo, bb]))[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity
+    #cc = np.where(hist_vals_lamb[imo, bb] == max(hist_vals_lamb[imo, bb]))[0] # which emissivity index corresponds to peak emissivity
+    #dd = np.where(aa > cc[0])[0]
+    emis_lim_lamb[imo] = corresp_emis_lamb[imo, bb][aa][0]
 
 # plot monthly evolution of emissivity spec threshold used for delimitation
@@ -208 +209 @@
 legend(loc = 2)
 xlim(-1, M)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emis_lim_frequ_LAMB_AMSUA30_2009.png')
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUB/emis_lim_frequ_LAMB_AMSUB' + str(frequ) + '_2009.png')
 
 
@@ -233 +234 @@
 ylabel('frequency of occurence')
 legend(prop = fontP)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_rate_AMSUA30_JANUARY-JUNE_2009.png')
+#plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_rate_AMSUA'+str(frequ)+'_JANUARY-JUNE_2009.png')
 ## plot six following months of spec emissivity histograms ##
 figure()
@@ -244 +245 @@
 ylabel('frequency of occurence')
 legend(loc = 9, prop = fontP)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_rate_AMSUA30_JULY-DECEMBER_2009.png')
+#plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_rate_AMSUA'+str(frequ)+'_JULY-DECEMBER_2009.png')
 
 # no definition of threshold for ICE / NO_ICE delimitation because no apparent distinction signal // use of this signal for internal ice classification
@@ -280 +281 @@
 legend(loc = 2, prop = fontP)
 xlim(-1, M)
-plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emis_lim_frequ_spec_and_lamb_AMSUA30_2009.png')
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUB/emis_lim_frequ_spec_and_lamb_AMSUB' + str(frequ) + '_2009.png')
 
 
@@ -304 +305 @@
             else:
                 if (emis_spec_month[imo, ilat, ilon] <= emis_lim_spec[imo]): # use the monthly SPEC emissivity threshold 
-                    ice_zone_spec[imo, ilat, ilon] = 0.
+                    ice_zone_spec[imo, ilat, ilon] = NaN
                 else:
-                    ice_zone_spec[imo, ilat, ilon] = 1.   
-    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_spec[imo, :, :], -1, 2., 1., colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_SPEC > ' + str(emis_lim_spec[imo])[0:6] + ')')
-    title(str(month[imo]) + ' 2009 - AMSUA 30GHz')
-    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/AMSUA30_ice_emis_spec_thresh' + '_' + month[imo] + '2009.png')
+                    ice_zone_spec[imo, ilat, ilon] = emis_spec_month[imo, ilat, ilon]   
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_spec[8, :, :], 0.5, 1.02, 0.01,  cm.s3pcpn_l_r, 'Sea ice emissivity spec (threshold : emis_SPEC > ' + str(emis_lim_spec[imo])[0:6] + ')')
+    title(str(month[imo]) + ' 2009 - AMSUA ' + str(frequ) + 'GHz')
+    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/AMSUA' + str(frequ) + '_ice_emis_spec_thresh' + '_' + month[imo] + '2009.png')
 
 #### LAMB emiss ####
@@ -322 +323 @@
             else:
                 if (emis_lamb_month[imo, ilat, ilon] <= emis_lim_lamb[imo]): # use the monthly SPEC emissivity threshold 
-                    ice_zone_lamb[imo, ilat, ilon] = 0.
+                    ice_zone_lamb[imo, ilat, ilon] = NaN
                 else:
-                    ice_zone_lamb[imo, ilat, ilon] = 1.   
-    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_lamb[imo, :, :], -1, 2., 1., colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_LAMB > ' + str(emis_lim_lamb[imo])[0:6] + ')')
-    title(str(month[imo]) + ' 2009 - AMSUA 30GHz')
-    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/AMSUA30_ice_emis_lamb_thresh' + '_' + month[imo] + '2009.png')
+                    ice_zone_lamb[imo, ilat, ilon] = emis_lamb_month[imo, ilat, ilon]     
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone_lamb[8, :, :], 0.5, 1.02, 0.01,  cm.s3pcpn_l_r, 'Sea ice emissivity lamb (threshold : emis_LAMB > ' + str(emis_lim_lamb[imo])[0:6] + ')')
+    title(str(month[imo]) + ' 2009 - AMSUA ' + str(frequ) + 'GHz')
+    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/AMSUA' + str(frequ) + '_ice_emis_lamb_thresh' + '_' + month[imo] + '2009.png')
 
 
@@ -343 +344 @@
 for imo in range (0, M):
     ice_spec = reshape(ice_zone_spec[imo, :, :], size(ice_zone_spec[imo, :, :]))
-    nb_pix_spec[imo] = len(np.where(ice_spec == 1.)[0])
+    nb_pix_spec[imo] = len(np.where(isnan(ice_spec) == False)[0])
     total_ice_area_spec[imo] = (pix_area * nb_pix_spec[imo]) + (926. * pix_area)
 
@@ -353 +354 @@
 for imo in range (0, M):
     ice_lamb = reshape(ice_zone_lamb[imo, :, :], size(ice_zone_lamb[imo, :, :]))
-    nb_pix_lamb[imo] = len(np.where(ice_lamb == 1.)[0])
+    nb_pix_lamb[imo] = len(np.where(isnan(ice_lamb) == False)[0])
     total_ice_area_lamb[imo] = (pix_area * nb_pix_lamb[imo]) + (926. * pix_area)
 
@@ -364 +365 @@
 ######################
 print 'start stacking in .dat file'
-data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/AMSUA30_data_classification_parameters_ice_no-ice_2009.dat', 'a')
+data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/AMSUB'+str(frequ)+'_data_classification_parameters_ice_no-ice_2009.dat', 'a')
 for imo in range (0, M):
     for ii in range (0, 50):
@@ -391 +392 @@
 for imo in range (0, M):
     print 'stack in file month ' + str(month[imo])
-    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + month[imo] + '2009_AMSUA30_spec_lamb_thresholds.nc', 'w', format='NETCDF3_CLASSIC')
+    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/cartesian_grid_map_ice_no-ice_' + month[imo] + '2009_AMSUB' + str(frequ) + '_spec_lamb_thresholds.nc', 'w', format='NETCDF3_CLASSIC')
     rootgrp.createDimension('longitude', len(xdist))
     rootgrp.createDimension('latitude', len(ydist))

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

@@ -5 +5 @@
 import matplotlib.pyplot as plt
 from pylab import *
-from mpl_toolkits.basemap import Basemap
-from mpl_toolkits.basemap import shiftgrid, cm
 from netCDF4 import Dataset
-import arctic_map # function to regrid coast limits
-import cartesian_grid_test # function to convert grid from polar to cartesian
-import scipy.special
-import ffgrid2
-import map_ffgrid
+#import arctic_map # function to regrid coast limits
+#import cartesian_grid_test # function to convert grid from polar to cartesian
 from matplotlib.font_manager import FontProperties
-import map_cartesian_grid
+#import map_cartesian_grid
 import ice_class_delimit_AMSU_data
 
 
 
-frequ = 23
-emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec = ice_class_delimit_AMSU_data.filtering(frequ)
-spec_23 = emis_spec_f
-lamb_23 = emis_lamb_f
-ratio_23 = emis_ratio_f
-diff_23 = emis_diff_f
-L_spec_23 = L_spec
-
+
+MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'])
+month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER'])
+month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
+M = len(month)
+
+
+
+frequ = 89
+
+################################################################################
+# compute filtered points of emissivity SPEC, LAMB, rate, difference LAMB-SPEC #
+################################################################################
+emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec, X, Y, hist_val_spec, hist_val_lamb, hist_val_ratio, hist_val_diff, corresp_val_spec, corresp_val_lamb, corresp_val_ratio, corresp_val_diff = ice_class_delimit_AMSU_data.filtering(frequ)
+spec = emis_spec_f
+lamb = emis_lamb_f
+ratio = emis_ratio_f
+diff = emis_diff_f
+L_spec = L_spec
+X = X
+Y = Y
+hist_spec = hist_val_spec
+hist_lamb = hist_val_lamb
+hist_ratio = hist_val_ratio
+hist_diff = hist_val_diff
+corresp_spec = corresp_val_spec
+corresp_lamb = corresp_val_lamb
+corresp_ratio = corresp_val_ratio
+corresp_diff = corresp_val_diff
+
+
+
+########
+# plot #
+########
+ion()
+c = np.array(['r', 'b', 'c', 'm', 'y', 'g'])
+fontP = FontProperties()
+fontP.set_size('small')
+#### SPEC ####
+figure()
+for imo in range (0, 6):
+    plot(corresp_spec[:, imo], hist_spec[:, imo], c = str(c[imo]), label = str(month[imo]))
+
+grid()
+xlim(corresp_spec.min() - 0.02, corresp_spec.max() + 0.02)
+xlabel('emissivity spec')
+ylabel('frequency of occurence')
+legend(prop = fontP, loc = 2)
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/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_spec[:, imo], hist_spec[:, imo], c = str(c[imo - 6]), label = str(month[imo]))
+
+grid()
+xlim(corresp_spec.min() - 0.02, corresp_spec.max() + 0.02)
+xlabel('emissivity spec')
+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/sub_classification/emiss_spec_AMSUA' + str(frequ) + '_JULY-DECEMBER_2009.png')
+
+
+
+#### LAMB ####
+figure()
+for imo in range (0, 6):
+    plot(corresp_lamb[:, imo], hist_lamb[:, imo], c = str(c[imo]), label = str(month[imo]))
+
+grid()
+xlim(corresp_lamb.min() - 0.02, corresp_lamb.max() + 0.02)
+xlabel('emissivity lamb')
+ylabel('frequency of occurence')
+fontP = FontProperties()
+fontP.set_size('small')
+legend(prop = fontP, loc = 2)
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/emiss_lamb_AMSUA' + str(frequ) + '_JANUARY-JUNE_2009.png')
+## plot six following months of spec emissivity histograms ##
+figure()
+for imo in range (6, M):
+    plot(corresp_lamb[:, imo], hist_lamb[:, imo], c = str(c[imo - 6]), label = str(month[imo]))
+
+grid()
+xlim(corresp_lamb.min() - 0.02, corresp_lamb.max() + 0.02)
+xlabel('emissivity lamb')
+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/sub_classification/emiss_lamb_AMSUA' + str(frequ) + '_JULY-DECEMBER_2009.png')
+
+#### RATE ####
+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() - 0.02, corresp_ratio.max() + 0.02)
+xlabel('emissivity ratio')
+ylabel('frequency of occurence')
+fontP = FontProperties()
+fontP.set_size('small')
+legend(prop = fontP, loc = 1)
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/emiss_ratio_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() - 0.02, corresp_ratio.max() + 0.02)
+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/sub_classification/emiss_ratio_AMSUA' + str(frequ) + '_JULY-DECEMBER_2009.png')
+
+#### DIFF ####
+figure()
+for imo in range (0, 6):
+    plot(corresp_diff[:, imo], hist_diff[:, imo], c = str(c[imo]), label = str(month[imo]))
+
+grid()
+xlim(corresp_diff.min() - 0.002, corresp_diff.max() + 0.002)
+xlabel('emissivity diff')
+ylabel('frequency of occurence')
+fontP = FontProperties()
+fontP.set_size('small')
+legend(prop = fontP, loc = 1)
+plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/emiss_diff_AMSUA' + str(frequ) + '_JANUARY-JUNE_2009.png')
+## plot six following months of spec emissivity histograms ##
+figure()
+for imo in range (6, M):
+    plot(corresp_diff[:, imo], hist_diff[:, imo], c = str(c[imo - 6]), label = str(month[imo]))
+
+grid()
+xlim(corresp_diff.min() - 0.002, corresp_diff.max() + 0.002)
+xlabel('emissivity diff')
+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/sub_classification/emiss_diff_AMSUA' + str(frequ) + '_JULY-DECEMBER_2009.png')
+
+
+
+
+'''
 frequ = 30
-emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec = ice_class_delimit_AMSU_data.filtering(frequ)
+emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec, X, Y, hist_val_spec, hist_val_lamb, hist_val_ratio, hist_val_diff, corresp_val_spec, corresp_val_lamb, corresp_val_ratio, corresp_val_diff = ice_class_delimit_AMSU_data.filtering(frequ)
 spec_30 = emis_spec_f
 lamb_30 = emis_lamb_f
@@ -34 +164 @@
 diff_30 = emis_diff_f
 L_spec_30 = L_spec
+X_30 = X
+Y_30 = Y
 
 frequ = 50
-emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec = ice_class_delimit_AMSU_data.filtering(frequ)
+emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec, X, Y, hist_val_spec, hist_val_lamb, hist_val_ratio, hist_val_diff, corresp_val_spec, corresp_val_lamb, corresp_val_ratio, corresp_val_diff = ice_class_delimit_AMSU_data.filtering(frequ)
 spec_50 = emis_spec_f
 lamb_50 = emis_lamb_f
@@ -42 +174 @@
 diff_50 = emis_diff_f
 L_spec_50 = L_spec
+X_50 = X
+Y_50 = Y
 
 frequ = 89
-emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec = ice_class_delimit_AMSU_data.filtering(frequ)
+emis_spec_f, emis_lamb_f, emis_ratio_f, emis_diff_f, L_spec, X, Y, hist_val_spec, hist_val_lamb, hist_val_ratio, hist_val_diff, corresp_val_spec, corresp_val_lamb, corresp_val_ratio, corresp_val_diff = ice_class_delimit_AMSU_data.filtering(frequ)
 spec_89 = emis_spec_f
 lamb_89 = emis_lamb_f
@@ -50 +184 @@
 diff_89 = emis_diff_f
 L_spec_89 = L_spec
-
-
-
-
-
-MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'])
-month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER'])
-month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
-M = len(month)
+X_89 = X
+Y_89 = Y
+
+
+XX = X_89[:, 0][nonzero(X_89[:, 0] != 0.)]
+YY = Y_89[:, 0][nonzero(Y_89[:, 0] != 0.)]
+L = len(XX)
+for ii in range (0, L):
+    emis_spec_moy[YY[ii], XX[ii]]
+
+
+
 
 a1 = np.zeros([M], float)
@@ -70 +207 @@
     a4[imo] = corrcoef(spec_89[0 : 3243, imo], diff_89[0 : 3243, imo])[0][1]
 
-correl_matrix = corrcoef(np.array([spec_89[0 : 3243, imo], lamb_89[0 : 3243, imo], ratio_89[0 : 3243, imo], diff_89[0 : 3243, imo], spec_50[0 : 3243, imo], lamb_50[0 : 3243, imo], ratio_50[0 : 3243, imo], diff_50[0 : 3243, imo], spec_30[0 : 3243, imo], lamb_30[0 : 3243, imo], ratio_30[0 : 3243, imo], diff_30[0 : 3243, imo], spec_23[0 : 3243, imo], lamb_23[0 : 3243, imo], ratio_23[0 : 3243, imo], diff_23[0 : 3243, imo]]))
+params = np.array([spec_89[0 : 3243, imo], lamb_89[0 : 3243, imo], ratio_89[0 : 3243, imo], diff_89[0 : 3243, imo], spec_50[0 : 3243, imo], lamb_50[0 : 3243, imo], ratio_50[0 : 3243, imo], diff_50[0 : 3243, imo], spec_30[0 : 3243, imo], lamb_30[0 : 3243, imo], ratio_30[0 : 3243, imo], diff_30[0 : 3243, imo], spec[0 : 3243, imo], lamb[0 : 3243, imo], ratio[0 : 3243, imo], diff[0 : 3243, imo]]))
 
 figure()
 pc = pcolor(correl_matrix)
 colorbar(pc)
-
+'''