Changeset 44
- Timestamp:
- 07/23/14 19:15:22 (10 years ago)
- Location:
- trunk/src/scripts_Laura/ARCTIC/Travail_CEN
- Files:
-
- 1 added
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/calcul_spec_lamb_nadir.py
r41 r44 26 26 for imo in range (0, M): 27 27 ####### open .dat file to stack data ####### 28 data_spec_lamb = open ('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA 23.dat', 'a')28 data_spec_lamb = open ('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA89.dat', 'a') 29 29 ############################################ 30 30 print 'read file' + month[imo] … … 69 69 fichier.close() 70 70 for i in range (0, nbtotal): 71 if (emis[ 0, i] < 0.):72 emis[ 0, i] = NaN73 if (emis[ 0, i] > 1.):74 emis[ 0, i] = 0.9971 if (emis[3, i] < 0.): 72 emis[3, i] = NaN 73 if (emis[3, i] > 1.): 74 emis[3, i] = 0.99 75 75 ############################ 76 76 # definition of parameters # … … 83 83 zen_angle = zen[bblat][bbzen] # zenith angle 84 84 field_view = fov[bblat][bbzen] # field of view 85 e = emis[ 0, :][bblat][bbzen] # emissivity85 e = emis[3, :][bblat][bbzen] # emissivity 86 86 t_surf = ts[bblat][bbzen] # surface temperature 87 t_up = tup[ 0, :][bblat][bbzen] # tup88 t_dn = tdn[ 0, :][bblat][bbzen] # tdown89 transm = trans[ 0, :][bblat][bbzen] # transmissivity90 t_b = tb[ 0, :][bblat][bbzen] # brightness temperature87 t_up = tup[3, :][bblat][bbzen] # tup 88 t_dn = tdn[3, :][bblat][bbzen] # tdown 89 transm = trans[3, :][bblat][bbzen] # transmissivity 90 t_b = tb[3, :][bblat][bbzen] # brightness temperature 91 91 L = len(llon) 92 92 ####################################### -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/discrim_rate_espec-elamb_espec.py
r40 r44 14 14 import map_ffgrid 15 15 from matplotlib import colors 16 import map_cartesian_grid 16 17 17 18 … … 24 25 25 26 26 ################### 27 # grid parameters # 28 ################### 29 dx=1. 30 dy=0.5 31 x0, x1 = -180., 180. 32 y0, y1 = 30., 90. 33 xvec = np.arange(x0, x1 + dx, dx) 34 yvec = np.arange(y0, y1 + dy, dy) 35 nx = len(xvec) 27 28 ######################## 29 # grid characteristics # 30 ######################## 31 x0 = -3000. # min limit of grid 32 x1 = 2500. # max limit of grid 33 dx = 40. 34 xvec = np.arange(x0, x1+dx, dx) 35 nx = len(xvec) 36 y0 = -3000. # min limit of grid 37 y1 = 3000. # max limit of grid 38 dy = 40. 39 yvec = np.arange(y0, y1+dy, dy) 36 40 ny = len(yvec) 37 38 39 41 40 42 … … 44 46 ############################################ 45 47 sl_rate = np.zeros([M, 31, ny, nx], float) 46 type_map = np.zeros([M, ny, nx], float)47 48 for imo in range (0, M): 48 49 print 'month: ' + month[imo] 49 fichier = Dataset('/mma/hermozol/Documents/Data/daily_OSISAF_SPEC_LAMB/daily_OSISAF_emis_spec_lamb_AMSUB_ffgrid_' + str(month[imo]) + '_2009.nc', 'r', format='NETCDF3_CLASSIC') 50 lon = fichier.variables['longitude'][:] 51 lat = fichier.variables['latitude'][:] 50 ### AMSUA23 ### 51 print 'open data file' 52 fichier = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA89_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 53 xdist = fichier.variables['longitude'][:] 54 ydist = fichier.variables['latitude'][:] 52 55 day = fichier.variables['days'][:] 53 osi_type = fichier.variables['osi_ice_type'][:] 54 type_map[imo, :, :] = osi_type[0, :, :] 55 amsub_spec = fichier.variables['spec_emis'][:] 56 amsub_lamb = fichier.variables['lamb_emis'][:] 56 emis_spec = fichier.variables['e_spec'][:] 57 emis_lamb = fichier.variables['e_lamb'][:] 57 58 fichier.close() 58 59 print 'calculate rate' 59 60 for ijr in range (0, month_day[imo]): 60 for ilat in range (0, len(lat)): 61 for ilon in range (0, len(lon)): 62 if (amsub_spec[ijr, ilat, ilon] == 0.): 63 amsub_spec[ijr, ilat, ilon] = NaN 64 if (amsub_lamb[ijr, ilat, ilon] == 0.): 65 amsub_lamb[ijr, ilat, ilon] = NaN 66 sl_rate[imo, ijr, ilat, ilon] = ((amsub_lamb[ijr, ilat, ilon] - amsub_spec[ijr, ilat, ilon]) / amsub_spec[ijr, ilat, ilon]) * 100. 61 for ilat in range (0, len(ydist)): 62 for ilon in range (0, len(xdist)): 63 sl_rate[imo, ijr, ilat, ilon] = ((emis_lamb[ilat, ilon, ijr] - emis_spec[ilat, ilon, ijr]) / emis_spec[ilat, ilon, ijr]) * 100. 67 64 68 65 … … 72 69 # compute monthly means of emissivity ratio # 73 70 ############################################# 71 print 'compute monthly mean of emissivity ratio' 74 72 monthly_ratio = np.zeros([M, ny, nx], float) 75 73 for imo in range (0, M): 76 for ilon in range (0, len(lon)): 77 for ilat in range (0, len(lat)): 74 print 'month ' + str(month[imo]) 75 for ilon in range (0, len(xdist)): 76 for ilat in range (0, len(ydist)): 78 77 monthly_ratio[imo, ilat, ilon] = mean(sl_rate[imo, 0 : month_day[imo], ilat, ilon][nonzero(isnan(sl_rate[imo, 0 : month_day[imo], ilat, ilon]) == False)]) 79 78 80 79 80 81 82 81 83 82 84 ######################################### 83 85 # map monthly means of emissivity ratio # 84 86 ######################################### 87 print 'start mapping' 88 ## parameters of coast coordinates for mapping ## 89 plt.ion() 90 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 91 x_coast = x_ind 92 y_coast = y_ind 93 z_coast = z_ind 94 colormap = cm.s3pcpn_l_r 95 ## start mapping ratio ## 85 96 for imo in range (0, M): 86 map_ffgrid.draw_map_l (lon, lat, 0., 10., 0.1, monthly_ratio[imo, :, :], cm.s3pcpn_l_r, 'emissivity ratio [(LAMB - SPEC)/SPEC]*100') 87 title(month[imo] + ' 2009 - AMSUB') 88 plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/space_evolution/EMIS/emis_rate_vs_ice_type/emis_rate_spec-lamb_AMSUB_' + month[imo] + '2009.png') 89 #map_ffgrid.draw_map_l (lon, lat, 0, 5, 1, type_map[imo, :, :], colors.ListedColormap(['1.', '0.25', '0.50', '0.75']), 'Ice Type') 90 #title(month[imo] + ' 2009 - OSISAF') 91 #plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/space_evolution/EMIS/emis_rate_vs_ice_type/ice_type_OSISAF_01' + month[imo] + '2009.png') 97 print 'map month ' +str(month[imo]) 98 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, monthly_ratio[imo, :, :], 0., 10., 0.1, colormap, 'emissivity ratio [(LAMB - SPEC)/SPEC]*100') 99 title(month[imo] + ' 2009 - AMSUA 89GHz') 100 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_ratio/emis_ratio_spec-lamb_AMSUA89_' + month[imo] + '2009.png') 101 102 103 92 104 93 105 … … 96 108 # stack monthly means of emissivity ratio in netcdf files # 97 109 ########################################################### 110 print 'start stacking' 111 ## AMSUA ratio ## 98 112 for imo in range (0, M): 99 print ' file for ' + month[imo]100 rootgrp = Dataset('/ mma/hermozol/Documents/Data/monthly_GLACE/gridded_data/grid_monthly_lamb-spec_ratio_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC')101 rootgrp.createDimension('longitude', len( lon))102 rootgrp.createDimension('latitude', len( lat))113 print 'stack in file month ' + str(month[imo]) 114 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 115 rootgrp.createDimension('longitude', len(xdist)) 116 rootgrp.createDimension('latitude', len(ydist)) 103 117 nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 104 118 nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 105 119 nc_ratio = rootgrp.createVariable('emis_ratio', 'f', ('latitude', 'longitude')) 106 nc_lon[:] = lon107 nc_lat[:] = lat120 nc_lon[:] = xdist 121 nc_lat[:] = ydist 108 122 nc_ratio[:] = monthly_ratio[imo, :, :] 109 123 rootgrp.close() 110 111 112 113 114 115 116 117 118 119 120 121 122 123 #############################124 # delimitation ice - no_ice #125 #############################126 ice_zone = np.zeros([M, ny, nx], float)127 for imo in range (0, M):128 print 'month: ' + month[imo]129 for ilat in range (0, len(lat)):130 for ilon in range (0, len(lon)):131 if (isnan(monthly_ratio[imo, ilat, ilon]) == True):132 ice_zone[imo, ilat, ilon] = NaN133 else:134 if (monthly_ratio[imo, ilat, ilon] >= 3):135 ice_zone[imo, ilat, ilon] = 0.136 else:137 ice_zone[imo, ilat, ilon] = 1.138 # map_ffgrid.draw_map_l (lon, lat, 0., 1.5, 0.5, ice_zone[imo, :, :], colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : ratio < 3.5 %)')139 # title(month[imo] + ' 2009 - AMSUB')140 # plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/space_evolution/EMIS/AMSUB_emis_class/AMSUB_emis_classif_ratio-inf-3_' + month[imo] + '2009.png')141 142 143 144 ##########################################145 # calculation of number of pixels of ice #146 ##########################################147 Rt = 6372.148 149 dist_along_lon = np.zeros([len(lat), len(lon)], float)150 for ilat in range (0, len(lat) - 1):151 for ilon in range (0, len(lon)):152 # calculate distances along each longitude slice153 arc_along_lon = distance_lon_lat.distance_on_unit_sphere(lat[ilat], lon[ilon], lat[ilat + 1], lon[ilon])154 dist_along_lon[ilat, ilon] = arc_along_lon * Rt155 156 dist_along_lat = np.zeros([len(lat), len(lon)], float)157 for ilon in range (0, len(lon) - 1):158 for ilat in range (0, len(lat)):159 # calculate distances along each latitude slice160 arc_along_lat = distance_lon_lat.distance_on_unit_sphere(lat[ilat], lon[ilon], lat[ilat], lon[ilon + 1])161 dist_along_lat[ilat, ilon] = arc_along_lat * Rt162 163 dist_along_lon = 55.164 pix_area = np.zeros([M, len(lat), len(lon)], float)165 total_ice_area = np.zeros([M], float)166 for imo in range (0, M):167 for ilon in range (0, len(lon)):168 for ilat in range (0, len(lat)):169 if (ice_zone[imo, ilat, ilon] == 0.):170 pix_area[imo, ilat, ilon] = NaN171 else:172 if (isnan(ice_zone[imo, ilat, ilon]) == True):173 pix_area[imo, ilat, ilon] = NaN174 else:175 pix_area[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon]176 total_ice_area[imo] = np.sum(reshape(pix_area[imo, :, :], size(pix_area[imo, :, :]))[nonzero(isnan(reshape(pix_area[imo, :, :], size(pix_area[imo, :, :]))) == False)])177 178 179 figure()180 plot(total_ice_area, '.b')181 map_ffgrid.draw_map_l (lon, lat, 0., 5300., 10., pix_area, cm.s3pcpn_l_r, 'dist along lon')182 183 184 185 186 187 188 189 190 191 192 ######## tests ###########193 np.transpose(dist_along_lon)194 195 ice_pixels = np.zeros([M], float)196 for imo in range (0, M):197 ice_vec = reshape(ice_zone[imo, :, :], size(ice_zone[imo, :, :]))198 ice_pixels[imo] = len(np.where(ice_vec == 1.)[0])199 200 201 202 radius = np.zeros([len(lat)], float)203 dx_km = np.zeros([len(lat)], float)204 for ilat in range (0, len(lat)):205 radius[ilat] = Rt * cos(lat[ilat] * c) # calculate radius of earth at a given latitude206 dx_km[ilat] = c * dx * radius[ilat] # calculate distance in km of longitude step dx (dx = 1.)207 208 209 210 211 212 213 214 215 216 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/ice_delimit_AMSUB_data.py
r40 r44 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 cartesian10 import arctic_map # function to regrid coast limits 11 import cartesian_grid_test # function to convert grid from polar to cartesian 12 12 import scipy.special 13 13 import ffgrid2 14 14 import map_ffgrid 15 15 from matplotlib import colors 16 import distance_lon_lat 17 16 from matplotlib.font_manager import FontProperties 17 import map_cartesian_grid 18 18 19 19 MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) … … 23 23 24 24 25 26 ################### 27 # grid parameters # 28 ################### 29 dx=1. 30 dy=0.5 31 x0, x1 = -180., 180. 32 y0, y1 = 30., 90. 33 xvec = np.arange(x0, x1 + dx, dx) 34 yvec = np.arange(y0, y1 + dy, dy) 35 nx = len(xvec) 25 ######################## 26 # grid characteristics # 27 ######################## 28 x0 = -3000. # min limit of grid 29 x1 = 2500. # max limit of grid 30 dx = 40. 31 xvec = np.arange(x0, x1+dx, dx) 32 nx = len(xvec) 33 y0 = -3000. # min limit of grid 34 y1 = 3000. # max limit of grid 35 dy = 40. 36 yvec = np.arange(y0, y1+dy, dy) 36 37 ny = len(yvec) 37 38 … … 41 42 # read NETCDF files # 42 43 ##################### 43 a = 361 * 51 # size of the grid with lats > 65 deg (lon * lat) 44 SP = np.zeros([M, 51, 361], float) 45 LB = np.zeros([M, 51, 361], float) 46 ER = np.zeros([M, 51, 361], float) 47 r_vec = np.zeros([M, a], float) 48 s_vec = np.zeros([M, a], float) 49 l_vec = np.zeros([M, a], float) 50 sl_vec = np.zeros([M, a], float) 51 for imo in range (0, M): 52 print 'file for ' + month[imo] 53 ## read netcdf file of emissivity ratio data ## 54 print 'ratio file' 55 ratio_file = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/grid_monthly_lamb-spec_ratio_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 56 lon_ratio = ratio_file.variables['longitude'][:] 57 lat_ratio = ratio_file.variables['latitude'][:] 58 emis_ratio = ratio_file.variables['emis_ratio'][:] 59 ratio_file.close() 60 ## read netcdf file of emissivity LAMB and SPEC data ## 61 print 'SPEC LAMB file' 62 spec_lamb_file = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/grid_monthly_data_lamb_spec_near_nadir_AMSUB_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 63 lon_sl = spec_lamb_file.variables['longitude'][:] 64 lat_sl = spec_lamb_file.variables['latitude'][:] 65 e_spec = spec_lamb_file.variables['e_spec'][:] 66 e_lamb = spec_lamb_file.variables['e_lamb'][:] 67 e_diff = spec_lamb_file.variables['e_spec_lamb'][:] 68 spec_lamb_file.close() 69 ## keep same dimension between ration and emissivity ## 70 dim_equal = np.where(lat_ratio >= 65.)[0] 71 ldim = len(dim_equal) 72 new_lat = lat_ratio[dim_equal[0] : dim_equal[ldim-1] + 1] # only keep latitude > 65 deg 73 ratio = emis_ratio[dim_equal[0] : dim_equal[ldim-1] + 1, :]# only keep corresponding emissivity ratio 74 SP[imo, :, :] = e_spec 75 LB[imo, :, :] = e_lamb 76 ER[imo, :, :] = ratio 77 r_vec[imo, :] = reshape(ratio, size(ratio)) 78 s_vec[imo, :] = reshape(e_spec, size(e_spec)) 79 l_vec[imo, :] = reshape(e_lamb, size(e_lamb)) 80 sl_vec[imo, :] = reshape(e_diff, size(e_diff)) 81 82 83 84 85 nan_index = np.zeros([M], float) 86 for imo in range (0, M): 87 nan_index[imo] = np.where(isnan(SP[imo, :, 180]) == True)[0][0] - 1 # select index of latitude vector from which data = NaN because of selection of near nadir angles only // above this index (excluded), all NaN values 88 89 90 91 # map test 92 map_ffgrid.draw_map_l (lon_sl, lat_sl, 0.6, 1.02, 0.01, SP[0, :, 180], cm.s3pcpn_l_r, 'emissivity SPEC') 93 94 95 96 ############################################################ 97 # distribution of emissivity rate emissivity spec and lamb # 98 ############################################################ 99 ## histogram plots 100 101 ##### emissivity rate #### 44 emis_spec_month = np.zeros([M, ny, nx], float) 45 emis_lamb_month = np.zeros([M, ny, nx], float) 46 ratio_emis = np.zeros([M, ny, nx], float) 47 ratio_emis_vec = np.zeros([M, ny * nx], float) 48 spec_emis_vec = np.zeros([M, ny * nx], float) 49 lamb_emis_vec = np.zeros([M, ny * nx], float) 50 for imo in range (0, M): 51 print 'month: ' + month[imo] 52 ### emissivity ### 53 print 'open emissivity file' 54 fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 55 xdist = fichier_emis.variables['longitude'][:] 56 ydist = fichier_emis.variables['latitude'][:] 57 day = fichier_emis.variables['days'][:] 58 emis_spec = fichier_emis.variables['e_spec'][:] 59 emis_lamb = fichier_emis.variables['e_lamb'][:] 60 fichier_emis.close() 61 ##### calculate monthly mean of emis ##### 62 for ilat in range(0, len(ydist)): 63 for ilon in range (0, len(xdist)): 64 emis_spec_month[imo, ilat, ilon] = mean(emis_spec[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_spec[ilat, ilon, 0 : month_day[imo]]) == False)]) 65 emis_lamb_month[imo, ilat, ilon] = mean(emis_lamb[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(emis_lamb[ilat, ilon, 0 : month_day[imo]]) == False)]) 66 ### monthly emissivity ratio ### 67 print 'open emissivity ratio file' 68 fichier_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA23_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC') 69 xdist = fichier_ratio.variables['longitude'][:] 70 ydist = fichier_ratio.variables['latitude'][:] 71 ratio_emis[imo, :, :] = fichier_ratio.variables['emis_ratio'][:] 72 fichier_ratio.close() 73 ### reshape in monthly vectors all data in 2D arrays ### 74 print 'reshape data' 75 ratio_emis_vec[imo, :] = reshape(ratio_emis[imo, :, :], size(ratio_emis[imo, :, :])) 76 spec_emis_vec[imo, :] = reshape(emis_spec_month[imo, :, :], size(emis_spec_month[imo, :, :])) 77 lamb_emis_vec[imo, :] = reshape(emis_lamb_month[imo, :, :], size(emis_lamb_month[imo, :, :])) 78 79 80 ''' 81 # test: 82 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 83 x_coast = x_ind 84 y_coast = y_ind 85 z_coast = z_ind 86 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, emis_spec[:, :, 0], 0.4, 1.2, 0.01, cm.s3pcpn_l_r, 'test') 87 ''' 88 89 102 90 c = np.array(['r', 'b', 'c', 'm', 'y', 'g']) 91 ################################### 92 # distribution of SPEC emissivity # 93 ################################### 94 hist_vals_spec = np.zeros([M, 50], float) 95 corresp_emis_spec = np.zeros([M, 50], float) 96 for imo in range (0, M): 97 hist_vals_spec[imo, :] = hist(spec_emis_vec[imo, :][nonzero(isnan(spec_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 98 for ibin in range (0, 50): 99 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]) 100 101 ## plot first six months of spec emissivity histograms ## 102 ion() 103 103 figure() 104 104 for imo in range (0, 6): 105 hist(r_vec[imo, :][nonzero(isnan(r_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step', color = str(c[imo]), label = str(month[imo]))[0] 106 107 legend() 108 grid() 109 xlim(0, 10) 110 xticks(np.arange(0., 12., 1.), np.array(['0%', '1%', '2%', '3%', '4%', '5%', '6%', '7%', '8%', '9%', '10%'])) 105 plot(corresp_emis_spec[imo], hist_vals_spec[imo], c = str(c[imo]), label = str(month[imo])) 106 107 grid() 108 xlim(corresp_emis_spec.min() - 0.02, 1.02) 109 xlabel('emissivity SPEC') 110 ylabel('frequency of occurence') 111 fontP = FontProperties() 112 fontP.set_size('small') 113 legend(loc = 2, prop = fontP) 114 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA23_JANUARY-JUNE_2009.png') 115 ## plot six following months of spec emissivity histograms ## 116 figure() 117 for imo in range (6, M): 118 plot(corresp_emis_spec[imo], hist_vals_spec[imo], c = str(c[imo - 6]), label = str(month[imo])) 119 120 grid() 121 xlim(corresp_emis_spec.min() - 0.02, 1.02) 122 xlabel('emissivity SPEC') 123 ylabel('frequency of occurence') 124 legend(loc = 9, prop = fontP) 125 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_spec_AMSUA23_JULY-DECEMBER_2009.png') 126 127 # find the emissivity value corresponding to the threshold of ice/no_ice limit 128 emis_lim_spec = np.zeros([M], float) 129 for imo in range (0, M): 130 print 'month ' + str(month[imo]) 131 bb = np.where((corresp_emis_spec[imo, :] < 0.6))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 132 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 133 cc = np.where(hist_vals_spec[imo, bb] == max(hist_vals_spec[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 134 dd = np.where(aa > cc[0])[0] 135 emis_lim_spec[imo] = corresp_emis_spec[imo, bb][aa][dd[0]] 136 137 # plot monthly evolution of emissivity spec threshold used for delimitation 138 figure() 139 plot(np.arange(0, M, 1), emis_lim_spec, 'b-+', label = 'emis SPEC') 140 ylabel('emissivity SPEC threshold - AMSUA 23GHz') 141 grid() 142 xticks(np.arange(0, M, 1), month, rotation = 25) 143 legend(loc = 2) 144 xlim(-1, M) 145 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_SPEC_AMSUA23_2009.png') 146 147 148 ################################### 149 # distribution of LAMB emissivity # 150 ################################### 151 hist_vals_lamb = np.zeros([M, 50], float) 152 corresp_emis_lamb = np.zeros([M, 50], float) 153 for imo in range (0, M): 154 hist_vals_lamb[imo, :] = hist(lamb_emis_vec[imo, :][nonzero(isnan(lamb_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 155 for ibin in range (0, 50): 156 corresp_emis_lamb[imo, ibin] = mean(hist(lamb_emis_vec[imo, :][nonzero(isnan(lamb_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 157 158 ## plot first six months of spec emissivity histograms ## 159 figure() 160 for imo in range (0, 6): 161 plot(corresp_emis_lamb[imo], hist_vals_lamb[imo], c = str(c[imo]), label = str(month[imo])) 162 163 grid() 164 xlim(corresp_emis_lamb.min() - 0.02, 1.02) 165 xlabel('emissivity LAMB') 166 ylabel('frequency of occurence') 167 legend(loc = 2, prop = fontP) 168 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA23_JANUARY-JUNE_2009.png') 169 ## plot six following months of spec emissivity histograms ## 170 figure() 171 for imo in range (6, M): 172 plot(corresp_emis_lamb[imo], hist_vals_lamb[imo], c = str(c[imo - 6]), label = str(month[imo])) 173 174 grid() 175 xlim(corresp_emis_lamb.min() - 0.02, 1.02) 176 xlabel('emissivity LAMB') 177 ylabel('frequency of occurence') 178 legend(loc = 9, prop = fontP) 179 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_lamb_AMSUA23_JULY-DECEMBER_2009.png') 180 181 # find the emissivity value corresponding to the threshold of ice/no_ice limit 182 emis_lim_lamb = np.zeros([M], float) 183 for imo in range (0, M): 184 print 'mont ' + str(month[imo]) 185 bb = np.where((corresp_emis_lamb[imo, :] < 0.6))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 186 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 187 cc = np.where(hist_vals_lamb[imo, bb] == max(hist_vals_lamb[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 188 dd = np.where(aa > cc[0])[0] 189 emis_lim_lamb[imo] = corresp_emis_lamb[imo, bb][aa][dd[0]] 190 191 # plot monthly evolution of emissivity spec threshold used for delimitation 192 figure() 193 plot(np.arange(0, M, 1), emis_lim_lamb, 'b-+', label = 'emis LAMB') 194 ylabel('emissivity LAMB threshold - AMSUB 23GHz') 195 grid() 196 xticks(np.arange(0, M, 1), month, rotation = 25) 197 legend(loc = 2) 198 xlim(-1, M) 199 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_LAMB_AMSUA23_2009.png') 200 201 202 203 ################################### 204 # distribution of emissivity rate # 205 ################################### 206 hist_vals_rate = np.zeros([M, 50], float) 207 corresp_emis_rate = np.zeros([M, 50], float) 208 for imo in range (0, M): 209 hist_vals_rate[imo, :] = hist(ratio_emis_vec[imo, :][nonzero(isnan(ratio_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 210 for ibin in range (0, 50): 211 corresp_emis_rate[imo, ibin] = mean(hist(ratio_emis_vec[imo, :][nonzero(isnan(ratio_emis_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 212 213 ## plot first six months of spec emissivity histograms ## 214 figure() 215 for imo in range (0, 6): 216 plot(corresp_emis_rate[imo], hist_vals_rate[imo], c = str(c[imo]), label = str(month[imo])) 217 218 grid() 219 xlim(corresp_emis_rate.min() - 0.02, corresp_emis_rate.max() + 0.02) 111 220 xlabel('emissivity rate (%)') 112 ylabel('frequency of appearence') 113 plt.savefig('/mma/hermozol/Documents/figure_output/ice_class_AMSUB/emiss_rate_AMSUB_JANUARY-JUNE_2009.png') 114 221 ylabel('frequency of occurence') 222 legend(prop = fontP) 223 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA23_JANUARY-JUNE_2009.png') 224 ## plot six following months of spec emissivity histograms ## 115 225 figure() 116 226 for imo in range (6, M): 117 hist(r_vec[imo, :][nonzero(isnan(r_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step', color = str(c[imo - 6]), label = str(month[imo])) 118 119 legend() 120 grid() 121 xlim(0, 10) 122 xticks(np.arange(0., 12., 1.), np.array(['0%', '1%', '2%', '3%', '4%', '5%', '6%', '7%', '8%', '9%', '10%'])) 227 plot(corresp_emis_rate[imo], hist_vals_rate[imo], c = str(c[imo - 6]), label = str(month[imo])) 228 229 grid() 230 xlim(corresp_emis_rate.min() - 0.02, corresp_emis_rate.max() + 0.02) 123 231 xlabel('emissivity rate (%)') 124 ylabel('frequency of appearence') 125 plt.savefig('/mma/hermozol/Documents/figure_output/ice_class_AMSUB/emiss_rate_AMSUB_JULY-DECEMBER_2009.png') 126 127 ##### SPEC emissivity #### 128 hist_vals = np.zeros([M, 50], float) 129 corresp_emis = np.zeros([M, 50], float) 130 for imo in range (0, M): 131 hist_vals[imo, :] = hist(s_vec[imo, :][nonzero(isnan(s_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 132 for ibin in range (0, 50): 133 corresp_emis[imo, ibin] = mean(hist(s_vec[imo, :][nonzero(isnan(s_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 134 135 figure() 136 for imo in range (0, 6): 137 plot(corresp_emis[imo], hist_vals[imo], c = str(c[imo]), label = str(month[imo])) 138 139 grid() 140 xlim(0.5, 1.02) 141 xlabel('emissivity SPEC') 142 ylabel('frequency of occurence') 232 ylabel('frequency of occurence') 233 legend(loc = 9, prop = fontP) 234 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emiss_rate_AMSUA23_JULY-DECEMBER_2009.png') 235 236 # no definition of threshold for ICE / NO_ICE delimitation because no apparent distinction signal // use of this signal for internal ice classification 237 ''' 238 # find the emissivity value corresponding to the threshold of ice/no_ice limit 239 emis_lim_rate = np.zeros([M], float) 240 for imo in range (0, M): 241 bb = np.where((corresp_emis_rate[imo, :] < 0.6))[0] # only consider low emissivity values (open sea) // CHANGE EMISSIVITY VALUE ACCORDING TO FREQUENCY 242 aa = np.where(hist_vals_rate[imo, bb] < max(hist_vals_rate[imo, bb]) / 4.)[0] # of these latter values, only consider emissivity values lower than 1/4*peak emissivity 243 cc = np.where(hist_vals_rate[imo, bb] == max(hist_vals_rate[imo, bb]))[0] # which emissivity index corresponds to peak emissivity 244 dd = np.where(aa > cc[0])[0] 245 emis_lim_rate[imo] = corresp_emis_rate[imo, bb][aa][dd[0]] 246 247 # plot monthly evolution of emissivity spec threshold used for delimitation 248 figure() 249 plot(np.arange(0, M, 1), emis_lim_rate, 'b-+', label = 'emis LAMB') 250 ylabel('emissivity rate threshold - AMSUB 89GHz') 251 grid() 252 xticks(np.arange(0, M, 1), month, rotation = 25) 143 253 legend(loc = 2) 144 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_spec_AMSUB_JANUARY-JUNE_2009.png') 145 146 147 figure() 148 for imo in range (6, M): 149 plot(corresp_emis[imo], hist_vals[imo], c = str(c[imo - 6]), label = str(month[imo])) 150 151 grid() 152 xlim(0.5, 1.02) 153 xlabel('emissivity SPEC') 154 ylabel('frequency of occurence') 155 legend(loc = 2) 156 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_spec_AMSUB_JULY-DECEMBER_2009.png') 157 158 159 # find the emissivity value corresponding to the threshold of ice/no_ice limit 160 # after many tests, we choose an emissivity SPEC value which frequency of occurence is minimal after the peak of frequency around 0.7 161 emis_lim = np.zeros([M], float) 162 for imo in range (0, M): 163 bb = np.where((corresp_emis[imo, :] > 0.7) & (corresp_emis[imo, :] < 0.8))[0] 164 aa = np.where(hist_vals[imo, bb] == min(hist_vals[imo, bb]))[0] 165 emis_lim[imo] = corresp_emis[imo, bb][aa] 166 167 figure() 168 plot(np.arange(0, M, 1), emis_lim, 'b-+', label = 'emis SPEC') 169 ylabel('emissivity SPEC threshold - AMSUB 89GHz') 170 grid() 171 xticks(np.arange(0, M, 1), month, rotation = 25) 172 legend(loc = 2) 173 xlim(-1, M) 174 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emis_lim_frequ_SPEC_AMSUB_2009.png') 175 176 177 ##### LAMB emissivity #### 178 hist_vals_l = np.zeros([M, 50], float) 179 corresp_emis_l = np.zeros([M, 50], float) 180 for imo in range (0, M): 181 hist_vals_l[imo, :] = hist(l_vec[imo, :][nonzero(isnan(l_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[0] 182 for ibin in range (0, 50): 183 corresp_emis_l[imo, ibin] = mean(hist(l_vec[imo, :][nonzero(isnan(l_vec[imo,:]) == False)], bins = 50, normed = True, histtype='step')[1][ibin : ibin + 2]) 184 185 figure() 186 for imo in range (0, 6): 187 plot(corresp_emis_l[imo], hist_vals_l[imo], c = str(c[imo]), label = str(month[imo])) 188 189 grid() 190 xlim(0.5, 1.02) 191 xlabel('emissivity LAMB') 192 ylabel('frequency of occurence') 193 legend(loc = 2) 194 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_lamb_AMSUB_JANUARY-JUNE_2009.png') 195 196 197 figure() 198 for imo in range (6, M): 199 plot(corresp_emis_l[imo], hist_vals_l[imo], c = str(c[imo - 6]), label = str(month[imo])) 200 201 grid() 202 xlim(0.5, 1.02) 203 xlabel('emissivity LAMB') 204 ylabel('frequency of occurence') 205 legend(loc = 2) 206 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emiss_lamb_AMSUB_JULY-DECEMBER_2009.png') 207 208 209 # find the emissivity value corresponding to the threshold of ice/no_ice limit 210 # after many tests, we choose an emissivity LAMB value which frequency of occurence is minimal after the peak of frequency around 0.73 211 emis_lim_l = np.zeros([M], float) 212 for imo in range (0, M): 213 bb = np.where((corresp_emis_l[imo, :] > 0.72) & (corresp_emis_l[imo, :] < 0.8))[0] 214 aa = np.where(hist_vals_l[imo, bb] == min(hist_vals_l[imo, bb]))[0] 215 emis_lim_l[imo] = corresp_emis_l[imo, bb][aa] 216 217 figure() 218 plot(np.arange(0, M, 1), emis_lim_l, 'b-+', label = 'emis LAMB') 219 ylabel('emissivity LAMB threshold - AMSUB 89GHz') 220 grid() 221 xticks(np.arange(0, M, 1), month, rotation = 25) 222 legend(loc = 2) 223 xlim(-1, M) 224 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUB/emis_lim_frequ_LAMB_AMSUB_2009.png') 225 226 227 228 ####################################################### 229 # delimitation ice - no_ice with different conditions # 230 ####################################################### 231 ## condition : threshold of SPEC emissivity ('SP') = f[month] ==> take out all emissivities <= f[month] (open sea) // f[month] depends on the month, around 0.7 232 ice_zone_l = np.zeros([M, 51, 361], float) 233 for imo in range (0, M): 234 print 'month: ' + month[imo] 235 for ilat in range (0, 51): 236 for ilon in range (0, 361): 237 if (isnan(LB[imo, ilat, ilon]) == True): 238 ice_zone_l[imo, ilat, ilon] = NaN 254 xlim(-1, M) 255 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_rate_AMSUA89_2009.png') 256 ''' 257 258 259 ############################################################################## 260 # plot comparison of emissivity thresholds between spec and lamb assumptions # 261 ############################################################################## 262 figure() 263 plot(np.arange(0, M, 1), emis_lim_spec, 'b-+', label = 'emis SPEC') 264 plot(np.arange(0, M, 1), emis_lim_lamb, 'g-+', label = 'emis LAMB') 265 ylabel('emissivity threshold - AMSUA 23GHz') 266 grid() 267 xticks(np.arange(0, M, 1), month, rotation = 25) 268 legend(loc = 2, prop = fontP) 269 xlim(-1, M) 270 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/emis_lim_frequ_spec_and_lamb_AMSUA23_2009.png') 271 272 273 274 275 ################################################################ 276 # delimitation ice - no_ice with emissivity or ratio threshold # 277 ################################################################ 278 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 279 x_coast = x_ind 280 y_coast = y_ind 281 z_coast = z_ind 282 ice_zone = np.zeros([M, 151, 139], float) 283 for imo in range (0, M): 284 print 'month: ' + str(month[imo]) 285 for ilat in range (0, 151): 286 for ilon in range (0, 139): 287 if (isnan(emis_spec_month[imo, ilat, ilon]) == True): 288 ice_zone[imo, ilat, ilon] = NaN 239 289 else: 240 if ( LB[imo, ilat, ilon] <= emis_lim_l[imo]): # use the monthly SPEC emissivity threshold241 ice_zone _l[imo, ilat, ilon] = 0.290 if (emis_spec_month[imo, ilat, ilon] <= emis_lim_spec[imo]): # use the monthly SPEC emissivity threshold 291 ice_zone[imo, ilat, ilon] = 0. 242 292 else: 243 ice_zone_l[imo, ilat, ilon] = 1. 244 map_ffgrid.draw_map_l (lon_sl, lat_sl, -1., 2., 1., ice_zone_l[imo, :, :], colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_LAMB > ' + str(emis_lim_l[imo]) + ')') 245 title(month[imo] + ' 2009 - AMSUB') 246 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/AMSUB_emis_class/AMSUB_emis_classif_elamb_thresh' + '_' + month[imo] + '2009.png') 247 248 249 250 251 252 ########################################## 253 # calculation of number of pixels of ice # 254 ########################################## 255 Rt = 6372. 256 257 '''dist_along_lon = np.zeros([len(lat_sl), len(lon_sl)], float) 258 for ilat in range (0, len(lat_sl) - 1): # leave last lat index (lat = 90.) 259 for ilon in range (0, len(lon_sl)): 260 # calculate distances along each longitude slice 261 arc_along_lon = distance_lon_lat.distance_on_unit_sphere(lat_sl[ilat], lon_sl[ilon], lat_sl[ilat + 1], lon_sl[ilon]) 262 dist_along_lon[ilat, ilon] = arc_along_lon * Rt''' 263 264 dist_along_lon = 55.60618997 265 dist_along_lat = np.zeros([len(lat_sl), len(lon_sl)], float) 266 for ilon in range (0, len(lon_sl) - 1): # leave last lon index (lon = 180.) 267 for ilat in range (0, len(lat_sl)): 268 # calculate distances along each latitude slice 269 arc_along_lat = distance_lon_lat.distance_on_unit_sphere(lat_sl[ilat], lon_sl[ilon], lat_sl[ilat], lon_sl[ilon + 1]) 270 dist_along_lat[ilat, ilon] = arc_along_lat * Rt 271 272 # deal with last lon index (lon = 180.) 273 for ilat in range (0, len(lat_sl)): 274 dist_along_lat[ilat, 360] = dist_along_lat[ilat, 359] 275 276 277 ## calculate the area of each pixel AMSUB ## 278 pix_area_l = np.zeros([M, len(lat_sl), len(lon_sl)], float) 279 total_ice_area_l = np.zeros([M], float) 280 for imo in range (0, M): 281 for ilon in range (0, len(lon_sl)): 282 for ilat in range (0, len(lat_sl)): 283 if (ice_zone_l[imo, ilat, ilon] == 0.): 284 pix_area_l[imo, ilat, ilon] = NaN 285 else: 286 if ((lat_sl[ilat] <= 83.5) & (isnan(ice_zone_l[imo, ilat, ilon]) == True)): 287 pix_area_l[imo, ilat, ilon] = NaN 288 else: 289 if ((lat_sl[ilat] > 83.5)): 290 pix_area_l[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon] 291 else: 292 pix_area_l[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon] 293 total_ice_area_l[imo] = np.sum(reshape(pix_area_l[imo, :, :], size(pix_area_l[imo, :, :]))[nonzero(isnan(reshape(pix_area_l[imo, :, :], size(pix_area_l[imo, :, :]))) == False)]) 294 map_ffgrid.draw_map_l (lon_sl, lat_sl, 0., 2614., 10., pix_area_l[imo,:,:], cm.s3pcpn_l_r, 'pixel area') 295 title(str(month[imo]) + '2009 - AMSUB') 296 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_pixel_area_AMSUB_' + month[imo] + '2009.png') 297 298 299 300 301 302 303 ################################################## 304 # Read OSISAF ice types to have total ice extent # 305 ################################################## 306 osi_type = np.zeros([M, 31, 51, 361], float) 307 for imo in range (0, M): 308 print 'month ' + month[imo] 309 # osisaf 310 print 'read OSISAF' 311 fichier_osi = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_ffgrid/OSISAF_ice_types_ffgrid_' + month[imo] + '_2009.nc', 'r', format='NETCDF3_CLASSIC') 312 lon_osi = fichier_osi.variables['longitude'][:] 313 lat_osi = fichier_osi.variables['latitude'][:] 314 days_osi = fichier_osi.variables['days'][:] 315 osi_ice_type = fichier_osi.variables['osi_ice_type'][:] 316 osi_type[imo, 0 : month_day[imo], :, :] = osi_ice_type[:, np.where(lat_osi == 65.)[0][0] : len(lat_osi), :] # keep the same dimension in lat (start at 65 deg) 317 fichier_osi.close() 318 319 320 osi_ice = np.zeros([M, 31, 51, 361], float) 321 for imo in range (0, M): 322 for ilon in range (0, len(lon_sl)): 323 for ilat in range (0, len(lat_sl)): 293 ice_zone[imo, ilat, ilon] = 1. 294 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, ice_zone[imo, :, :], -1, 2., 1., colors.ListedColormap(['0.5', 'blue']), 'Ice limit (threshold : emis_SPEC > ' + str(emis_lim_spec[imo])[0:6] + ')') 295 title(str(month[imo]) + ' 2009 - AMSUA 23GHz') 296 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_class_AMSUA/AMSUA23_ice_emis_spec_thresh' + '_' + month[imo] + '2009.png') 297 298 299 300 301 302 ########################### 303 # calculation of ice area # 304 ########################### 305 # nb of pixels near pole = 926 306 pix_area = 40. * 40. 307 nb_pix = np.zeros([M], float) 308 total_ice_area = np.zeros([M], float) 309 for imo in range (0, M): 310 ice = reshape(ice_zone[imo, :, :], size(ice_zone[imo, :, :])) 311 nb_pix[imo] = len(np.where(ice == 1.)[0]) 312 total_ice_area[imo] = (pix_area * nb_pix[imo]) + (926. * pix_area) 313 314 315 spec_ice_area = total_ice_area 316 317 318 figure() 319 plot(lamb_ice_area, 'g', label = 'lamb') 320 plot(spec_ice_area, 'b', label = 'spec') 321 plot(total_monthly_ice_area_osisaf, 'r', label = 'OSISAF') 322 ylabel('total ice area (square km)') 323 grid() 324 xticks(np.arange(0, M, 1), month, rotation = 25) 325 legend(loc = 3, prop = fontP) 326 title('AMSUA 23GHz') 327 xlim(-1, M) 328 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/total_ice_area_AMSUA23_OSISAF_2009.png') 329 330 ########################## 331 # stack of ice area data # 332 ########################## 333 # stack in .dat file OSISAF data 334 print 'start stacking in .dat file' 335 data_osisaf = open ('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_daily_ice_no-ice_2009.dat', 'a') 336 for imo in range (0, M): 337 for ijr in range (0, month_day[imo]): 338 data_osisaf.write(('%(months)10s %(days)2.0f %(total_monthly_ice_area_osisaf)10.5f %(total_ice_area_osisaf)10.5f \n' % { 339 'months':month[imo], 340 'days':np.arange(1, month_day[imo] + 1)[ijr], 341 'total_monthly_ice_area_osisaf':total_monthly_ice_area_osisaf[imo], 342 'total_ice_area_osisaf':total_ice_area_osisaf[imo, ijr], 343 })) 344 345 data_osisaf.close() 346 347 348 # stack in .dat file 349 print 'start stacking in .dat file' 350 data_classif = open ('/net/argos/data/parvati/lahlod/ARCTIC/2009.dat', 'a') 351 for imo in range (0, 50): 352 for ii in range (0, month_day[imo]): 353 data_classif.write(('%(months)10s %(hist_vals_spec)10.5f %(corresp_emis_spec)10.5f %(hist_vals_lamb)10.5f %(orresp_emis_lamb)10.5f %(hist_vals_rate)10.5f %(corresp_emis_rate)10.5f %(emis_lim_spec)10.5f (emis_lim_lamb)10.5f (spec_ice_area)10.5f (lamb_ice_area)10.5f \n' % { 354 'months':month[imo], 355 'hist_vals_spec':hist_vals_spec[imo, ii], 356 'corresp_emis_spec':corresp_emis_spec[imo, ii], 357 'hist_vals_lamb':hist_vals_lamb[imo, ii], 358 'corresp_emis_lamb':corresp_emis_lamb[imo, ii], 359 'hist_vals_rate':hist_vals_rate[imo, ii], 360 'emis_lim_spec':emis_lim_spec[imo], 361 'emis_lim_lamb':emis_lim_lamb[imo], 362 'spec_ice_area':spec_ice_area[imo], 363 'lamb_ice_area':lamb_ice_area[imo], 364 })) 365 366 data_classif.close() 367 368 # stack in netcdf file 369 print 'start stacking' 370 for imo in range (0, M): 371 print 'stack in file month ' + str(month[imo]) 372 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + month[imo] + '2009_AMSUA23.nc', 'w', format='NETCDF3_CLASSIC') 373 rootgrp.createDimension('longitude', len(xdist)) 374 rootgrp.createDimension('latitude', len(ydist)) 375 nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) 376 nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 377 nc_ice = rootgrp.createVariable('ice_area', 'f', ('latitude', 'longitude')) 378 nc_lon[:] = xdist 379 nc_lat[:] = ydist 380 nc_ice[:] = ice_zone[imo, :, :] 381 rootgrp.close() 382 383 384 ############################################################################################################# 385 386 ######################################### 387 # comparison with OSISAF sea ice extent # 388 ######################################### 389 osi_ice = np.zeros([M, 31, 151, 139]) 390 for imo in range (0, M): 391 print 'open OSISAF file month ' + str(month[imo]) 392 fichier_osisaf = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_cartesian_grid/OSISAF_ice_types_cartes_grid_' + str(month[imo]) + '_2009.nc', 'r', format='NETCDF3_CLASSIC') 393 xdist = fichier_osisaf.variables['longitude'][:] 394 ydist = fichier_osisaf.variables['latitude'][:] 395 day = fichier_osisaf.variables['days'][:] 396 osi_type = fichier_osisaf.variables['osi_ice_type'][:] 397 fichier_osisaf.close() 398 for ilon in range (0, len(xdist)): 399 for ilat in range (0, len(ydist)): 324 400 for ijr in range (0, month_day[imo]): 325 if ( osi_type[imo, ijr, ilat, ilon] >= 2.):401 if ((isnan(osi_type[ijr, ilat, ilon]) == True) & (xdist[ilon] >= -312.5) & (xdist[ilon] <= 328.5) & (ydist[ilat] >=-334.3) & (ydist[ilat] <= 338.5)): 326 402 osi_ice[imo, ijr, ilat, ilon] = 1. 327 403 else: 328 if ((lat_sl[ilat] <= 86.) & (isnan(osi_type[imo, ijr, ilat, ilon]) == True)): 329 osi_ice[imo, ijr, ilat, ilon] = NaN 330 else: 331 if (lat_sl[ilat] > 86.): 332 osi_ice[imo, ijr, ilat, ilon] = 1. 333 else: 404 if (osi_type[ijr, ilat, ilon] >= 2.): 405 osi_ice[imo, ijr, ilat, ilon] = 1. 406 else: 407 if (isnan(osi_type[ijr, ilat, ilon]) == True): 334 408 osi_ice[imo, ijr, ilat, ilon] = NaN 335 336 337 338 osi_ice_month = np.zeros([M, 51, 361], float) 339 for imo in range (0, M): 340 for ilon in range (0, len(lon_sl)): 341 for ilat in range (0, len(lat_sl)): 342 how_many_ice_pix = len(np.where(osi_ice[imo, :, ilat, ilon] == 1.)[0]) 343 if (how_many_ice_pix >= 15): 344 osi_ice_month[imo, ilat, ilon] = 1. 345 else: 346 osi_ice_month[imo, ilat, ilon] = NaN 347 348 349 350 ########################################### 351 # calculate the area of each pixel OSISAF # 352 ########################################### 353 pix_area_osi = np.zeros([M, len(lat_sl), len(lon_sl)], float) 354 total_ice_area_osi = np.zeros([M], float) 355 for imo in range (0, M): 356 for ilon in range (0, len(lon_sl)): 357 for ilat in range (0, len(lat_sl)): 358 if (osi_ice_month[imo, ilat, ilon] == 1.): 359 pix_area_osi[imo, ilat, ilon] = dist_along_lon * dist_along_lat[ilat, ilon] 360 else: 361 pix_area_osi[imo, ilat, ilon] = NaN 362 total_ice_area_osi[imo] = np.sum(reshape(pix_area_osi[imo, :, :], size(pix_area_osi[imo, :, :]))[nonzero(isnan(reshape(pix_area_osi[imo, :, :], size(pix_area_osi[imo, :, :]))) == False)]) 363 map_ffgrid.draw_map_l (lon_sl, lat_sl, 0., 2614., 10., pix_area_osi[imo,:,:], cm.s3pcpn_l_r, 'pixel area') 364 title(str(month[imo]) + '2009 - OSISAF') 365 plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/ice_pix_area/ice_pixel_area_OSISAF_' + month[imo] + '2009.png') 366 367 368 369 409 else: 410 osi_ice[imo, ijr, ilat, ilon] = 0. 411 412 413 # test: 414 colormap = colors.ListedColormap(['cyan','blue']) 415 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_ice[10, 0, :, :], -1, 2, 1, colormap, 'ice_type') 416 417 418 ################################## 419 # calculation of OSISAF ice area # 420 ################################## 421 pix_area = 40. * 40. 422 nb_pix_osisaf = np.zeros([M, 31], float) 423 total_ice_area_osisaf = np.zeros([M, 31], float) 424 total_monthly_ice_area_osisaf = np.zeros([M], float) 425 for imo in range (0, M): 426 for ijr in range (0, month_day[imo]): 427 ice = reshape(osi_ice[imo, ijr, :, :], size(osi_ice[imo, ijr, :, :])) 428 nb_pix_osisaf[imo, ijr] = len(np.where(ice == 1.)[0]) 429 total_ice_area_osisaf[imo, ijr] = pix_area * nb_pix_osisaf[imo, ijr] 430 total_monthly_ice_area_osisaf[imo] = mean(total_ice_area_osisaf[imo, 0 : month_day[imo]][nonzero(isnan(total_ice_area_osisaf[imo, 0 : month_day[imo]]) == False)]) + (926. * pix_area) 431 432 433 434 435 436 ################################### 437 # plot time evolution of ice area # 438 ################################### 370 439 figure() 371 440 plot(total_ice_area, '-+b', label = 'AMSUB SPEC') … … 381 450 382 451 383 map_ffgrid.draw_map_l (lon, lat, 0., 5300., 10., pix_area, cm.s3pcpn_l_r, 'dist along lon') 384 385 386 387 388 452 453 454 455 456 457 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_cartesian_gridded_AMSUA23.py
r42 r44 53 53 for imo in range (0, M): 54 54 print 'open file ' + str(month[imo]) 55 file_amsua 23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC')56 lon_amsua = file_amsua 23.variables['longitude'][:]57 lat_amsua = file_amsua 23.variables['latitude'][:]58 days = file_amsua 23.variables['days'][:]59 es_a = file_amsua 23.variables['e_spec'][:]60 el_a = file_amsua 23.variables['e_lamb'][:]61 esl_a = file_amsua 23.variables['e_spec_lamb'][:]62 file_amsua 23.close()55 file_amsua = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') 56 lon_amsua = file_amsua.variables['longitude'][:] 57 lat_amsua = file_amsua.variables['latitude'][:] 58 days = file_amsua.variables['days'][:] 59 es_a = file_amsua.variables['e_spec'][:] 60 el_a = file_amsua.variables['e_lamb'][:] 61 esl_a = file_amsua.variables['e_spec_lamb'][:] 62 file_amsua.close() 63 63 print 'compute monthly mean of data for map' 64 64 for ilon in range (0, len(lon_amsua)): … … 69 69 70 70 71 71 print 'map of data' 72 ion() 73 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() 74 x_coast = x_ind 75 y_coast = y_ind 76 z_coast = z_ind 77 colormap = cm.s3pcpn_l_r 72 78 for imo in range (0, M): 73 x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()74 x_coast = x_ind75 y_coast = y_ind76 z_coast = z_ind77 colormap = cm.s3pcpn_l_r78 79 # emiss SPEC 79 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, es_a_month[imo, :, :], 0. 45, 1.02, 0.01, colormap, 'emissivity SPEC')80 title(month[imo] + ' 2009 - AMSUA 23GHz')81 savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/ emis_spec_AMSUA23_' + month[imo] + '2009.png')80 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, es_a_month[imo, :, :], 0.6, 1.02, 0.01, colormap, 'emissivity SPEC') 81 title(month[imo] + ' 2009 - AMSUA 89GHz') 82 savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec_AMSUA89_' + month[imo] + '2009.png') 82 83 # emis SPEC-LAMB 83 84 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, esl_a_month[imo, :, :], 0., 0.04, 0.001, colormap, 'emissivity LAMB - SPEC') 84 title(month[imo] + ' 2009 - AMSUA 23GHz')85 savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/ emis_spec-lamb_AMSUA23_' + month[imo] + '2009.png')85 title(month[imo] + ' 2009 - AMSUA 89GHz') 86 savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_spec-lamb_AMSUA89_' + month[imo] + '2009.png') 86 87 87 88 88 89 89 ###################90 '''################### 90 91 # Read AMSUB data # 91 92 ################### … … 125 126 map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, esl_b_month[imo, :, :], 0., 0.05, 0.001, colormap, 'emissivity LAMB - SPEC') 126 127 title(month[imo] + ' 2009 - AMSUB 89GHz') 127 savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/emis_spec-lamb_AMSUB89_' + month[imo] + '2009.png') 128 savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/comp_lamb_spec/space_evolution/EMIS/emis_spec-lamb_AMSUB89_' + month[imo] + '2009.png')''' 128 129 129 130 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_spec_lamb_nadir.py
r42 r44 59 59 for imo in range (0, M): 60 60 print 'month: ' + month[imo] 61 fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009 .dat','r')61 fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA89.dat','r') 62 62 numlines = 0 63 63 for line in fichier: numlines += 1 64 64 fichier.close 65 fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009 .dat','r')65 fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA89.dat','r') 66 66 nbtotal = numlines-1 67 67 iligne = 0 … … 185 185 e_sl_100_day = zgrid_output 186 186 esl100[:, :, 0 : month_day[imo], imo] = e_sl_100_day[:, :, :]''' 187 188 189 190 ############################################### 191 # stack gridded data spec lamb in NETCDF file # 192 ############################################### 193 print 'stacking of gridded data' 194 for imo in range (0, M): 195 print 'file for ' + month[imo] 196 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 187 ############################################### 188 # stack gridded data spec lamb in NETCDF file # 189 ############################################### 190 print 'stacking of gridded data' 191 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 197 192 rootgrp.createDimension('longitude', nx) 198 193 rootgrp.createDimension('latitude', ny) -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_spec_lamb_nadir_bis.py
r41 r44 45 45 # grid data from .dat files # 46 46 ############################# 47 tu = np.zeros([ny, nx, 31, M], float)47 '''tu = np.zeros([ny, nx, 31, M], float) 48 48 td = np.zeros([ny, nx, 31, M], float) 49 49 tbs = np.zeros([ny, nx, 31, M], float) 50 tbl = np.zeros([ny, nx, 31, M], float) 50 tbl = np.zeros([ny, nx, 31, M], float)''' 51 51 es = np.zeros([ny, nx, 31, M], float) 52 52 el = np.zeros([ny, nx, 31, M], float) 53 53 esl = np.zeros([ny, nx, 31, M], float) 54 esl00 = np.zeros([ny, nx, 31, M], float)54 '''esl00 = np.zeros([ny, nx, 31, M], float) 55 55 esl25 = np.zeros([ny, nx, 31, M], float) 56 56 esl50 = np.zeros([ny, nx, 31, M], float) 57 57 esl75 = np.zeros([ny, nx, 31, M], float) 58 esl100 = np.zeros([ny, nx, 31, M], float) 58 esl100 = np.zeros([ny, nx, 31, M], float)''' 59 59 for imo in range (0, M): 60 60 print 'month: ' + month[imo] 61 fichier=open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA 23.dat','r')61 fichier=open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA50.dat','r') 62 62 numlines = 0 63 63 for line in fichier: numlines += 1 64 64 fichier.close 65 fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA 23.dat','r')65 fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009_AMSUA50.dat','r') 66 66 nbtotal = numlines-1 67 67 iligne = 0 … … 70 70 lon = np.zeros([nbtotal],float) 71 71 e = np.zeros([nbtotal],float) 72 ts = np.zeros([nbtotal],float)72 '''ts = np.zeros([nbtotal],float) 73 73 tup = np.zeros([nbtotal],float) 74 74 tdn = np.zeros([nbtotal],float) … … 78 78 tdn_lamb = np.zeros([nbtotal],float) 79 79 tb_spec = np.zeros([nbtotal],float) 80 tb_lamb = np.zeros([nbtotal],float) 80 tb_lamb = np.zeros([nbtotal],float)''' 81 81 e_spec = np.zeros([nbtotal],float) 82 82 e_lamb = np.zeros([nbtotal],float) 83 83 e_spec_lamb = np.zeros([nbtotal],float) 84 e_sl_00 = np.zeros([nbtotal],float)84 '''e_sl_00 = np.zeros([nbtotal],float) 85 85 e_sl_25 = np.zeros([nbtotal],float) 86 86 e_sl_50 = np.zeros([nbtotal],float) 87 87 e_sl_75 = np.zeros([nbtotal],float) 88 e_sl_100 = np.zeros([nbtotal],float) 88 e_sl_100 = np.zeros([nbtotal],float)''' 89 89 while (iligne < nbtotal) : 90 90 line=fichier.readline() … … 94 94 lon[iligne] = float(liste[1]) 95 95 e[iligne] = float(liste[5]) 96 ts[iligne] = float(liste[6])96 '''ts[iligne] = float(liste[6]) 97 97 tup[iligne] = float(liste[7]) 98 98 tdn[iligne] = float(liste[8]) … … 102 102 tdn_lamb[iligne] = float(liste[13]) 103 103 tb_spec[iligne] = float(liste[14]) 104 tb_lamb[iligne] = float(liste[15]) 104 tb_lamb[iligne] = float(liste[15])''' 105 105 e_spec[iligne] = float(liste[16]) 106 106 e_lamb[iligne] = float(liste[17]) 107 107 e_spec_lamb[iligne] = float(liste[18]) 108 e_sl_00[iligne] = float(liste[19])108 '''e_sl_00[iligne] = float(liste[19]) 109 109 e_sl_25[iligne] = float(liste[20]) 110 110 e_sl_50[iligne] = float(liste[21]) 111 111 e_sl_75[iligne] = float(liste[22]) 112 e_sl_100[iligne] = float(liste[23]) 112 e_sl_100[iligne] = float(liste[23])''' 113 113 iligne=iligne+1 114 114 fichier.close() 115 print 'tup'115 '''print 'tup' 116 116 z0 = tup.min() 117 117 z1 = tup.max() … … 136 136 zgrid_output, ngrid_output, z2grid_output, sigmagrid_output, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(month_day[imo], jjr, month[imo], lon, lat, tb_lamb, z0, z1, dx, dy) 137 137 tb_lamb_day = zgrid_output 138 tbl[:, :, 0 : month_day[imo], imo] = tb_lamb_day[:, :, :] 138 tbl[:, :, 0 : month_day[imo], imo] = tb_lamb_day[:, :, :]''' 139 139 print 'emis spec' 140 140 z0 = e_spec.min() … … 155 155 e_spec_lamb_day = zgrid_output 156 156 esl[:, :, 0 : month_day[imo], imo] = e_spec_lamb_day[:, :, :] 157 print 'emis spec lamb 00'157 '''print 'emis spec lamb 00' 158 158 z0 = e_sl_00.min() 159 159 z1 = e_sl_00.max() … … 184 184 zgrid_output, ngrid_output, z2grid_output, sigmagrid_output, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(month_day[imo], jjr, month[imo], lon, lat, e_sl_100, z0, z1, dx, dy) 185 185 e_sl_100_day = zgrid_output 186 esl100[:, :, 0 : month_day[imo], imo] = e_sl_100_day[:, :, :] 186 esl100[:, :, 0 : month_day[imo], imo] = e_sl_100_day[:, :, :]''' 187 187 ############################################### 188 188 # stack gridded data spec lamb in NETCDF file # … … 190 190 print 'stacking of gridded data' 191 191 print 'file for ' + month[imo] 192 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA 23_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC')192 rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA50_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC') 193 193 rootgrp.createDimension('longitude', nx) 194 194 rootgrp.createDimension('latitude', ny) … … 197 197 nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) 198 198 nc_day = rootgrp.createVariable('days', 'f', ('days',)) 199 nc_tup = rootgrp.createVariable('tup', 'f', ('latitude', 'longitude', 'days'))199 '''nc_tup = rootgrp.createVariable('tup', 'f', ('latitude', 'longitude', 'days')) 200 200 nc_tdn = rootgrp.createVariable('tdn', 'f', ('latitude', 'longitude', 'days')) 201 201 nc_tb_spec = rootgrp.createVariable('tb_spec', 'f', ('latitude', 'longitude', 'days')) 202 nc_tb_lamb = rootgrp.createVariable('tb_lamb', 'f', ('latitude', 'longitude', 'days')) 202 nc_tb_lamb = rootgrp.createVariable('tb_lamb', 'f', ('latitude', 'longitude', 'days'))''' 203 203 nc_e_spec = rootgrp.createVariable('e_spec', 'f', ('latitude', 'longitude', 'days')) 204 204 nc_e_lamb = rootgrp.createVariable('e_lamb', 'f', ('latitude', 'longitude', 'days')) 205 205 nc_e_spec_lamb = rootgrp.createVariable('e_spec_lamb', 'f', ('latitude', 'longitude', 'days')) 206 nc_e_sl_00 = rootgrp.createVariable('e_mixed_s00', 'f', ('latitude', 'longitude', 'days'))206 '''nc_e_sl_00 = rootgrp.createVariable('e_mixed_s00', 'f', ('latitude', 'longitude', 'days')) 207 207 nc_e_sl_25 = rootgrp.createVariable('e_mixed_s25', 'f', ('latitude', 'longitude', 'days')) 208 208 nc_e_sl_50 = rootgrp.createVariable('e_mixed_s50', 'f', ('latitude', 'longitude', 'days')) 209 209 nc_e_sl_75 = rootgrp.createVariable('e_mixed_s75', 'f', ('latitude', 'longitude', 'days')) 210 nc_e_sl_100 = rootgrp.createVariable('e_mixed_s100', 'f', ('latitude', 'longitude', 'days')) 210 nc_e_sl_100 = rootgrp.createVariable('e_mixed_s100', 'f', ('latitude', 'longitude', 'days'))''' 211 211 nc_lon[:] = xvec 212 212 nc_lat[:] = yvec 213 nc_tup[:] = tu[:, :, 0 : month_day[imo], imo]213 '''nc_tup[:] = tu[:, :, 0 : month_day[imo], imo] 214 214 nc_tdn[:] = td[:, :, 0 : month_day[imo], imo] 215 215 nc_tb_spec[:] = tbs[:, :, 0 : month_day[imo], imo] 216 nc_tb_lamb[:] = tbl[:, :, 0 : month_day[imo], imo] 216 nc_tb_lamb[:] = tbl[:, :, 0 : month_day[imo], imo]''' 217 217 nc_e_spec[:] = es[:, :, 0 : month_day[imo], imo] 218 218 nc_e_lamb[:] = el[:, :, 0 : month_day[imo], imo] 219 219 nc_e_spec_lamb[:] = esl[:, :, 0 : month_day[imo], imo] 220 nc_e_sl_00[:] = esl00[:, :, 0 : month_day[imo], imo]220 '''nc_e_sl_00[:] = esl00[:, :, 0 : month_day[imo], imo] 221 221 nc_e_sl_25[:] = esl25[:, :, 0 : month_day[imo], imo] 222 222 nc_e_sl_50[:] = esl50[:, :, 0 : month_day[imo], imo] 223 223 nc_e_sl_75[:] = esl75[:, :, 0 : month_day[imo], imo] 224 nc_e_sl_100[:] = esl100[:, :, 0 : month_day[imo], imo] 224 nc_e_sl_100[:] = esl100[:, :, 0 : month_day[imo], imo]''' 225 225 rootgrp.close() 226 226 -
trunk/src/scripts_Laura/ARCTIC/Travail_CEN/temporal_evolution_AMSUA23-AMSUB89.py
r43 r44 100 100 101 101 102 #90, 51 103 #95, 82 104 105 spec_amsua = np.append(es_a[0, 95, 82, 0 : 31], es_a[1, 95, 82, 0 : 28]) 106 spec_amsub = np.append(es_b[0, 95, 82, 0 : 31], es_b[1, 95, 82, 0 : 28]) 107 ambig = np.append(amb_ice[0, 95, 82, 0 : 31], amb_ice[1, 95, 82, 0 : 28]) 108 for imo in range (2, M): 109 spec_amsua = np.append(spec_amsua, es_a[imo, 95, 82, 0 : month_day[imo]]) 110 spec_amsub = np.append(spec_amsub, es_b[imo, 95, 82, 0 : month_day[imo]]) 111 ambig = np.append(ambig, amb_ice[imo, 95, 82, 0 : month_day[imo]]) 102 112 103 113 104 spec_amsua = np.append(es_a[0, 119, 51, 0 : 31], es_a[1, 119, 51, 0 : 31]) 105 spec_amsub = np.append(es_b[0, 119, 51, 0 : 31], es_b[1, 119, 51, 0 : 31]) 106 ambig = np.append(amb_ice[0, 119, 51, 0 : 31], amb_ice[1, 119, 51, 0 : 31]) 107 for imo in range (2, M): 108 spec_amsua = np.append(spec_amsua, es_a[imo, 119, 51, 0 : month_day[imo]]) 109 spec_amsub = np.append(spec_amsub, es_b[imo, 119, 51, 0 : month_day[imo]]) 110 ambig = np.append(ambig, amb_ice[imo, 119, 51, 0 : month_day[imo]]) 111 114 # seperately AMSUA and AMSUB 112 115 ion() 113 116 figure() 114 plot( spec_amsua[nonzero(spec_amsua != 0.)], 'b', label = 'AMSUA 23GHz')115 plot( spec_amsub[nonzero(spec_amsub != 0.)], 'g', label = 'AMSUB 89GHz')116 plot(np.arange(0, 36 8, 1)[nonzero(ambig == 1.)], spec_amsua[nonzero(ambig == 1.)], 'ob', label = 'obs inambiguous area')117 plot(np.arange(0, 36 8, 1)[nonzero(ambig == 1.)], spec_amsub[nonzero(ambig == 1.)], 'og', label = 'obs inambiguous area')117 plot(np.arange(0, 365, 1), spec_amsua, 'b', label = 'AMSUA 23GHz') 118 plot(np.arange(0, 365, 1), spec_amsub, 'g', label = 'AMSUB 89GHz') 119 plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_amsua[nonzero(ambig == 1.)], 'ob', label = 'ambiguous area') 120 plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_amsub[nonzero(ambig == 1.)], 'og', label = 'ambiguous area') 118 121 legend(loc = 3) 122 xlim(-3, 368) 123 xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 334]), month, rotation = 20) 124 grid() 125 126 # bias AMSUA-AMSUB 127 figure() 128 plot(np.arange(0, 365, 1), spec_amsua - spec_amsub, 'r', label = 'AMSUA23 - AMSUB89') 129 plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], (spec_amsua - spec_amsub)[nonzero(ambig == 1.)], 'or', label = 'ambiguous area') 130 plot(np.arange(0, 365, 1), np.zeros([365], float), 'k') 131 legend(loc = 3) 132 xlim(-3, 368) 133 xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 335]), month, rotation = 20) 134 135 136 # anomaly with a 15 day time range 137 spec_anom = np.zeros([365], float) 138 spec_diff = spec_amsua - spec_amsub 139 for ii in range (15, 365 - 15): 140 spec_anom[ii] = (spec_diff[ii]) - mean(spec_diff[ii - 15 : ii + 15][nonzero(isnan(spec_diff[ii - 15 : ii + 15]) == False)]) 141 142 figure() 143 plot(np.arange(0, 365, 1), spec_anom, 'r') 144 plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_anom[nonzero(ambig == 1.)], 'or', label = 'ambiguous area') 145 plot(np.arange(0, 365, 1), np.zeros([365], float), 'k') 146 legend(loc = 3) 147 xlim(-3, 368) 148 xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 335]), month, rotation = 20) 149 150 151 152 119 153 120 154 # plot of anomaly (espec amsub - espec amsua) on a 15 day sliding base (15 day correct to simulate season transition)
Note: See TracChangeset
for help on using the changeset viewer.