Changeset 48

Timestamp:

08/01/14 18:56:25 (10 years ago)

Author:

lahlod

Message:

new scripts

Location:

trunk/src/scripts_Laura/ARCTIC/Travail_CEN

Files:

• choose_new_classified_points.py (modified) (2 diffs)
• compare_ice_area_different_data.py (modified) (1 diff)
• read_sub_classification_cross_correlation.py (added)
• regrid_distrib.py (added)

trunk/src/scripts_Laura/ARCTIC/Travail_CEN/choose_new_classified_points.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
+from matplotlib import colors
 from matplotlib.font_manager import FontProperties
 import map_cartesian_grid
-import ice_class_delimit_AMSU_data
-
-
-
-
-
-
-
-
-fichier1 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB89_SEPTEMBER2009.nc', 'r', format='NETCDF3_CLASSIC')
-x89 = fichier1.variables['longitude'][:]
-y89 = fichier1.variables['latitude'][:]
-day89 = fichier1.variables['days'][:]
-emis_spec_89 = fichier1.variables['e_spec'][:]
-emis_lamb_89 = fichier1.variables['e_lamb'][:]
-fichier1.close()
-
-
-
-fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_SEPTEMBER2009_AMSUA30_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC')
-xdist = fichier_emis.variables['longitude'][:]
-ydist = fichier_emis.variables['latitude'][:]
-emis_spec = fichier_emis.variables['spec_ice_area'][:]
-emis_lamb = fichier_emis.variables['lamb_ice_area'][:]
-fichier_emis.close()
-spec_89 = np.zeros([len(ydist), len(xdist)], float)
-for ilon in range (0, len(xdist)):
-    for ilat in range (0, len(ydist)):
-        if (isnan(emis_spec[ilat, ilon]) == True):
-            spec_89[ilat, ilon] = NaN
-        else:
-            spec_89[ilat, ilon] = emis_spec_89[ilat, ilon, 0]
-
-
-
-
+
+
+###############################
+# time period characteristics #
+###############################
+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)
+
+
+########################
+# grid characteristics #
+########################
+x0 = -3000. # min limit of grid
+x1 = 2500. # max limit of grid
+dx = 40.
+xvec = np.arange(x0, x1+dx, dx)
+nx = len(xvec) 
+y0 = -3000. # min limit of grid
+y1 = 3000. # max limit of grid
+dy = 40.
+yvec = np.arange(y0, y1+dy, dy)
+ny = len(yvec)
+
+
+##################################################################################################################
+# We devide the loop in two steps : 
+# - first loop concerns all months except for AUGUST and SEPTEMBER - use of AMSUA23GHz SPEC emissivity to seperate ice from no-ice zones 
+# - second loop concerns AUGUST and SEPTEMBER - use of AMSUA30GHz SPEC emissivity to seperate ice from no_ice zones
+##################################################################################################################
+frequ = 89 # apply threshold on this frequency
+# open .dat file to stack data (see end of loop)
+#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')
+bin = 50
+
+
+# monthly mean parameter (2D-array) on ARCTIC area
+spec_month = np.zeros([M, ny, nx], float)
+lamb_month = np.zeros([M, ny, nx], float)
+diff_month = np.zeros([M, ny, nx], float)
+ratio_month = np.zeros([M, ny, nx], float)
+# monthly mean parameter (2D-array) on ARCTIC SEA ICE area
+spec_ice = np.zeros([M, ny, nx], float)
+lamb_ice = np.zeros([M, ny, nx], float)
+diff_ice = np.zeros([M, ny, nx], float)
+ratio_ice = np.zeros([M, ny, nx], float)
+# monthly mean parameter (1D-array) on ARCTIC SEA ICE area transformed into vector
+spec_vec = np.zeros([M, ny * nx], float)
+lamb_vec = np.zeros([M, ny * nx], float)
+diff_vec = np.zeros([M, ny * nx], float)
+ratio_vec = np.zeros([M, ny * nx], float)
+# histogram distribution (intensity of occurence) of parameter in SEA ICE area (1D-array, bins = 200)
+hist_vals_spec = np.zeros([M, bin], float)
+hist_vals_lamb = np.zeros([M, bin], float)
+hist_vals_diff = np.zeros([M, bin], float)
+hist_vals_ratio = np.zeros([M, bin], float)
+# histogram distribution (emissivity value) of parameter in SEA ICE area (1D-array, bins = 200)
+corresp_emis_spec = np.zeros([M, bin], float)
+corresp_emis_lamb = np.zeros([M, bin], float)
+corresp_emis_diff = np.zeros([M, bin], float)
+corresp_emis_ratio = np.zeros([M, bin], float)
+months1 = np.array([0, 1, 2, 3, 4, 5, 6, 9, 10, 11]) # use AMSUA 23GHz to delimit ice/no_ice for all months except for AUGUST and SEPTEMBER
+for imo in months1:
+    print 'month ' + month[imo]
+    ##################################################################################
+    # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study #
+    ##################################################################################
+    print 'open threshold file'
+    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')
+    spec_lim = fichier_emis.variables['spec_ice_area'][:]
+    #lamb_lim = fichier_emis.variables['lamb_ice_area'][:]
+    fichier_emis.close()
+    #########################################################
+    # application of AMSUA 23GHz delimitation to other data #
+    #########################################################
+    print 'open emissivity to sub_classify file'
+    ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB)
+    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')
+    day = fichier_e.variables['days'][:]
+    emis_spec = fichier_e.variables['e_spec'][:]
+    emis_lamb = fichier_e.variables['e_lamb'][:]
+    emis_diff = fichier_e.variables['e_spec_lamb'][:]
+    fichier_e.close()
+    for ilon in range (0, nx):
+        for ilat in range (0, ny):
+            spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, :][nonzero(isnan(emis_spec[ilat, ilon, :]) == False)])
+            lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, :][nonzero(isnan(emis_lamb[ilat, ilon, :]) == False)])
+            diff_month[imo, ilat, ilon] = mean(emis_diff[ilat, ilon, :][nonzero(isnan(emis_diff[ilat, ilon, :]) == False)])
+    ## data of emis RATIO
+    fichier_r = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
+    ratio_month[imo, :, :] = fichier_r.variables['emis_ratio'][:]
+    fichier_r.close()
+    print 'compute matrix of parameter on SEA ICE area'
+    for ilon in range (0, nx):
+        for ilat in range (0, ny):
+            if (isnan(spec_lim[ilat, ilon]) == True):
+                spec_ice[imo, ilat, ilon] = NaN
+                lamb_ice[imo, ilat, ilon] = NaN
+                diff_ice[imo, ilat, ilon] = NaN
+                ratio_ice[imo, ilat, ilon] = NaN
+            else:
+                spec_ice[imo, ilat, ilon] = spec_month[imo, ilat, ilon]
+                lamb_ice[imo, ilat, ilon] = lamb_month[imo, ilat, ilon]
+                diff_ice[imo, ilat, ilon] = diff_month[imo, ilat, ilon]
+                ratio_ice[imo, ilat, ilon] = ratio_month[imo, ilat, ilon]
+    print 'compute SPEC distribution'
+    ########
+    # SPEC #
+    ########
+    cs = reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))[nonzero(isnan(reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))) == False)]
+    spec_vec[imo, 0 : len(cs)] = cs
+    hist_vals_spec[imo, :] = hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_spec[imo, ibin] = mean(hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    print 'compute LAMB distribution'
+    ########
+    # LAMB #
+    ########
+    cl = reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))[nonzero(isnan(reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))) == False)]
+    lamb_vec[imo, 0 : len(cl)] = cl
+    hist_vals_lamb[imo, :] = hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_lamb[imo, ibin] = mean(hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    print 'compute DIFF distribution'
+    ########
+    # DIFF #
+    ########
+    cd = reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))[nonzero(isnan(reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))) == False)]
+    diff_vec[imo, 0 : len(cd)] = cd
+    hist_vals_diff[imo, :] = hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_diff[imo, ibin] = mean(hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    print 'compute RATIO distribution'
+    #########
+    # RATIO #
+    #########
+    cr = reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))[nonzero(isnan(reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))) == False)]
+    ratio_vec[imo, 0 : len(cr)] = cr
+    hist_vals_ratio[imo, :] = hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_ratio[imo, ibin] = mean(hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    ######################
+    # stack in .dat file #
+    ######################
+    print 'start stacking in .dat file'
+    #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')
+    '''for ii in range (0, bin):
+        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' % {
+                                            'months':month[imo],
+                                            'hist_vals_spec':hist_vals_spec[imo, ii],
+                                            'corresp_emis_spec':corresp_emis_spec[imo, ii],
+                                            'hist_vals_lamb':hist_vals_lamb[imo, ii],
+                                            'corresp_emis_lamb':corresp_emis_lamb[imo, ii],
+                                            'hist_vals_diff':hist_vals_diff[imo, ii],
+                                            'corresp_emis_diff':corresp_emis_diff[imo, ii],
+                                            'hist_vals_rate':hist_vals_ratio[imo, ii],
+                                            'corresp_emis_rate':corresp_emis_ratio[imo, ii],
+                                            }))'''
+    ########################
+    # stack in netcdf file #
+    ######################## 
+    print 'stack in file month ' + str(month[imo])
+    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-30_' + month[imo] + '2009_AMSUA' + str(frequ) + '_spec_lamb_thresholds.nc', 'w', format='NETCDF3_CLASSIC')
+    rootgrp.createDimension('longitude', len(xvec))
+    rootgrp.createDimension('latitude', len(yvec))
+    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
+    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
+    nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude'))
+    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude'))
+    nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude'))
+    nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude'))
+    nc_lon[:] = xvec
+    nc_lat[:] = yvec
+    nc_ice_spec[:] = spec_ice[imo, :, :]
+    nc_ice_lamb[:] = lamb_ice[imo, :, :]
+    nc_ice_diff[:] = diff_ice[imo, :, :]
+    nc_ice_ratio[:] = ratio_ice[imo, :, :]
+    rootgrp.close()
+
+
+months2 = np.array([7, 8])# use AMSUA 30GHz to delimit ice/no_ice for AUGUST and SEPTEMBER
+for imo in months2:
+    print 'month ' + month[imo]
+    ##################################################################################
+    # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study #
+    ##################################################################################
+    print 'open threshold file'
+    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')
+    spec_lim = fichier_emis.variables['spec_ice_area'][:]
+    #lamb_lim = fichier_emis.variables['lamb_ice_area'][:]
+    fichier_emis.close()
+    #########################################################
+    # application of AMSUA 23GHz delimitation to other data #
+    #########################################################
+    print 'open emissivity to sub_classify file'
+    ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB)
+    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')
+    day = fichier_e.variables['days'][:]
+    emis_spec = fichier_e.variables['e_spec'][:]
+    emis_lamb = fichier_e.variables['e_lamb'][:]
+    emis_diff = fichier_e.variables['e_spec_lamb'][:]
+    fichier_e.close()
+    for ilon in range (0, nx):
+        for ilat in range (0, ny):
+            spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, :][nonzero(isnan(emis_spec[ilat, ilon, :]) == False)])
+            lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, :][nonzero(isnan(emis_lamb[ilat, ilon, :]) == False)])
+            diff_month[imo, ilat, ilon] = mean(emis_diff[ilat, ilon, :][nonzero(isnan(emis_diff[ilat, ilon, :]) == False)])
+    ## data of emis RATIO
+    fichier_r = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
+    ratio_month[imo, :, :] = fichier_r.variables['emis_ratio'][:]
+    fichier_r.close()
+    print 'compute matrix of parameter on SEA ICE area'
+    for ilon in range (0, nx):
+        for ilat in range (0, ny):
+            if (isnan(spec_lim[ilat, ilon]) == True):
+                spec_ice[imo, ilat, ilon] = NaN
+                lamb_ice[imo, ilat, ilon] = NaN
+                diff_ice[imo, ilat, ilon] = NaN
+                ratio_ice[imo, ilat, ilon] = NaN
+            else:
+                spec_ice[imo, ilat, ilon] = spec_month[imo, ilat, ilon]
+                lamb_ice[imo, ilat, ilon] = lamb_month[imo, ilat, ilon]
+                diff_ice[imo, ilat, ilon] = diff_month[imo, ilat, ilon]
+                ratio_ice[imo, ilat, ilon] = ratio_month[imo, ilat, ilon]
+    print 'compute SPEC distribution'
+    ########
+    # SPEC #
+    ########
+    cs = reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))[nonzero(isnan(reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))) == False)]
+    spec_vec[imo, 0 : len(cs)] = cs
+    hist_vals_spec[imo, :] = hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_spec[imo, ibin] = mean(hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    print 'compute LAMB distribution'
+    ########
+    # LAMB #
+    ########
+    cl = reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))[nonzero(isnan(reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))) == False)]
+    lamb_vec[imo, 0 : len(cl)] = cl
+    hist_vals_lamb[imo, :] = hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_lamb[imo, ibin] = mean(hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    print 'compute DIFF distribution'
+    ########
+    # DIFF #
+    ########
+    cd = reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))[nonzero(isnan(reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))) == False)]
+    diff_vec[imo, 0 : len(cd)] = cd
+    hist_vals_diff[imo, :] = hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_diff[imo, ibin] = mean(hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    print 'compute RATIO distribution'
+    #########
+    # RATIO #
+    #########
+    cr = reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))[nonzero(isnan(reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))) == False)]
+    ratio_vec[imo, 0 : len(cr)] = cr
+    hist_vals_ratio[imo, :] = hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[0]
+    for ibin in range (0, bin):
+        corresp_emis_ratio[imo, ibin] = mean(hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])
+    ######################
+    # stack in .dat file #
+    ######################
+    print 'start stacking in .dat file'
+    #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')
+    '''for ii in range (0, bin):
+        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' % {
+                                            'months':month[imo],
+                                            'hist_vals_spec':hist_vals_spec[imo, ii],
+                                            'corresp_emis_spec':corresp_emis_spec[imo, ii],
+                                            'hist_vals_lamb':hist_vals_lamb[imo, ii],
+                                            'corresp_emis_lamb':corresp_emis_lamb[imo, ii],
+                                            'hist_vals_diff':hist_vals_diff[imo, ii],
+                                            'corresp_emis_diff':corresp_emis_diff[imo, ii],
+                                            'hist_vals_rate':hist_vals_ratio[imo, ii],
+                                            'corresp_emis_rate':corresp_emis_ratio[imo, ii],
+                                            }))'''
+    ########################
+    # stack in netcdf file #
+    ######################## 
+    print 'stack in file month ' + str(month[imo])
+    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-30_' + month[imo] + '2009_AMSUA' + str(frequ) + '_spec_lamb_thresholds.nc', 'w', format='NETCDF3_CLASSIC')
+    rootgrp.createDimension('longitude', len(xvec))
+    rootgrp.createDimension('latitude', len(yvec))
+    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
+    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
+    nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude'))
+    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude'))
+    nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude'))
+    nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude'))
+    nc_lon[:] = xvec
+    nc_lat[:] = yvec
+    nc_ice_spec[:] = spec_ice[imo, :, :]
+    nc_ice_lamb[:] = lamb_ice[imo, :, :]
+    nc_ice_diff[:] = diff_ice[imo, :, :]
+    nc_ice_ratio[:] = ratio_ice[imo, :, :]
+    rootgrp.close()
+
+'''
+data_classif.close()
+'''
+
+
+
+'''
 # test:
 ion()
@@ -52 +328 @@
 y_coast = y_ind
 z_coast = z_ind
-map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, emis_spec[:, :], 0.45, 1.02, 0.01, cm.s3pcpn_l_r, 'test')
+for imo in range (0, M):
+    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ratio_ice[imo, :, :], 0., 4., 0.01, cm.s3pcpn_l_r, 'Sea ice extent - emissivity RATIO')
+    title('AMSUA ' + str(frequ) + ' - ' + 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/emiss_ratio_map_AMSUA'+str(frequ)+'_with_AMSUA23-and-30-spec_'+str(month[imo])+'_2009.png')
+'''
+
+

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

@@ -164 +164 @@
 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]
+    a[imo] = (AS[0, imo] - area_osi[imo]) / area_osi[imo]
+    b[imo] = (AL[0, imo] - area_osi[imo]) / area_osi[imo]
+    c[imo] = (AS[1, imo] - area_osi[imo]) / area_osi[imo]
+    d[imo] = (AL[1, imo] - area_osi[imo]) / area_osi[imo]
+    e[imo] = (AS[2, imo] - area_osi[imo]) / area_osi[imo]
+    f[imo] = (AL[2, imo] - area_osi[imo]) / area_osi[imo]
+    g[imo] = (AS[3, imo] - area_osi[imo]) / area_osi[imo]
+    h[imo] = (AL[3, imo] - area_osi[imo]) / area_osi[imo]
+    i[imo] = (area_s_B[imo] - area_osi[imo]) / area_osi[imo]
+    j[imo] = (area_l_B[imo] - area_osi[imo]) / area_osi[imo]