- Timestamp:
- 08/22/14 15:06:02 (10 years ago)
- Location:
- trunk/src/scripts_Laura/ARCTIC/Travail_CEN
- Files:
-
- 5 added
- 4 deleted
- 20 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/arctic_map.py
r41 r56 11 11 12 12 13 ####################################################### 14 # Defines every point from bathymetry NETCDF file # 15 # that corresponds to altitude = 0 (coast) and stacks # 16 # these points in cartesian coordinate vectors x_ind # 17 # and y_ind # 18 ####################################################### 13 19 14 20 def arctic_map_lat50(): … … 40 46 z_ind = Z_LIM[indices] 41 47 42 volume = z_ind / 50. 48 volume = z_ind / 50. # size of the points defining the coast (with z_ind = 1.) 43 49 #plt.ion() 44 50 #fig = plt.figure() -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/cartesian_grid_test.py
r46 r56 32 32 # definition of the new cartesian grid (x, y) # 33 33 ############################################### 34 x0 = -3000. # min limit of grid35 x1 = 2500. # max limit of grid34 x0 = -3000. # min x limit of grid (in km) 35 x1 = 2500. # max x limit of grid (in km) 36 36 xvec = np.arange(x0, x1+dx, dx) 37 37 nx = len(xvec) 38 y0 = -3000. # min limit of grid39 y1 = 3000. # max limit of grid38 y0 = -3000. # min y limit of grid (in km) 39 y1 = 3000. # max y limit of grid (in km) 40 40 yvec = np.arange(y0, y1+dy, dy) 41 41 ny = len(yvec) 42 xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole )42 xgrid_cart, ygrid_cart= np.meshgrid(xvec, yvec) # new cartesian grid (centered on North pole [x, y] = [0, 0]) 43 43 44 44 … … 46 46 # calculation of the new gridding # 47 47 ################################### 48 zgrid_output = np.zeros([ny, nx, nb_days], float) 49 ngrid_output = np.zeros([ny, nx, nb_days], float) 50 z2grid_output = np.zeros([ny, nx, nb_days], float) 51 sigmagrid_output = np.zeros([ny, nx, nb_days], float) 52 bblat = nonzero(lati >= 65.) # only select high latitude values of z48 zgrid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of data in new cartesian grid 49 ngrid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of number of data in each grid cell of new cartesian grid 50 z2grid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of data*data in new cartesian grid 51 sigmagrid_output = np.zeros([ny, nx, nb_days], float) # 3D-array of std of data in each grid cell of new cartesian grid 52 bblat = nonzero(lati >= 65.) # only select high latitude observations (above 65° lat) 53 53 for ijr in range (0, nb_days): # loop on time - here time = days 54 54 bbjour = nonzero(jour[bblat] == ijr + 1.)[0] … … 61 61 # origin of the (x, y) grid : (x=0, y=0) <=> (lon=0, lat=0) 62 62 L = len(z) 63 Rt = 6371. 64 alpha = (pi*Rt)/180. 63 Rt = 6371. # radius of earth 64 alpha = (pi*Rt)/180. # convertion from deg to rad angles 65 65 beta = pi/180. 66 x = np.zeros([L], float) 67 y = np.zeros([L], float) 66 x = np.zeros([L], float) # x coordinates of new cartesian grid 67 y = np.zeros([L], float) # y coordinates of new cartesian grid 68 68 for k in range (0, L): 69 69 if ((longitude[k] >= 0.) & (longitude[k] <= 90.)): # 4eme quadrant … … 90 90 ix[i] = 0 91 91 else: 92 ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x va kue to a cell number92 ix[i] = math.ceil((x[i] - x0) / dx) - 1 # associates each x value to a cell number 93 93 i = 0 94 94 iy = np.zeros([L], int) … … 101 101 # calculation of distances between point M(x,y) and 4 grid points of its belonging cell # 102 102 ######################################################################################### 103 close_point = np.zeros([L, 2], int) 103 close_point = np.zeros([L, 2], int) # for x- and y-axis, 2D-array of closest grid cell of new cartesian grid to observation 104 104 k = 0 105 105 for k in range (0, L): # boucle sur les points M (x et y) … … 108 108 d3 = sqrt(((x[k] - xvec[ix[k]]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2)) 109 109 d4 = sqrt(((x[k] - xvec[ix[k] + 1]) ** 2) + ((y[k] - yvec[iy[k] + 1]) ** 2)) 110 d_vec = np.array([d1, d2, d3, d4]) 110 d_vec = np.array([d1, d2, d3, d4]) # vector of distances between observation and each closer grid_cell of new cartesian grid 111 111 ind = np.where(d_vec == min(d_vec)) # finds the smallest distance between the 4 points of the grid 112 112 point_vec = np.array([(ix[k], iy[k]), (ix[k] + 1, iy[k]), (ix[k], iy[k] + 1), (ix[k] + 1, iy[k] + 1)]) -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/choose_new_classified_points.py
r54 r56 47 47 # - second loop concerns MAY, JUNE, JULY, AUGUST - use of AMSUA89GHz SPEC emissivity to seperate ice from no_ice zones 48 48 ################################################################################################################## 49 50 49 51 frequ = 89 # apply threshold on this frequency 50 '''51 #open .dat file to stack data (see end of loop)52 data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/AMSUA'+str(frequ)+'_data_classification_parameters_ice_no-ice_with_AMSUA23-and-30-spec_2009.dat', 'a')53 bin = 5054 '''55 52 56 53 … … 67 64 daily_ratio_ice = np.zeros([M, ny, nx, 31], float) 68 65 69 ''' 70 # monthly mean parameter (1D-array) on ARCTIC SEA ICE area transformed into vector 71 spec_vec = np.zeros([M, ny * nx], float) 72 lamb_vec = np.zeros([M, ny * nx], float) 73 diff_vec = np.zeros([M, ny * nx], float) 74 ratio_vec = np.zeros([M, ny * nx], float) 75 76 # histogram distribution (intensity of occurence) of parameter in SEA ICE area (1D-array, bins = 200) 77 hist_vals_spec = np.zeros([M, bin], float) 78 hist_vals_lamb = np.zeros([M, bin], float) 79 hist_vals_diff = np.zeros([M, bin], float) 80 hist_vals_ratio = np.zeros([M, bin], float) 81 82 # histogram distribution (emissivity value) of parameter in SEA ICE area (1D-array, bins = 200) 83 corresp_emis_spec = np.zeros([M, bin], float) 84 corresp_emis_lamb = np.zeros([M, bin], float) 85 corresp_emis_diff = np.zeros([M, bin], float) 86 corresp_emis_ratio = np.zeros([M, bin], float) 87 ''' 66 88 67 months1 = np.array([0, 1, 2, 3, 8, 9, 10, 11]) # use AMSUA 23GHz to delimit ice/no_ice for JANUARY, FEBRUARY, MARCH, APRIL, SEPTEMBER, OCTOBER, NOVEMBER, DECEMBER 89 68 for imo in months1: … … 94 73 print 'open threshold file' 95 74 fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA23_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 75 spec_lim = fichier_emis.variables['spec_ice_area'][:] # sea ice pixels defined with spec emis at 23GHz 76 #lamb_lim = fichier_emis.variables['lamb_ice_area'][:] 77 fichier_emis.close() 78 ######################################################### 79 # application of AMSUA 23GHz delimitation to other data # 80 ######################################################### 81 print 'open emissivity to sub_classify file' 82 ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB) 83 fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 84 day = fichier_e.variables['days'][:] 85 emis_spec[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec'][:] 86 emis_lamb[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_lamb'][:] 87 emis_diff[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec_lamb'][:] 88 fichier_e.close() 89 # calculate emis ratio daily 90 for ijr in range (0, month_day[imo]): 91 for ilon in range (0, nx): 92 for ilat in range (0, ny): 93 emis_ratio[imo, ilat, ilon, ijr] = ((emis_lamb[imo, ilat, ilon, ijr] - emis_spec[imo, ilat, ilon, ijr]) / emis_spec[imo, ilat, ilon, ijr]) * 100. 94 # create 2D-array of emissivity spec, lamb, diff and ratio on sea ice extent only, defined by AMSUA 23GHz spec emiss threshold 95 if (isnan(spec_lim[ilat, ilon]) == True): # if pixel of sea ice extent defined with spec_emiss_23_threshold corresponds to 'no_ice', then compute NaN in pixel 96 daily_spec_ice[imo, ilat, ilon, ijr] = NaN 97 daily_lamb_ice[imo, ilat, ilon, ijr] = NaN 98 daily_diff_ice[imo, ilat, ilon, ijr] = NaN 99 daily_ratio_ice[imo, ilat, ilon, ijr] = NaN 100 else: # if pixel of sea ice extent defined with spec_emiss_23_threshold corresponds to 'ice', then compute value of emis spec, emis lamb or emis diff in pixel 101 daily_spec_ice[imo, ilat, ilon, ijr] = emis_spec[imo, ilat, ilon, ijr] 102 daily_lamb_ice[imo, ilat, ilon, ijr] = emis_lamb[imo, ilat, ilon, ijr] 103 daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr] 104 daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr] 105 ######################## 106 # stack in netcdf file # 107 ######################## 108 print 'stack in file month ' + str(month[imo]) 109 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUB' + str(frequ) + '_spec_thresholds.nc', 'w', format='NETCDF3_CLASSIC') 110 rootgrp.createDimension('longitude', nx) 111 rootgrp.createDimension('latitude', ny) 112 rootgrp.createDimension('days', month_day[imo]) 113 nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 114 nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 115 nc_days = rootgrp.createVariable('days', 'f', ('days',)) 116 nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude', 'days')) 117 nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude', 'days')) 118 nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude', 'days')) 119 nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude', 'days')) 120 nc_lon[:] = xvec 121 nc_lat[:] = yvec 122 nc_days[:] = np.arange(0, month_day[imo]) 123 nc_ice_spec[:] = daily_spec_ice[imo, :, :, 0 : month_day[imo]] 124 nc_ice_lamb[:] = daily_lamb_ice[imo, :, :, 0 : month_day[imo]] 125 nc_ice_diff[:] = daily_diff_ice[imo, :, :, 0 : month_day[imo]] 126 nc_ice_ratio[:] = daily_ratio_ice[imo, :, :, 0 : month_day[imo]] 127 rootgrp.close() 128 129 130 131 132 months2 = np.array([4, 5, 6, 7])# use AMSUA 89GHz to delimit ice/no_ice for MAY, JUNE, JULY, AUGUST 133 for imo in months2: 134 print 'month ' + month[imo] 135 ################################################################################## 136 # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study # 137 ################################################################################## 138 print 'open threshold file' 139 fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA89_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 96 140 spec_lim = fichier_emis.variables['spec_ice_area'][:] 97 141 #lamb_lim = fichier_emis.variables['lamb_ice_area'][:] … … 123 167 daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr] 124 168 daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr] 125 '''126 # ATTENTION : previous part of script has been modified ==> cannot be applied directly to this following part of script (modification of arrays and names....127 print 'compute SPEC distribution'128 ########129 # SPEC #130 ########131 cs = reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))[nonzero(isnan(reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))) == False)]132 spec_vec[imo, 0 : len(cs)] = cs133 hist_vals_spec[imo, :] = hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[0]134 for ibin in range (0, bin):135 corresp_emis_spec[imo, ibin] = mean(hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])136 print 'compute LAMB distribution'137 ########138 # LAMB #139 ########140 cl = reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))[nonzero(isnan(reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))) == False)]141 lamb_vec[imo, 0 : len(cl)] = cl142 hist_vals_lamb[imo, :] = hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[0]143 for ibin in range (0, bin):144 corresp_emis_lamb[imo, ibin] = mean(hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])145 print 'compute DIFF distribution'146 ########147 # DIFF #148 ########149 cd = reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))[nonzero(isnan(reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))) == False)]150 diff_vec[imo, 0 : len(cd)] = cd151 hist_vals_diff[imo, :] = hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[0]152 for ibin in range (0, bin):153 corresp_emis_diff[imo, ibin] = mean(hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])154 print 'compute RATIO distribution'155 #########156 # RATIO #157 #########158 cr = reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))[nonzero(isnan(reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))) == False)]159 ratio_vec[imo, 0 : len(cr)] = cr160 hist_vals_ratio[imo, :] = hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[0]161 for ibin in range (0, bin):162 corresp_emis_ratio[imo, ibin] = mean(hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])163 ######################164 # stack in .dat file #165 ######################166 print 'start stacking in .dat file'167 #data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/AMSUB'+str(frequ)+'_data_classification_parameters_ice_no-ice_with_AMSUA23-spec_2009.dat', 'a')168 for ii in range (0, bin):169 data_classif.write(('%(months)10s %(hist_vals_spec)10.5f %(corresp_emis_spec)10.5f %(hist_vals_lamb)10.5f %(corresp_emis_lamb)10.5f %(hist_vals_diff)10.5f %(corresp_emis_diff)10.5f %(hist_vals_rate)10.5f %(corresp_emis_rate)10.5f \n' % {170 'months':month[imo],171 'hist_vals_spec':hist_vals_spec[imo, ii],172 'corresp_emis_spec':corresp_emis_spec[imo, ii],173 'hist_vals_lamb':hist_vals_lamb[imo, ii],174 'corresp_emis_lamb':corresp_emis_lamb[imo, ii],175 'hist_vals_diff':hist_vals_diff[imo, ii],176 'corresp_emis_diff':corresp_emis_diff[imo, ii],177 'hist_vals_rate':hist_vals_ratio[imo, ii],178 'corresp_emis_rate':corresp_emis_ratio[imo, ii],179 }))180 '''181 169 ######################## 182 170 # stack in netcdf file # … … 206 194 207 195 208 months2 = np.array([4, 5, 6, 7])# use AMSUA 89GHz to delimit ice/no_ice for MAY, JUNE, JULY, AUGUST209 for imo in months2:210 print 'month ' + month[imo]211 ##################################################################################212 # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study #213 ##################################################################################214 print 'open threshold file'215 fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA89_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC')216 spec_lim = fichier_emis.variables['spec_ice_area'][:]217 #lamb_lim = fichier_emis.variables['lamb_ice_area'][:]218 fichier_emis.close()219 #########################################################220 # application of AMSUA 23GHz delimitation to other data #221 #########################################################222 print 'open emissivity to sub_classify file'223 ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB)224 fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')225 day = fichier_e.variables['days'][:]226 emis_spec[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec'][:]227 emis_lamb[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_lamb'][:]228 emis_diff[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec_lamb'][:]229 fichier_e.close()230 # calculate emis ratio daily231 for ijr in range (0, month_day[imo]):232 for ilon in range (0, nx):233 for ilat in range (0, ny):234 emis_ratio[imo, ilat, ilon, ijr] = ((emis_lamb[imo, ilat, ilon, ijr] - emis_spec[imo, ilat, ilon, ijr]) / emis_spec[imo, ilat, ilon, ijr]) * 100.235 if (isnan(spec_lim[ilat, ilon]) == True):236 daily_spec_ice[imo, ilat, ilon, ijr] = NaN237 daily_lamb_ice[imo, ilat, ilon, ijr] = NaN238 daily_diff_ice[imo, ilat, ilon, ijr] = NaN239 daily_ratio_ice[imo, ilat, ilon, ijr] = NaN240 else:241 daily_spec_ice[imo, ilat, ilon, ijr] = emis_spec[imo, ilat, ilon, ijr]242 daily_lamb_ice[imo, ilat, ilon, ijr] = emis_lamb[imo, ilat, ilon, ijr]243 daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr]244 daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr]245 '''246 # ATTENTION : previous part of script has been modified ==> cannot be applied directly to this following part of script (modification of arrays and names....247 print 'compute SPEC distribution'248 ########249 # SPEC #250 ########251 cs = reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))[nonzero(isnan(reshape(spec_ice[imo, :, :], size(spec_ice[imo, :, :]))) == False)]252 spec_vec[imo, 0 : len(cs)] = cs253 hist_vals_spec[imo, :] = hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[0]254 for ibin in range (0, bin):255 corresp_emis_spec[imo, ibin] = mean(hist(spec_vec[imo, 0 : len(cs)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])256 print 'compute LAMB distribution'257 ########258 # LAMB #259 ########260 cl = reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))[nonzero(isnan(reshape(lamb_ice[imo, :, :], size(lamb_ice[imo, :, :]))) == False)]261 lamb_vec[imo, 0 : len(cl)] = cl262 hist_vals_lamb[imo, :] = hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[0]263 for ibin in range (0, bin):264 corresp_emis_lamb[imo, ibin] = mean(hist(lamb_vec[imo, 0 : len(cl)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])265 print 'compute DIFF distribution'266 ########267 # DIFF #268 ########269 cd = reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))[nonzero(isnan(reshape(diff_ice[imo, :, :], size(diff_ice[imo, :, :]))) == False)]270 diff_vec[imo, 0 : len(cd)] = cd271 hist_vals_diff[imo, :] = hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[0]272 for ibin in range (0, bin):273 corresp_emis_diff[imo, ibin] = mean(hist(diff_vec[imo, 0 : len(cd)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])274 print 'compute RATIO distribution'275 #########276 # RATIO #277 #########278 cr = reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))[nonzero(isnan(reshape(ratio_ice[imo, :, :], size(ratio_ice[imo, :, :]))) == False)]279 ratio_vec[imo, 0 : len(cr)] = cr280 hist_vals_ratio[imo, :] = hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[0]281 for ibin in range (0, bin):282 corresp_emis_ratio[imo, ibin] = mean(hist(ratio_vec[imo, 0 : len(cr)], bins = bin, normed = True, histtype='step')[1][ibin : ibin + 2])283 ######################284 # stack in .dat file #285 ######################286 print 'start stacking in .dat file'287 #data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/AMSUB'+str(frequ)+'_data_classification_parameters_ice_no-ice_with_AMSUA23-spec_2009.dat', 'a')288 for ii in range (0, bin):289 data_classif.write(('%(months)10s %(hist_vals_spec)10.5f %(corresp_emis_spec)10.5f %(hist_vals_lamb)10.5f %(corresp_emis_lamb)10.5f %(hist_vals_diff)10.5f %(corresp_emis_diff)10.5f %(hist_vals_rate)10.5f %(corresp_emis_rate)10.5f \n' % {290 'months':month[imo],291 'hist_vals_spec':hist_vals_spec[imo, ii],292 'corresp_emis_spec':corresp_emis_spec[imo, ii],293 'hist_vals_lamb':hist_vals_lamb[imo, ii],294 'corresp_emis_lamb':corresp_emis_lamb[imo, ii],295 'hist_vals_diff':hist_vals_diff[imo, ii],296 'corresp_emis_diff':corresp_emis_diff[imo, ii],297 'hist_vals_rate':hist_vals_ratio[imo, ii],298 'corresp_emis_rate':corresp_emis_ratio[imo, ii],299 }))300 '''301 ########################302 # stack in netcdf file #303 ########################304 print 'stack in file month ' + str(month[imo])305 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUB' + str(frequ) + '_spec_thresholds.nc', 'w', format='NETCDF3_CLASSIC')306 rootgrp.createDimension('longitude', nx)307 rootgrp.createDimension('latitude', ny)308 rootgrp.createDimension('days', month_day[imo])309 nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))310 nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))311 nc_days = rootgrp.createVariable('days', 'f', ('days',))312 nc_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude', 'days'))313 nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude', 'days'))314 nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude', 'days'))315 nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude', 'days'))316 nc_lon[:] = xvec317 nc_lat[:] = yvec318 nc_days[:] = np.arange(0, month_day[imo])319 nc_ice_spec[:] = daily_spec_ice[imo, :, :, 0 : month_day[imo]]320 nc_ice_lamb[:] = daily_lamb_ice[imo, :, :, 0 : month_day[imo]]321 nc_ice_diff[:] = daily_diff_ice[imo, :, :, 0 : month_day[imo]]322 nc_ice_ratio[:] = daily_ratio_ice[imo, :, :, 0 : month_day[imo]]323 rootgrp.close()324 196 325 197 ''' 326 data_classif.close()327 '''328 329 330 '''331 fichier = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUA' + str(frequ) + '_spec_thresholds.nc', 'r', format='NETCDF3_CLASSIC')332 ice_spec = fichier.variables['spec_ice_area'][:]333 ice_lamb = fichier.variables['lamb_ice_area'][:]334 ice_ratio = fichier.variables['ratio_ice_area'][:]335 fichier.close()336 mean_ratio = np.zeros([ny, nx], float)337 for ilon in range (0, nx):338 for ilat in range (0, ny):339 mean_ratio[ilat, ilon] = mean(ice_ratio[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_ratio[ilat, ilon, 0 : month_day[imo]]) == False)])340 341 342 ion()343 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()344 x_coast = x_ind345 y_coast = y_ind346 z_coast = z_ind347 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, mean_ratio[:, :], -3., 5., 0.1, cm.s3pcpn_l_r, 'test')348 349 350 351 352 198 # test: 353 199 ion() -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/compare_ice_area_different_data.py
r54 r56 19 19 20 20 21 21 ######################## 22 # Time characteristics # 23 ######################## 22 24 MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) 23 25 month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER']) … … 25 27 M = len(month) 26 28 27 frequ = np.array([23, 30, 50, 89]) 29 30 frequ = np.array([23, 30, 50, 89]) # frequencies of AMSUA used for the study 28 31 29 32 … … 48 51 liste = line.split() 49 52 # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es) 50 tot_area_spec[iligne] = float(liste[9]) 51 tot_area_lamb[iligne] = float(liste[10]) 53 tot_area_spec[iligne] = float(liste[9]) # total sea ice area defined with emis_spec threshold 54 tot_area_lamb[iligne] = float(liste[10]) # total sea ice area defined with emis_lamb threshold 52 55 iligne=iligne+1 53 56 fichier.close() 54 57 vec = np.arange(0, nbtotal + 51, 50) 55 area_s = np.zeros([M], float) 56 area_l = np.zeros([M], float) 58 area_s = np.zeros([M], float) # vector of monthly mean total ice area defined with emis_spec threshold 59 area_l = np.zeros([M], float) # vector of monthly mean total ice area defined with emis_spec threshold 57 60 for imo in range (0, M): 58 61 area_s[imo] = tot_area_spec[imo + vec[imo]] 59 62 area_l[imo] = tot_area_lamb[imo + vec[imo]] 60 AS[ifr, :] = area_s[:] 61 AL[ifr, :] = area_l[:] 63 AS[ifr, :] = area_s[:] # create 2D-array of monthly mean total sea ice area defined with emis_spec threshold, for each frequency 64 AL[ifr, :] = area_l[:] # create 2D-array of monthly mean total sea ice area defined with emis_lamb threshold, for each frequency 62 65 63 66 … … 127 130 # plot time evolution of ice area # 128 131 ################################### 132 # plot various total sea ice areas for each data, monthly in 2009 129 133 ion() 130 134 figure() … … 139 143 plot(area_s_B[:], 'b', label = 'AMSUB spec_89GHz') 140 144 plot(area_l_B[:], 'b--', label = 'AMSUB lamb_89GHz') 141 plot(area_osi + 1500000, 'k', label = 'OSISAF + 1.5e7 (correction)') 145 plot(area_osi + 1500000, 'k', label = 'OSISAF + 1.5e7 (correction)') # '+1500000' corresponds to translation of OSISAF curve due to constant bias between OSISAF and AMSU total sea ice areas - quantity of translation is arbitrary 142 146 fontP = FontProperties() 143 147 fontP.set_size('small') … … 150 154 151 155 152 153 ############################### 154 # calculation of bias and std#155 ############################### 156 ''' 157 ################################################################################## 158 # calculation of bias between each AMSU total ice area and OSISAF total ice area # 159 ################################################################################## 156 160 a = np.zeros([M], float) 157 161 b = np.zeros([M], float) … … 197 201 grid() 198 202 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') 199 200 201 202 203 203 ''' 204 205 206 207 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/daily_emis_AMSUA_AMSUB_89.py
r55 r56 30 30 # grid characteristics # 31 31 ######################## 32 ### X AXIS 32 33 x0 = -3000. # min limit of grid 33 34 x1 = 2500. # max limit of grid … … 41 42 nx_b = len(xvec_b) 42 43 44 ### Y AXIS 43 45 y0 = -3000. # min limit of grid 44 46 y1 = 3000. # max limit of grid … … 85 87 86 88 87 89 ''' 90 ##################### 91 # map area of study # 92 ##################### 88 93 ion() 89 94 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() … … 92 97 z_coast = z_ind 93 98 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') 94 99 ''' 95 100 96 101 … … 171 176 172 177 178 179 180 ######################################################### 181 # Join daily evolution of each month in a yearly vector # 182 ######################################################### 173 183 mean_year_spec_a = np.append(mean_spec_a_zone[0, 0 : month_day[0]], mean_spec_a_zone[1, 0 : month_day[1]]) 174 184 mean_year_spec_a23 = np.append(mean_spec_a_zone23[0, 0 : month_day[0]], mean_spec_a_zone23[1, 0 : month_day[1]]) … … 204 214 # plot daily parameters # 205 215 ######################### 216 ### plot difference between the two zonal std (of amsua and amsub at 89GHz) 206 217 vec_months = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) 207 218 ion() … … 222 233 223 234 224 235 ### plot zonal mean of emissivity spec of AMSUA and B at 89GHz and the difference between them (subplot) 225 236 vec_months = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) 226 237 ion() … … 264 275 #title('Chukchi Sea') 265 276 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_AMSUA23_bias_std_mean_zone_Chukchi_Sea.png') 266 '''267 268 269 270 271 272 273 274 ########################################275 # monthly mean and std in whole arctic #276 ########################################277 std_spec_a89 = np.zeros([M, ny_a, nx_a], float)278 std_spec_a23 = np.zeros([M, ny_a, nx_a], float)279 std_spec_b = np.zeros([M, ny_a, nx_a], float)280 std_lamb_a89 = np.zeros([M, ny_a, nx_a], float)281 std_lamb_a23 = np.zeros([M, ny_a, nx_a], float)282 std_lamb_b = np.zeros([M, ny_a, nx_a], float)283 '''bias_spec = np.zeros([M, ny_a, nx_a], float)284 bias_anom = np.zeros([M, ny_a, nx_a], float)285 bias_mean = np.zeros([M], float)'''286 for imo in range (0, M):287 print 'month ' + month[imo]288 print 'read daily emis spec lamb and diff AMSUA and AMSUB'289 # AMSUA 89290 fichier_a89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')291 spec_a89 = fichier_a89.variables['mean_spec'][:]292 lamb_a89 = fichier_a89.variables['mean_lamb'][:]293 std_spec_a89[imo, :, :] = fichier_a89.variables['std_spec'][:]294 std_lamb_a89[imo, :, :] = fichier_a89.variables['std_lamb'][:]295 fichier_a89.close()296 # AMSUA 23297 fichier_b89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')298 spec_a23 = fichier_a23.variables['mean_spec'][:]299 lamb_a23 = fichier_a23.variables['mean_lamb'][:]300 std_spec_a23[imo, :, :] = fichier_a23.variables['std_spec'][:]301 std_lamb_a23[imo, :, :] = fichier_a23.variables['std_lamb'][:]302 fichier_a23.close()303 # AMSUB304 fichier_b = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')305 spec_b = fichier_b.variables['mean_spec'][:]306 lamb_b = fichier_b.variables['mean_lamb'][:]307 std_spec_b[imo, :, :] = fichier_b.variables['std_spec'][:]308 std_lamb_b[imo, :, :] = fichier_b.variables['std_lamb'][:]309 fichier_b.close()310 '''bias_spec[imo, :, :] = spec_a - spec_b311 bias_mean[imo] = mean(bias_spec[imo, :, :][nonzero(isnan(bias_spec[imo, :, :]) == False)])312 for ilon in range (0, nx_a):313 for ilat in range (0, ny_a):314 bias_anom[imo, ilat, ilon] = bias_spec[imo, ilat, ilon] - bias_mean[imo]'''315 316 # calculate difference of std between spec and lamb for amsua and amsub317 c = np.zeros([M], float)318 f = np.zeros([M], float)319 for imo in range (0, M):320 a = mean(std_spec_a89[imo, :, :][nonzero(isnan(std_spec_a89[imo, :, :])==False)])321 b = mean(std_lamb_a89[imo, :, :][nonzero(isnan(std_lamb_a89[imo, :, :])==False)])322 c[imo] = a - b323 d = mean(std_spec_b[imo, :, :][nonzero(isnan(std_spec_b[imo, :, :])==False)])324 e = mean(std_lamb_b[imo, :, :][nonzero(isnan(std_lamb_b[imo, :, :])==False)])325 f[imo] = d - e326 327 # calculate difference of std between spec amsua and spec amsub328 k = np.zeros([M], float)329 for imo in range (0, M):330 g = mean(std_spec_a89[imo, :, :][nonzero(isnan(std_spec_a89[imo, :, :])==False)])331 h = mean(std_spec_b[imo, :, :][nonzero(isnan(std_spec_b[imo, :, :])==False)])332 k[imo] = abs(g - h)333 334 335 336 ##################################337 # map bias spec_AMSUA-spec_AMSUB #338 ##################################339 #ion()340 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()341 x_coast = x_ind342 y_coast = y_ind343 z_coast = z_ind344 for imo in range (0, M):345 print 'map for month ' + month[imo]346 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec_a, yvec_a, std_lamb_a89[imo, :, :], 0., 0.12, 0.001, cm.s3pcpn_l_r, 'std emis lamb AMSUA89')347 title(month[imo] + ' 2009')348 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/monthly_std/std_emis_lamb_AMSUA89_' + month[imo] + '2009.png')349 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/daily_evolution_anomaly_grad_ratio.py
r52 r56 43 43 44 44 45 ''' 45 46 46 ############################################ 47 47 # time evolution (monthly) in a given zone # … … 66 66 len(yvec[yyi:yyf+1]) 67 67 68 mean_grad_ratio_zone = np.zeros([M, 31], float) 69 std_grad_ratio_zone = np.zeros([M, 31], float) 70 71 mean_spec_anom_23_zone = np.zeros([M, 31], float) 72 std_spec_anom_23_zone = np.zeros([M, 31], float) 73 74 mean_spec_anom_89_zone = np.zeros([M, 31], float) 75 std_spec_anom_89_zone = np.zeros([M, 31], float) 76 77 mean_ratio_anom_89_zone = np.zeros([M, 31], float) 78 std_ratio_anom_89_zone = np.zeros([M, 31], float) 79 80 S = np.zeros([M, 31], float) 68 69 mean_grad_ratio_zone = np.zeros([M, 31], float) # 2D-array of zonal mean gradient ratio for each day in each month 70 std_grad_ratio_zone = np.zeros([M, 31], float) # 2D-array of zonal std of gradient ratio for each day in each month 71 72 mean_spec_anom_23_zone = np.zeros([M, 31], float) # 2D-array of zonal mean emis_spec spatial anomaly for each day in each month (at 23GHz) 73 std_spec_anom_23_zone = np.zeros([M, 31], float) # 2D-array of zonal std of emis_spec spatial anomaly for each day in each month (at 23GHz) 74 75 mean_spec_anom_89_zone = np.zeros([M, 31], float ) # 2D-array of zonal mean emis_spec spatial anomaly for each day in each month (at 89GHz) 76 std_spec_anom_89_zone = np.zeros([M, 31], float) # 2D-array of zonal std of emis_spec spatial anomaly for each day in each month (at 89GHz) 77 78 mean_ratio_anom_89_zone = np.zeros([M, 31], float) # 2D-array of zonal mean emis_ratio spatial anomaly for each day in each month (at 89GHz) 79 std_ratio_anom_89_zone = np.zeros([M, 31], float) # 2D-array of zonal std of emis_ratio spatial anomaly for each day in each month (at 89GHz) 80 81 S = np.zeros([M, 31], float) # number of data pixels in selected zone, per day and per month 81 82 82 83 for imo in range (0, M): … … 165 166 166 167 167 ############################################################################ 168 # calculate standard deviation of emissivity parameters over the year 2009 #169 ############################################################################ 168 ###################################################################################### 169 # calculate standard deviation of emissivity parameters over the year 2009 (364 days)# 170 ###################################################################################### 170 171 print 'calculate standard deviation of emissivity parameters over the year 2009' 171 172 time_std1 = sqrt((1./364.)*sum((mean_year_grad_ratio[nonzero(isnan(mean_year_grad_ratio) == False)] - mean(mean_year_grad_ratio[nonzero(isnan(mean_year_grad_ratio) == False)]))**2)) … … 244 245 title('area of study - Chukchi Sea') 245 246 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') 246 ''' 247 248 249 ################################ 250 # mapparameters in all Arctic #251 ################################ 247 248 249 250 ############################################ 251 # read emissivity parameters in all Arctic # 252 ############################################ 252 253 cumul_params = np.zeros([M, ny, nx], float) 253 254 grad_ratio = np.zeros([M, ny, nx], float) … … 275 276 for ilon in range (0, nx): 276 277 for ilat in range (0, ny): 278 # calculate monthly mean of emissivity parameters 277 279 grad_ratio[imo, ilat, ilon] = mean(gr[:, ilat, ilon][nonzero(isnan(gr[:, ilat, ilon]) == False)]) 278 280 spec_anom_23[imo, ilat, ilon] = mean(sa_23[:, ilat, ilon][nonzero(isnan(sa_23[:, ilat, ilon]) == False)]) 279 281 spec_anom_89[imo, ilat, ilon] = mean(sa_89[:, ilat, ilon][nonzero(isnan(sa_89[:, ilat, ilon]) == False)]) 280 282 ratio_anom_89[imo, ilat, ilon] = mean(ra_89[:, ilat, ilon][nonzero(isnan(ra_89[:, ilat, ilon]) == False)]) 283 # calculate monthly cumulation index = (sum(abs(all emissivity parameters))/maximum of this sum)*100 (to have a percentage) 281 284 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]) 282 285 CP[imo, :, :] = (cumul_params[imo, :, :]/max(cumul_params[imo, :, :][nonzero(isnan(cumul_params[imo, :, :]) == False)])) * 100. 283 286 284 287 288 ####################### 289 # map cuulation index # 290 ####################### 285 291 #ion() 286 292 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/daily_gradient_ratio.py
r51 r56 49 49 # 23GHz 50 50 print 'read daily emis spec 23GHz' 51 fichier_23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') 51 fichier_23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') 52 52 emis_spec_23 = fichier_23.variables['e_spec'][:] 53 53 fichier_23.close() -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/distance_lon_lat.py
r40 r56 15 15 from matplotlib import colors 16 16 17 17 ############################################# 18 # This function gives the length of the arc # 19 # which links two points (1 and 2) defined # 20 # by their polar lon/lat coordinates # 21 ############################################# 18 22 19 23 def distance_on_unit_sphere(lat1, lon1, lat2, lon2): -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/evolution_espec_elamb_fov.py
r40 r56 17 17 18 18 ##################### 19 # read fi chier.DAT #19 # read file .DAT # 20 20 ##################### 21 21 fichier=open('/mma/hermozol/Documents/Data_tests/GLACE/GLACE_AMSUB89_EMIS_SEPTEMBER2009.DAT','r') … … 95 95 ####################################### 96 96 r = pi / 180. 97 opac = - cos(zen_angle * r) * np.log(transm) 98 E3 = scipy.special.expn(3, opac) 99 theta_eff = np.zeros([L], float) 100 tdn_lamb = np.zeros([L], float) 101 tb_spec = np.zeros([L], float) 102 tb_lamb = np.zeros([L], float) 103 e_spec = np.zeros([L], float) 104 e_lamb = np.zeros([L], float) 105 e_spec_lamb = np.zeros([L], float) 97 opac = - cos(zen_angle * r) * np.log(transm) # opacity 98 E3 = scipy.special.expn(3, opac) # exponentielle d ordre 3 99 theta_eff = np.zeros([L], float) # effective incidence angle from Matzler method 100 tdn_lamb = np.zeros([L], float) # downwelling radiation 101 tb_spec = np.zeros([L], float) # specular brightness temperature 102 tb_lamb = np.zeros([L], float) # lambertian brightness temperature 103 e_spec = np.zeros([L], float) # specular emissivity 104 e_lamb = np.zeros([L], float) # lambertian emissivity 105 e_spec_lamb = np.zeros([L], float) # difference between lambertian and specular emissivity 106 106 for ii in range (0, L): 107 107 theta_eff[ii] = math.acos(- opac[ii] / (np.log(2. * E3[ii]))) * (1. / r) -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/ffgrid2.py
r40 r56 1 # griddata.py - 2010-07-11 ccampo2 1 import numpy as np 3 2 from numpy import * -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/ffgrid_evolution_espec_elamb_fov.py
r40 r56 15 15 16 16 17 18 19 ####################################################################### 20 # Open .dat file to read emis spec and emis lamb in September and fov # 21 ####################################################################### 17 22 fichier=open('/mma/hermozol/Documents/Data_tests/monthly_GLACE/lamb_spec_param_all_angles_SEPTEMBER2009.dat','r') 18 23 numlines = 0 … … 21 26 22 27 fichier.close 23 24 28 fichier = open('/mma/hermozol/Documents/Data_tests/monthly_GLACE/lamb_spec_param_all_angles_SEPTEMBER2009.dat','r') 25 29 nbtotal = numlines-1 … … 40 44 iligne=iligne+1 41 45 42 43 44 46 fichier.close() 45 47 46 47 48 # compute fov as x axis in scatter plot 48 49 x0 = fov.min() 49 50 x1 = fov.max() 51 # compute diff emis_lamb-emis_spec as y axis in scatter plot 50 52 y0 = e_spec_lamb.min() 51 53 y1 = e_spec_lamb.max() 54 # to each observation of coord x and y, add z = value of emis_spec 52 55 z0 = e_spec.min() 53 56 z1 = e_spec.max() 57 # step of x and y axis 54 58 dx = 1. 55 59 dy = 0.001 60 # use ffgrid to create a new grid of fov/diff(lamb-spec) 56 61 zgrid, xvec, yvec = ffgrid2.ffgrid(fov, e_spec_lamb, e_spec, dx, dy, x0, x1, y0, y1, z0, z1) 57 62 xii, yii = np.meshgrid(xvec, yvec) 63 64 ############# 65 # draw grid # 66 ############# 58 67 plt.ion() 59 68 figure() -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/ice_delimit_AMSUB_data.py
r54 r56 41 41 42 42 43 frequ = 89 43 frequ = 89 # define which frequency to use for the study 44 44 45 45 ##################### 46 46 # read NETCDF files # 47 47 ##################### 48 emis_spec_month = np.zeros([M, ny, nx], float) 49 emis_lamb_month = np.zeros([M, ny, nx], float) 50 ratio_emis = np.zeros([M, ny, nx], float) 48 49 emis_spec_month = np.zeros([M, ny, nx], float) # monthly mean of emis_spec 50 emis_lamb_month = np.zeros([M, ny, nx], float) # monthly mean of emis_lamb 51 ratio_emis = np.zeros([M, ny, nx], float) # monthly mean of emis_ratio 52 53 # conversion of previous data into 1D arrays 51 54 ratio_emis_vec = np.zeros([M, ny * nx], float) 52 55 spec_emis_vec = np.zeros([M, ny * nx], float) 53 56 lamb_emis_vec = np.zeros([M, ny * nx], float) 57 54 58 print 'read data from files' 55 59 for imo in range (0, M): … … 93 97 ''' 94 98 95 96 c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) 99 ############################################################## 100 # definition of parameters of the histograms of distribution # 101 ############################################################## 102 c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) # definition of one color per month 97 103 #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 98 104 #limit_coef_lamb = np.array([0.6, 0.6, 0.8, 0.8, 0.77]) # i1 = AMSUA23GHz / i2 = AMSUA30GHz / i3 = AMSUA50GHz / i4 = AMSUA89GHz / i5 = AMSUB89GHz 99 #idata = 0 100 fontP = FontProperties() 105 #idata = 0 # defines which frequency we use (0=A23, 1=A30, 2=A50, 3=A89, 4=B89) 106 fontP = FontProperties() # concerns legend font in plots 101 107 fontP.set_size('small') 108 109 102 110 print 'distribution of emissivity values' 103 111 ################################### … … 108 116 corresp_emis_spec = np.zeros([M, 50], float) 109 117 for imo in range (0, M): 110 hist_vals_spec[imo, :] = hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 118 hist_vals_spec[imo, :] = hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] # calculate histogram distribution of emissivity 111 119 for ibin in range (0, 50): 112 corresp_emis_spec[imo, ibin] = mean(hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 120 corresp_emis_spec[imo, ibin] = mean(hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) # calculate corresponding emissivity value to histogram distribution (in frequency of occurence) 113 121 114 122 ## plot first six months of spec emissivity histograms ## -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/map_data_lamb_spec_near_nadir_netcdf.py
r40 r56 8 8 from mpl_toolkits.basemap import shiftgrid, cm 9 9 from netCDF4 import Dataset 10 import arctic_map # function to regrid coast limits11 import cartesian_grid_test # function to convert grid from polar to cartesian12 10 import scipy.special 13 11 import ffgrid2 … … 22 20 23 21 24 dx=1. 25 dy=0.5 22 ###################################### 23 # polar lon/lat grid characteristics # 24 ###################################### 25 dx=1. # longitude every 1 deg 26 dy=0.5 # latitude every 0.5 deg 26 27 x0, x1 = -180., 180. 27 y0, y1 = 65., 90. 28 y0, y1 = 65., 90. # only high latitudes 28 29 xvec = np.arange(x0, x1 + dx, dx) 29 30 yvec = np.arange(y0, y1 + dy, dy) … … 32 33 33 34 34 e_spec_lamb_month = np.zeros([ny, nx, M], float) 35 #################################################################################### 36 # read emis_spec, emis_lamb and emis_diff data in monthly NETCDF files for mapping # 37 #################################################################################### 38 39 e_spec_lamb_month = np.zeros([ny, nx, M], float) # monthly mean of diff(emis_lamb-emis_spec) 35 40 e_spec_month = np.zeros([ny, nx, M], float) 41 36 42 for imo in range (0, M): 37 43 fichier = Dataset('/mma/hermozol/Documents/Data_tests/monthly_GLACE/gridded_data/grid_monthly_data_lamb_spec_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') … … 45 51 e_lamb = fichier.variables['e_lamb'][:] 46 52 e_spec_lamb = fichier.variables['e_spec_lamb'][:] 47 e_sl_00 = fichier.variables['e_mixed_s00'][:] 48 e_sl_25 = fichier.variables['e_mixed_s25'][:] 49 e_sl_50 = fichier.variables['e_mixed_s50'][:] 50 e_sl_75 = fichier.variables['e_mixed_s75'][:] 51 e_sl_100 = fichier.variables['e_mixed_s100'][:] 53 e_sl_00 = fichier.variables['e_mixed_s00'][:] # emissivity specular and lambertian with specular coefficient of 0% (only lambertian emissivity) 54 e_sl_25 = fichier.variables['e_mixed_s25'][:] # emissivity specular and lambertian with specular coefficient of 25% 55 e_sl_50 = fichier.variables['e_mixed_s50'][:] # emissivity specular and lambertian with specular coefficient of 50% 56 e_sl_75 = fichier.variables['e_mixed_s75'][:] # emissivity specular and lambertian with specular coefficient of 75% 57 e_sl_100 = fichier.variables['e_mixed_s100'][:] # emissivity specular and lambertian with specular coefficient of 100% (only specular emissivity) 52 58 fichier.close() 53 e_spec_month[:, :, imo] = e_spec[:, :] 54 e_spec_lamb_month[:, :, imo] = e_spec_lamb[:, :] 59 e_spec_month[:, :, imo] = e_spec[:, :] # monthly mean of emis_spec 60 e_spec_lamb_month[:, :, imo] = e_spec_lamb[:, :] # monthly mean of diff(emis_lamb-emis_spec) 55 61 56 62 57 63 ####### 64 # map # 65 ####### 66 print 'start mapping' 58 67 cm = plt.cm.get_cmap('RdBu') 59 68 #cmap = cm.s3pcpn_l_r 60 69 plt.ion() 61 70 for imo in range (0, M): 71 print 'month ' + month[imo] 62 72 map_ffgrid.draw_map_l(lon, lat, 0.6, 1.02, 0.01, e_spec_month[:,:,imo], cm.s3pcpn_l_r, 'emissivity SPEC') 63 73 title(month[imo] + ' 2009 - AMSUB') -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/map_emis_spec_diff_all_arctic.py
r55 r56 41 41 42 42 43 frequ = np.array(['23', '30', '50', '89']) 43 frequ = np.array(['23', '30', '50', '89']) # choose studied frequency 44 44 45 46 ####################################################################### 47 # Reads monthly mean and std of emis spec, lamb and diff NETCDF files # 48 ####################################################################### 45 49 mean_spec = np.zeros([4, M, ny, nx], float) 46 50 mean_lamb = np.zeros([4, M, ny, nx], float) 47 51 mean_lamb_spec = np.zeros([4, M, ny, nx], float) 52 48 53 std_spec = np.zeros([4, M, ny, nx], float) 49 54 std_lamb = np.zeros([4, M, ny, nx], float) 55 50 56 for ifr in range (0, 4): 51 57 print 'frequency ' + frequ[ifr] … … 61 67 62 68 63 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, mean_spec[0, 10, :, :], 0.45, 1.02, 0.001, colormap, 'emissivity LAMB')64 69 65 70 66 71 ################################################# 72 # map of spec, lamb or diff of emis std or mean # 73 ################################################# 67 74 ion() 68 75 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/map_ffgrid.py
r40 r56 10 10 11 11 12 12 ###################################################### 13 # Function that maps 2D-array data in a lon/lat grid # 14 # computed in ffgrid function # 15 ###################################################### 13 16 14 17 15 18 def draw_map_l (xcoord, ycoord, c0, c1, dc, zvar, colormap, text): 16 19 20 ############################################################## 21 # xcoord = vector of longitudes # 22 # ycoord = vector of latitudes # 23 # c0, c1, dc = boundaries and step of scale of the colorbar # 24 # zvar = 2D-array of data to map # 25 # colormap = colors to use for mapping od data # 26 # text = title of the colorbar # 27 ############################################################## 28 17 29 figure() 18 30 m = Basemap(projection = 'nplaea', boundinglat = 65., lon_0 = 0., resolution = 'l', fix_aspect = True) -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/map_statistic_parameters.py
r51 r56 63 63 print 'stack in file month ' + str(month[imo]) 64 64 print 'open anomaly file' 65 fichier_anom = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and- 30_param_anomaly_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '.nc', 'r', format='NETCDF3_CLASSIC')65 fichier_anom = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_param_anomaly_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '.nc', 'r', format='NETCDF3_CLASSIC') 66 66 anom_spec[ifr, imo, :, :] = fichier_anom.variables['spec_anomaly'][:] 67 67 anom_lamb[ifr, imo, :, :] = fichier_anom.variables['lamb_anomay'][:] … … 70 70 fichier_anom.close() 71 71 print 'open ice file' 72 fichier_ice = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and- 30_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC')72 fichier_ice = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUA' + str(frequ[ifr]) + '_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 73 73 ice_spec[ifr, imo, :, :] = fichier_ice.variables['spec_ice_area'][:] 74 74 ice_lamb[ifr, imo, :, :] = fichier_ice.variables['lamb_ice_area'][:] … … 99 99 100 100 101 ################################################ 102 # calculate gradient ratio fortwo frequencies #103 ################################################ 101 #################################################### 102 # calculate gradient ratio between two frequencies # 103 #################################################### 104 104 print 'calculate gradient ratio' 105 105 grad_ratio = np.zeros([M, ny, nx], float) … … 111 111 for ilon in range (0, nx): 112 112 grad_ratio[imo, ilat, ilon] = (ice_spec[0, imo, ilat, ilon] - ice_spec[3, imo, ilat, ilon]) / (ice_spec[0, imo, ilat, ilon] + ice_spec[3, imo, ilat, ilon]) 113 ''' 114 # compute histogram distribution of gradient ratio to find characteristic threshold 115 grad_ratio_vec = reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))[nonzero(isnan(reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))) == False)] 116 hist_ratio[imo, :] = hist(grad_ratio_vec, bins = 50, normed = True, histtype='step')[0] 117 for ibin in range (0, 50): 118 corresp_ratio[imo, ibin] = mean(hist(grad_ratio_vec, bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 119 120 # plot histogram 121 c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) 122 figure() 123 for imo in range (0, 6): 124 plot(corresp_ratio[imo], hist_ratio[imo], c = str(c[imo]), label = str(month[imo])) 125 126 grid() 127 xlim(corresp_ratio.min(), corresp_ratio.max()) 128 xlabel('emissivity ratio') 129 ylabel('frequency of occurence') 130 fontP = FontProperties() 131 fontP.set_size('small') 132 legend(loc = 2, prop = fontP) 133 #plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_spec_AMSUA' + str(frequ) + '_JANUARY-JUNE_2009.png') 134 ## plot six following months of spec emissivity histograms ## 135 figure() 136 for imo in range (6, M): 137 plot(corresp_ratio[imo], hist_ratio[imo], c = str(c[imo - 6]), label = str(month[imo])) 138 139 grid() 140 xlim(corresp_ratio.min(), corresp_ratio.max()) 141 xlabel('emissivity ratio') 142 ylabel('frequency of occurence') 143 legend(loc = 1, prop = fontP) 144 #plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/emiss_spec_AMSUA' + str(frequ) + '_JULY-DECEMBER_2009.png') 145 ''' 113 146 114 147 115 ################################ … … 184 152 185 153 186 187 '''188 figure()189 scatter(x_coast, y_coast, c = z_coast, s = volume)190 #cmap = cm.s3pcpn_l_r191 levels = np.arange(-1, 11, 1)192 cs = contour(xvec, yvec, cumul[imo, :, :])193 cs = contour(xvec, yvec, new_anom_spec0[imo, :, :], levels)194 cs = contour(xvec, yvec, new_anom_ratio[imo, :, :], levels)195 cs = contour(xvec, yvec, new_grad_ratio[imo, :, :], levels)196 cbar = colorbar(cs, ticks = [-1, 0, 1, 2, 3, 4, 5])197 cbar.ax.set_yticklabels(['emis_spec_89', 'emis_spec_23', 'emis_ratio_89', 'grad_ratio_23-89', '-'])198 #cbar.set_label(text)199 xlim(-3500, 2700.)200 ylim(-4000, 2800.)201 '''202 203 204 205 206 207 208 209 '''210 ############################################211 # time evolution (monthly) in a given zone #212 ############################################213 #zone 1 : yi = 400. // yf = 720. // xi = -1080. // xf = -920. (Canadian Archipelago)214 #zone 2 : yi = 920. // yf = 1280. // xi = -720. // xf = -560. (Chukchi Sea, Bering Detroit)215 # select borders of zone216 yi = 920.217 yf = 1280.218 xi = -720.219 xf = -560.220 #find corresponding index in xvec and yvec221 xxi = np.where(xvec == xi)[0][0]222 xxf = np.where(xvec == xf)[0][0]223 yyi = np.where(yvec == yi)[0][0]224 yyf = np.where(yvec == yf)[0][0]225 226 # define vectors227 anom_spec0_zone = np.zeros([M], float)228 anom_spec3_zone = np.zeros([M], float)229 anom_ratio_zone = np.zeros([M], float)230 grad_ratio_zone = np.zeros([M], float)231 std_as0 = np.zeros([M], float)232 std_as3 = np.zeros([M], float)233 std_ar = np.zeros([M], float)234 std_gr = np.zeros([M], float)235 for imo in range (0, M):236 # anom spec emis from 23GHz237 anom_spec0_vec = reshape(anom_spec[0, imo, yyi : yyf + 1, xxi : xxf + 1], size(anom_spec[0, imo, yyi : yyf + 1, xxi : xxf + 1]))238 if (len(anom_spec0_vec[nonzero(isnan(anom_spec0_vec) == True)]) != 0):239 anom_spec0_zone[imo] = NaN240 std_as0[imo] = NaN241 anom_spec3_zone[imo] = NaN242 std_as3[imo] = NaN243 anom_ratio_zone[imo] = NaN244 std_ar[imo] = NaN245 grad_ratio_zone[imo] = NaN246 std_gr[imo] = NaN247 continue248 else:249 anom_spec0_zone[imo] = anom_spec0_vec.mean()250 std_as0[imo] = sqrt((1. / (size(anom_spec[0, imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((anom_spec0_vec[:] - anom_spec0_zone[imo])**2))251 # anom spec emis from 89GHz252 anom_spec3_vec = reshape(anom_spec[3, imo, yyi : yyf + 1, xxi : xxf + 1], size(anom_spec[3, imo, yyi : yyf + 1, xxi : xxf + 1]))253 anom_spec3_zone[imo] = anom_spec3_vec.mean()254 std_as3[imo] = sqrt((1. / (size(anom_spec[3, imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((anom_spec3_vec[:] - anom_spec3_zone[imo])**2))255 # anom emis ratio from 89GHz256 anom_ratio_vec = reshape(anom_ratio[3, imo, yyi : yyf + 1, xxi : xxf + 1], size(anom_ratio[3, imo, yyi : yyf + 1, xxi : xxf + 1]))257 anom_ratio_zone[imo] = anom_ratio_vec.mean()258 std_ar[imo] = sqrt((1. / (size(anom_ratio[3, imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((anom_ratio_vec[:] - anom_ratio_zone[imo])**2))259 # gradient ratio from 23 and 89GHz260 grad_ratio_vec = reshape(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1], size(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1]))261 grad_ratio_zone[imo] = grad_ratio_vec.mean()262 std_gr[imo] = sqrt((1. / (size(grad_ratio[imo, yyi : yyf + 1, xxi : xxf + 1]) - 1.)) * sum((grad_ratio_vec[:] - grad_ratio_zone[imo])**2))263 264 265 figure()266 fontP = FontProperties()267 fontP.set_size('small')268 subplot(2,1,1)269 errorbar(arange(0, M, 1), anom_spec0_zone, yerr = std_as0, label = 'emis spec anomaly 23GHz')270 errorbar(arange(0, M, 1), anom_spec3_zone, yerr = std_as3, label = 'emis spec anomaly 89GHz')271 errorbar(arange(0, M, 1), grad_ratio_zone, yerr = std_gr, label = 'gradient ratio 23-89GHz')272 plot(arange(0, M, 1), np.zeros([M], float), '--k')273 legend(loc = 3, prop = fontP)274 grid()275 xticks(arange(0, M, 1), month, rotation = 15)276 xlim(-1, M)277 subplot(2,1,2)278 errorbar(arange(0, M, 1), anom_ratio_zone, yerr = std_ar, label = 'emis ratio anomaly 89GHz')279 plot(arange(0, M, 1), np.zeros([M], float), '--k')280 legend(loc = 2, prop = fontP)281 grid()282 xticks(arange(0, M, 1), month, rotation = 15)283 xlim(-1, M)284 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/study_by_zones/monthly_evolution_emis_params_zone_North_Beaufort_Sea.png')285 ### map the selected zone ###286 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec[xxi : xxf], yvec[yyi : yyf], grad_ratio[imo, yyi : yyf, xxi : xxf], grad_ratio[imo, yyi : yyf, xxi : xxf][nonzero(isnan(grad_ratio[imo, yyi : yyf, xxi : xxf]) == False)].min(), grad_ratio[imo, yyi : yyf, xxi : xxf][nonzero(isnan(grad_ratio[imo, yyi : yyf, xxi : xxf]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'selected zone')287 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/study_by_zones/map_zone_North_Beaufort_Sea.png')288 '''289 290 291 292 '''293 154 ############################################################################ 294 155 # calculate correlation between anom ratio and gradient ratio (spec emiss) # … … 304 165 for jj in range (0, 4): 305 166 corr_mat[imo, ii, jj] = corrcoef(params[ii], params[jj])[0][1] 167 306 168 307 169 figure() … … 322 184 grid() 323 185 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/correlation_params_2009.png') 324 ''' 186 325 187 326 188 ''' -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/param_anomaly.py
r55 r56 57 57 ice_diff_anom = np.zeros([4, M, ny, nx], float) 58 58 ice_ratio_anom = np.zeros([4, M, ny, nx], float) 59 59 60 for ifr in range (3, 4): 60 61 print 'frequency ' + frequ[ifr] … … 69 70 for ilat in range (0, ny): 70 71 for ilon in range (0, nx): 72 # calculate monthly mean of emis spec, lamb, diff and ratio 71 73 mean_ice_spec[ifr, imo, ilat, ilon] = mean(ice_spec[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_spec[ilat, ilon, 0 : month_day[imo]]) == False)]) 72 74 mean_ice_lamb[ifr, imo, ilat, ilon] = mean(ice_lamb[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_lamb[ilat, ilon, 0 : month_day[imo]]) == False)]) 73 75 mean_ice_diff[ifr, imo, ilat, ilon] = mean(ice_diff[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_diff[ilat, ilon, 0 : month_day[imo]]) == False)]) 74 76 mean_ice_ratio[ifr, imo, ilat, ilon] = mean(ice_ratio[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(ice_ratio[ilat, ilon, 0 : month_day[imo]]) == False)]) 77 # transforms 2D-arrays of monthly mean of emissivity into vectors, in order to calculate the average emissivity over total sea ice extent, for each month 75 78 a = reshape(mean_ice_spec[ifr, imo, :, :], size(mean_ice_spec[ifr, imo, :, :])) 76 79 b = reshape(mean_ice_lamb[ifr, imo, :, :], size(mean_ice_lamb[ifr, imo, :, :])) … … 80 83 for ilat in range (0, ny): 81 84 for ilon in range (0, nx): 85 # calculate monthly spatial anomaly in each pixel 82 86 ice_spec_anom[ifr, imo, ilat, ilon] = mean_ice_spec[ifr, imo, ilat, ilon] - mean(a[nonzero(isnan(a) == False)]) 83 87 ice_lamb_anom[ifr, imo, ilat, ilon] = mean_ice_lamb[ifr, imo, ilat, ilon] - mean(b[nonzero(isnan(b) == False)]) … … 107 111 108 112 109 '''110 111 # test:112 ion()113 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()114 x_coast = x_ind115 y_coast = y_ind116 z_coast = z_ind117 for ifr in range (0, 2):118 for imo in range (0, M):119 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ice_ratio_anom[ifr, imo, :, :], ice_ratio_anom[nonzero(isnan(ice_ratio_anom) == False)].min(), ice_ratio_anom[nonzero(isnan(ice_ratio_anom) == False)].max(), 0.0001, cm.s3pcpn_l_r, 'Emissivity diff anomaly')120 title('AMSUA ' + str(frequ) + ' - ' + str(month[imo]) + ' 2009')121 plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/emiss_diff_anomaly_map_AMSUA'+str(frequ[ifr])+'_with_AMSUA23-and-30-spec_'+str(month[imo])+'_2009.png')122 #cm.s3pcpn_l_r123 124 '''125 113 126 114 127 128 129 ''' 130 frequ = 89 131 spec_ice = np.zeros([M, ny, nx], float) 132 lamb_ice = np.zeros([M, ny, nx], float) 133 diff_ice = np.zeros([M, ny, nx], float) 134 spec_month = np.zeros([M, ny, nx], float) 135 lamb_month = np.zeros([M, ny, nx], float) 136 diff_month = np.zeros([M, ny, nx], float) 137 std_spec = np.zeros([M, ny, nx], float) 138 std_lamb = np.zeros([M, ny, nx], float) 139 std_diff = np.zeros([M, ny, nx], float) 140 for imo in range (0, 1): 141 print 'month ' + month[imo] 142 ################################################################################## 143 # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study # 144 ################################################################################## 145 print 'open threshold file' 146 fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA23_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC') 147 spec_lim = fichier_emis.variables['spec_ice_area'][:] 148 fichier_emis.close() 149 fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 150 day = fichier_e.variables['days'][:] 151 emis_spec = fichier_e.variables['e_spec'][:] 152 emis_lamb = fichier_e.variables['e_lamb'][:] 153 emis_diff = fichier_e.variables['e_spec_lamb'][:] 154 fichier_e.close() 155 for ilon in range (0, nx): 156 for ilat in range (0, ny): 157 if (isnan(spec_lim[ilat, ilon]) == True): 158 spec_ice[imo, ilat, ilon] = NaN 159 lamb_ice[imo, ilat, ilon] = NaN 160 diff_ice[imo, ilat, ilon] = NaN 161 else: 162 spec_ice[imo, ilat, ilon] = spec_month[imo, ilat, ilon] 163 lamb_ice[imo, ilat, ilon] = lamb_month[imo, ilat, ilon] 164 diff_ice[imo, ilat, ilon] = diff_month[imo, ilat, ilon] 165 spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, :][nonzero(isnan(emis_spec[ilat, ilon, :]) == False)]) 166 lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, :][nonzero(isnan(emis_lamb[ilat, ilon, :]) == False)]) 167 diff_month[imo, ilat, ilon] = mean(emis_diff[ilat, ilon, :][nonzero(isnan(emis_diff[ilat, ilon, :]) == False)]) 168 Ss = sum((emis_spec[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_spec[ilat, ilon, 0 : month_day[imo]]) == False)] - spec_month[imo, ilat, ilon])**2) 169 std_spec[imo, ilat, ilon] = sqrt((1. / month_day[imo] - 1.) * Ss) 170 Sl = sum((emis_lamb[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_lamb[ilat, ilon, 0 : month_day[imo]]) == False)] - lamb_month[imo, ilat, ilon])**2) 171 std_lamb[imo, ilat, ilon] = sqrt((1. / month_day[imo] - 1.) * Sl) 172 173 174 175 176 # test: 115 ####################### 116 # map of emis anomaly # 117 ####################### 177 118 ion() 178 119 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() … … 186 127 #cm.s3pcpn_l_r 187 128 188 imo = 1 189 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec[imo, :, :], std_spec[imo,:,:][nonzero(isnan(std_spec[imo,:,:]) == False)].min(), std_spec[imo, :, :][nonzero(isnan(std_spec[imo , :, :]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'std spec') 190 ''' 129 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_ffgrid_osisaf_ice_types.py
r40 r56 38 38 39 39 40 osi_type = np.zeros([M, 31, ny, nx], float) 41 spec = np.zeros([M, 31, ny, nx], float) 42 lamb = np.zeros([M, 31, ny, nx], float) 40 osi_type = np.zeros([M, 31, ny, nx], float) # daily ice type in polar lon/lat grid 41 spec = np.zeros([M, 31, ny, nx], float) # daily emissivity spec in polar lon/lat grid 42 lamb = np.zeros([M, 31, ny, nx], float) # daily emissivity lamb in polar lon/lat grid 43 43 44 for imo in range (3, 4): 44 45 print 'month ' + month[imo] … … 76 77 iligne=iligne+1 77 78 fichier_amsub.close() 78 for ijr in range (0, month_day[imo]): 79 for ijr in range (0, month_day[imo]): # read daily data 80 # stop computing for this month if day > maximum day in month 79 81 if ((ijr + 1) > month_day[imo]): continue 80 82 else: … … 84 86 continue 85 87 else: 88 # re-grid in polar lon/lat daily emis spec and emis_lamb 86 89 print 'date ', jjr[bbjr][0] 87 90 lat_jr = lat[bbjr] … … 129 132 130 133 131 134 ''' 132 135 ######### 133 136 # tests # … … 140 143 cmap2 = colors.ListedColormap(['1.', '0.25', '0.50', '0.75']) 141 144 map_ffgrid.draw_map_l (xvec, yvec, 0, 5, 1, osi_type[11, 19, :, :] , cmap2, 'OSISAF ice type') 145 ''' 142 146 143 147 … … 148 152 149 153 150 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_spec_glace_amb89.py
r40 r56 17 17 month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) 18 18 M = len(month) 19 20 19 21 20 22 … … 95 97 96 98 97 bbjr = nonzero(jjr == 1.)98 bbzen = nonzero(zen[bbjr] > 20.)99 z0 = e[nonzero(isnan(e) == False)].min()100 z1 = e[nonzero(isnan(e) == False)].max()101 zgrid, xvec, yvec = ffgrid2.ffgrid(llon, llat, e, 1., 1., -180., 180., 65., 90., z0, z1)102 emis_amsub_89 = zgrid103 plt.ion()104 105 map_ffgrid.draw_map_l(, lat, emis.min(), emis.max(), 0.01, emis, 'emissivity')106 107 108 109 110 111 112 99 113 100 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_spec_lamb_nadir.py
r53 r56 31 31 x0 = -3000. # min limit of grid 32 32 x1 = 2500. # max limit of grid 33 dx = 100.33 dx = 40. # x resolution of grid 34 34 xvec = np.arange(x0, x1+dx, dx) 35 35 nx = len(xvec) 36 36 y0 = -3000. # min limit of grid 37 37 y1 = 3000. # max limit of grid 38 dy = 100.38 dy = 40. # y resolution of grid 39 39 yvec = np.arange(y0, y1+dy, dy) 40 40 ny = len(yvec)
Note: See TracChangeset
for help on using the changeset viewer.