[51] | 1 | #!/usr/bin/env python |
---|
| 2 | # -*- coding: utf-8 -*- |
---|
| 3 | import string |
---|
| 4 | import numpy as np |
---|
| 5 | import matplotlib.pyplot as plt |
---|
| 6 | from pylab import * |
---|
| 7 | from mpl_toolkits.basemap import Basemap |
---|
| 8 | from mpl_toolkits.basemap import shiftgrid, cm |
---|
| 9 | from netCDF4 import Dataset |
---|
| 10 | import arctic_map # function to regrid coast limits |
---|
| 11 | import cartesian_grid_test # function to convert grid from polar to cartesian |
---|
| 12 | import scipy.special |
---|
| 13 | import ffgrid2 |
---|
| 14 | import map_ffgrid |
---|
| 15 | from matplotlib import colors |
---|
| 16 | from matplotlib.font_manager import FontProperties |
---|
| 17 | import map_cartesian_grid |
---|
| 18 | |
---|
| 19 | |
---|
| 20 | ############################### |
---|
| 21 | # time period characteristics # |
---|
| 22 | ############################### |
---|
| 23 | MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) |
---|
| 24 | month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER']) |
---|
| 25 | month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) |
---|
| 26 | M = len(month) |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | ######################## |
---|
| 30 | # grid characteristics # |
---|
| 31 | ######################## |
---|
| 32 | x0 = -3000. # min limit of grid |
---|
| 33 | x1 = 2500. # max limit of grid |
---|
| 34 | dx = 40. |
---|
| 35 | xvec = np.arange(x0, x1+dx, dx) |
---|
| 36 | nx = len(xvec) |
---|
| 37 | y0 = -3000. # min limit of grid |
---|
| 38 | y1 = 3000. # max limit of grid |
---|
| 39 | dy = 40. |
---|
| 40 | yvec = np.arange(y0, y1+dy, dy) |
---|
| 41 | ny = len(yvec) |
---|
| 42 | |
---|
| 43 | |
---|
| 44 | |
---|
[56] | 45 | |
---|
[51] | 46 | ############################################ |
---|
| 47 | # time evolution (monthly) in a given zone # |
---|
| 48 | ############################################ |
---|
| 49 | #zone 1 (seasonal ice) : yi = 960. // yf = 1360. // xi = -680. // xf = -320. (North Beaufort Sea) |
---|
| 50 | #zone 2 (multiyear ice) : yi = 320. // yf = 720. // xi = -1080. // xf = -720. (North Canadian Archipelago) |
---|
| 51 | #zone 3 (young ice) : yi = 1880. // yf = 2280. // xi = -480. // xf = -120. (Chukchi Sea) |
---|
| 52 | |
---|
| 53 | # select borders of zone |
---|
| 54 | yi = 320. |
---|
| 55 | yf = 720. |
---|
| 56 | xi = -1080. |
---|
| 57 | xf = -720. |
---|
| 58 | |
---|
| 59 | #find corresponding index in xvec and yvec |
---|
| 60 | xxi = np.where(xvec == xi)[0][0] |
---|
| 61 | xxf = np.where(xvec == xf)[0][0] |
---|
| 62 | yyi = np.where(yvec == yi)[0][0] |
---|
| 63 | yyf = np.where(yvec == yf)[0][0] |
---|
| 64 | |
---|
| 65 | len(xvec[xxi:xxf+1]) |
---|
| 66 | len(yvec[yyi:yyf+1]) |
---|
| 67 | |
---|
| 68 | |
---|
[56] | 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 |
---|
[51] | 71 | |
---|
[56] | 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) |
---|
[51] | 74 | |
---|
[56] | 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) |
---|
[51] | 77 | |
---|
[56] | 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) |
---|
[51] | 80 | |
---|
[56] | 81 | S = np.zeros([M, 31], float) # number of data pixels in selected zone, per day and per month |
---|
| 82 | |
---|
[51] | 83 | for imo in range (0, M): |
---|
| 84 | # daily read gradient ratio file |
---|
| 85 | print 'read daily gradient ratio for month ' + month[imo] |
---|
| 86 | fichier_grad_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_grad_ratio_spec_23-89_near_nadir_AMSUA_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') |
---|
| 87 | gr = fichier_grad_ratio.variables['grad_ratio'][:] |
---|
| 88 | fichier_grad_ratio.close() |
---|
| 89 | # read daily emis anomaly file for 23GHz |
---|
| 90 | print 'read daily emis anomaly 23GHz for month ' + month[imo] |
---|
| 91 | fichier_anom23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') |
---|
| 92 | sa_23 = fichier_anom23.variables['spec_anomaly'][:] |
---|
| 93 | fichier_anom23.close() |
---|
| 94 | # read daily emis anomaly file for 89GHz |
---|
| 95 | print 'read daily emis anomaly 89GHz for month ' + month[imo] |
---|
| 96 | fichier_anom89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') |
---|
| 97 | sa_89 = fichier_anom89.variables['spec_anomaly'][:] |
---|
| 98 | ra_89 = fichier_anom89.variables['ratio_anomaly'][:] |
---|
| 99 | fichier_anom89.close() |
---|
| 100 | grad_ratio_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float) |
---|
| 101 | spec_anom_23_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float) |
---|
| 102 | spec_anom_89_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float) |
---|
| 103 | ratio_anom_89_vec = np.zeros([month_day[imo], len(xvec[xxi : xxf+1]) * len(yvec[yyi : yyf+1])], float) |
---|
| 104 | print 'calculate daily mean and std zonal gradient ratio' |
---|
| 105 | for ijr in range (0, month_day[imo]): |
---|
| 106 | print 'date ' + str(ijr+1) + ' ' + month[imo] + ' 2009' |
---|
| 107 | S[imo, ijr] = shape(gr[ijr, yyi : yyf + 1, xxi : xxf + 1][nonzero(isnan(gr[ijr, yyi : yyf + 1, xxi : xxf + 1]) == False)])[0] |
---|
| 108 | # gradient ratio |
---|
| 109 | grad_ratio_vec[ijr, :] = reshape(gr[ijr, yyi : yyf + 1, xxi : xxf + 1], size(gr[ijr, yyi : yyf + 1, xxi : xxf + 1])) |
---|
| 110 | mean_grad_ratio_zone[imo, ijr] = mean(grad_ratio_vec[ijr, :][nonzero(isnan(grad_ratio_vec[ijr, :]) == False)]) |
---|
| 111 | std_grad_ratio_zone[imo, ijr] = sqrt((1. / (size(gr[ijr, yyi : yyf+1, xxi : xxf+1]) - 1.)) * sum((grad_ratio_vec[ijr, :][nonzero(isnan(grad_ratio_vec[ijr, :]) == False)] - mean_grad_ratio_zone[imo, ijr])**2)) |
---|
| 112 | # spec anomaly 23GHz |
---|
| 113 | spec_anom_23_vec[ijr, :] = reshape(sa_23[ijr, yyi : yyf + 1, xxi : xxf + 1], size(sa_23[ijr, yyi : yyf + 1, xxi : xxf + 1])) |
---|
| 114 | mean_spec_anom_23_zone[imo, ijr] = mean(spec_anom_23_vec[ijr, :][nonzero(isnan(spec_anom_23_vec[ijr, :]) == False)]) |
---|
| 115 | std_spec_anom_23_zone[imo, ijr] = sqrt((1. / (size(sa_23[ijr, yyi : yyf+1, xxi : xxf+1]) - 1.)) * sum((spec_anom_23_vec[ijr, :][nonzero(isnan(spec_anom_23_vec[ijr, :]) == False)] - mean_spec_anom_23_zone[imo, ijr])**2)) |
---|
| 116 | # spec anomaly 89GHz |
---|
| 117 | spec_anom_89_vec[ijr, :] = reshape(sa_89[ijr, yyi : yyf + 1, xxi : xxf + 1], size(sa_89[ijr, yyi : yyf + 1, xxi : xxf + 1])) |
---|
| 118 | mean_spec_anom_89_zone[imo, ijr] = mean(spec_anom_89_vec[ijr, :][nonzero(isnan(spec_anom_89_vec[ijr, :]) == False)]) |
---|
| 119 | std_spec_anom_89_zone[imo, ijr] = sqrt((1. / (size(sa_89[ijr, yyi : yyf+1, xxi : xxf+1]) - 1.)) * sum((spec_anom_89_vec[ijr, :][nonzero(isnan(spec_anom_89_vec[ijr, :]) == False)] - mean_spec_anom_89_zone[imo, ijr])**2)) |
---|
| 120 | # ratio anomaly 89GHz |
---|
| 121 | ratio_anom_89_vec[ijr, :] = reshape(ra_89[ijr, yyi : yyf + 1, xxi : xxf + 1], size(ra_89[ijr, yyi : yyf + 1, xxi : xxf + 1])) |
---|
| 122 | mean_ratio_anom_89_zone[imo, ijr] = mean(ratio_anom_89_vec[ijr, :][nonzero(isnan(ratio_anom_89_vec[ijr, :]) == False)]) |
---|
| 123 | std_ratio_anom_89_zone[imo, ijr] = sqrt((1. / (size(ra_89[ijr, yyi : yyf+1, xxi : xxf+1]) - 1.)) * sum((ratio_anom_89_vec[ijr, :][nonzero(isnan(ratio_anom_89_vec[ijr, :]) == False)] - mean_ratio_anom_89_zone[imo, ijr])**2)) |
---|
| 124 | if (isnan(mean_grad_ratio_zone[imo, ijr]) == True): |
---|
| 125 | std_grad_ratio_zone[imo, ijr] = NaN |
---|
| 126 | if (isnan(mean_spec_anom_23_zone[imo, ijr]) == True): |
---|
| 127 | std_spec_anom_23_zone[imo, ijr] = NaN |
---|
| 128 | if (isnan(mean_spec_anom_89_zone[imo, ijr]) == True): |
---|
| 129 | std_spec_anom_89_zone[imo, ijr] = NaN |
---|
| 130 | if (isnan(mean_ratio_anom_89_zone[imo, ijr]) == True): |
---|
| 131 | std_ratio_anom_89_zone[imo, ijr] = NaN |
---|
| 132 | |
---|
| 133 | |
---|
| 134 | |
---|
| 135 | # append daily zonal gradient ratio for study of the whole year 2009 |
---|
| 136 | mean_year_grad_ratio = np.append(mean_grad_ratio_zone[0, 0 : month_day[0]], mean_grad_ratio_zone[1, 0 : month_day[1]]) |
---|
| 137 | std_year_grad_ratio = np.append(std_grad_ratio_zone[0, 0 : month_day[0]], std_grad_ratio_zone[1, 0 : month_day[1]]) |
---|
| 138 | |
---|
| 139 | mean_year_spec_anom_23 = np.append(mean_spec_anom_23_zone[0, 0 : month_day[0]], mean_spec_anom_23_zone[1, 0 : month_day[1]]) |
---|
| 140 | std_year_spec_anom_23 = np.append(std_spec_anom_23_zone[0, 0 : month_day[0]], std_spec_anom_23_zone[1, 0 : month_day[1]]) |
---|
| 141 | |
---|
| 142 | mean_year_spec_anom_89 = np.append(mean_spec_anom_89_zone[0, 0 : month_day[0]], mean_spec_anom_89_zone[1, 0 : month_day[1]]) |
---|
| 143 | std_year_spec_anom_89 = np.append(std_spec_anom_89_zone[0, 0 : month_day[0]], std_spec_anom_89_zone[1, 0 : month_day[1]]) |
---|
| 144 | |
---|
| 145 | mean_year_ratio_anom_89 = np.append(mean_ratio_anom_89_zone[0, 0 : month_day[0]], mean_ratio_anom_89_zone[1, 0 : month_day[1]]) |
---|
| 146 | std_year_ratio_anom_89 = np.append(std_ratio_anom_89_zone[0, 0 : month_day[0]], std_ratio_anom_89_zone[1, 0 : month_day[1]]) |
---|
| 147 | |
---|
| 148 | year_S = np.append(S[0, 0 : month_day[0]], S[1, 0 : month_day[1]]) |
---|
| 149 | |
---|
| 150 | for imo in range (2, M): |
---|
| 151 | # gradient ratio |
---|
| 152 | mean_year_grad_ratio = np.append(mean_year_grad_ratio, mean_grad_ratio_zone[imo, 0 : month_day[imo]]) |
---|
| 153 | std_year_grad_ratio = np.append(std_year_grad_ratio, std_grad_ratio_zone[imo, 0 : month_day[imo]]) |
---|
| 154 | # spec anomaly 23GHz |
---|
| 155 | mean_year_spec_anom_23 = np.append(mean_year_spec_anom_23, mean_spec_anom_23_zone[imo, 0 : month_day[imo]]) |
---|
| 156 | std_year_spec_anom_23 = np.append(std_year_spec_anom_23, std_spec_anom_23_zone[imo, 0 : month_day[imo]]) |
---|
| 157 | # spec anomaly 89GHz |
---|
| 158 | mean_year_spec_anom_89 = np.append(mean_year_spec_anom_89, mean_spec_anom_89_zone[imo, 0 : month_day[imo]]) |
---|
| 159 | std_year_spec_anom_89 = np.append(std_year_spec_anom_89, std_spec_anom_89_zone[imo, 0 : month_day[imo]]) |
---|
| 160 | # ratio anomaly 89GHz |
---|
| 161 | mean_year_ratio_anom_89 = np.append(mean_year_ratio_anom_89, mean_ratio_anom_89_zone[imo, 0 : month_day[imo]]) |
---|
| 162 | std_year_ratio_anom_89 = np.append(std_year_ratio_anom_89, std_ratio_anom_89_zone[imo, 0 : month_day[imo]]) |
---|
| 163 | # number of data points in area of study |
---|
| 164 | year_S = np.append(year_S, S[imo, 0 : month_day[imo]]) |
---|
| 165 | |
---|
| 166 | |
---|
| 167 | |
---|
[56] | 168 | ###################################################################################### |
---|
| 169 | # calculate standard deviation of emissivity parameters over the year 2009 (364 days)# |
---|
| 170 | ###################################################################################### |
---|
[51] | 171 | print 'calculate standard deviation of emissivity parameters over the year 2009' |
---|
| 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)) |
---|
| 173 | |
---|
| 174 | time_std2 = sqrt((1./364.)*sum((mean_year_spec_anom_23[nonzero(isnan(mean_year_spec_anom_23) == False)] - mean(mean_year_spec_anom_23[nonzero(isnan(mean_year_spec_anom_23) == False)]))**2)) |
---|
| 175 | |
---|
| 176 | time_std3 = sqrt((1./364.)*sum((mean_year_spec_anom_89[nonzero(isnan(mean_year_spec_anom_89) == False)] - mean(mean_year_spec_anom_89[nonzero(isnan(mean_year_spec_anom_89) == False)]))**2)) |
---|
| 177 | |
---|
| 178 | time_std4 = sqrt((1./364.)*sum((mean_year_ratio_anom_89[nonzero(isnan(mean_year_ratio_anom_89) == False)] - mean(mean_year_ratio_anom_89[nonzero(isnan(mean_year_ratio_anom_89) == False)]))**2)) |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | |
---|
| 183 | ######################### |
---|
| 184 | # plot daily parameters # |
---|
| 185 | ######################### |
---|
| 186 | vec_months = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) |
---|
| 187 | ion() |
---|
| 188 | figure() |
---|
| 189 | subplot(2, 1, 1) |
---|
| 190 | plot(mean_year_grad_ratio, 'b', label = 'grad ratio') |
---|
| 191 | plot(mean_year_spec_anom_23, 'r', label = 'spec anom 23GHz') |
---|
| 192 | plot(mean_year_spec_anom_89, 'g', label = 'spec anom 89GHz') |
---|
| 193 | plot(np.zeros([365], float), '--k') |
---|
| 194 | xticks(vec_months, month, rotation = 25) |
---|
| 195 | xlim(0, 365) |
---|
| 196 | fontP = FontProperties() |
---|
| 197 | fontP.set_size('small') |
---|
| 198 | legend(loc = 3, prop = fontP) |
---|
| 199 | grid() |
---|
| 200 | subplot(2, 1, 2) |
---|
| 201 | plot(mean_year_ratio_anom_89, 'y', label = 'ratio anom 89GHz') |
---|
| 202 | plot(np.zeros([365], float), '--k') |
---|
| 203 | xlim(0, 365) |
---|
| 204 | xticks(vec_months, month, rotation = 25) |
---|
| 205 | legend(loc = 1, prop = fontP) |
---|
| 206 | grid() |
---|
| 207 | 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_grad_ratio_emis_anomaly_params_zone_Chukchi_Sea.png') |
---|
| 208 | ############################################ |
---|
| 209 | # plot daily number of data points in area # |
---|
| 210 | ############################################ |
---|
| 211 | vec_months = np.array([0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]) |
---|
| 212 | figure() |
---|
| 213 | plot(year_S, label = 'mean=' + str(mean(year_S))[0:5] + ' ; std=' + str(std(year_S))[0:5]) |
---|
| 214 | xticks(vec_months, month, rotation = 25) |
---|
| 215 | xlim(0, 365) |
---|
| 216 | ylabel('Number of data points in area') |
---|
| 217 | fontP = FontProperties() |
---|
| 218 | fontP.set_size('small') |
---|
| 219 | legend(loc = 1, prop = fontP) |
---|
| 220 | title('North Canadian Archipelago') |
---|
| 221 | 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_nber_data_points_zone_North_Canadian_Archipelago.png') |
---|
| 222 | |
---|
| 223 | |
---|
[52] | 224 | |
---|
[51] | 225 | subplot(2, 1, 2) |
---|
| 226 | plot(std_year_grad_ratio, '--b', label = 'std grad ratio') |
---|
| 227 | plot(std_year_grad_ratio, '--r', label = 'std grad ratio') |
---|
| 228 | plot(std_year_grad_ratio, '--g', label = 'std grad ratio') |
---|
| 229 | plot(std_year_grad_ratio, '--y', label = 'std grad ratio') |
---|
| 230 | xticks(vec_months, month, rotation = 25) |
---|
| 231 | legend(loc = 2, prop = fontP) |
---|
| 232 | xlim(0, 365) |
---|
| 233 | |
---|
| 234 | |
---|
[52] | 235 | |
---|
[51] | 236 | #################### |
---|
| 237 | # map studied zone # |
---|
| 238 | #################### |
---|
| 239 | ion() |
---|
| 240 | x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() |
---|
| 241 | x_coast = x_ind |
---|
| 242 | y_coast = y_ind |
---|
| 243 | z_coast = z_ind |
---|
| 244 | map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec[xxi : xxf + 1], yvec[yyi : yyf + 1], gr[0, yyi : yyf + 1, xxi : xxf + 1], gr[0, :, :][nonzero(isnan(gr[0, :, :]) == False)].min(), gr[0, :, :][nonzero(isnan(gr[0, :, :]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'daily gradient ratio (01-01-2009)') |
---|
| 245 | title('area of study - Chukchi Sea') |
---|
| 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') |
---|
| 247 | |
---|
| 248 | |
---|
[56] | 249 | |
---|
| 250 | ############################################ |
---|
| 251 | # read emissivity parameters in all Arctic # |
---|
| 252 | ############################################ |
---|
[52] | 253 | cumul_params = np.zeros([M, ny, nx], float) |
---|
| 254 | grad_ratio = np.zeros([M, ny, nx], float) |
---|
| 255 | spec_anom_23 = np.zeros([M, ny, nx], float) |
---|
| 256 | spec_anom_89 = np.zeros([M, ny, nx], float) |
---|
| 257 | ratio_anom_89 = np.zeros([M, ny, nx], float) |
---|
| 258 | CP = np.zeros([M, ny, nx], float) |
---|
| 259 | for imo in range (0, M): |
---|
| 260 | # daily read gradient ratio file |
---|
| 261 | print 'read daily gradient ratio for month ' + month[imo] |
---|
| 262 | fichier_grad_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_grad_ratio_spec_23-89_near_nadir_AMSUA_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') |
---|
| 263 | gr = fichier_grad_ratio.variables['grad_ratio'][:] |
---|
| 264 | fichier_grad_ratio.close() |
---|
| 265 | # read daily emis anomaly file for 23GHz |
---|
| 266 | print 'read daily emis anomaly 23GHz for month ' + month[imo] |
---|
| 267 | fichier_anom23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') |
---|
| 268 | sa_23 = fichier_anom23.variables['spec_anomaly'][:] |
---|
| 269 | fichier_anom23.close() |
---|
| 270 | # read daily emis anomaly file for 89GHz |
---|
| 271 | print 'read daily emis anomaly 89GHz for month ' + month[imo] |
---|
| 272 | fichier_anom89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC') |
---|
| 273 | sa_89 = fichier_anom89.variables['spec_anomaly'][:] |
---|
| 274 | ra_89 = fichier_anom89.variables['ratio_anomaly'][:] |
---|
| 275 | fichier_anom89.close() |
---|
| 276 | for ilon in range (0, nx): |
---|
| 277 | for ilat in range (0, ny): |
---|
[56] | 278 | # calculate monthly mean of emissivity parameters |
---|
[52] | 279 | grad_ratio[imo, ilat, ilon] = mean(gr[:, ilat, ilon][nonzero(isnan(gr[:, ilat, ilon]) == False)]) |
---|
| 280 | spec_anom_23[imo, ilat, ilon] = mean(sa_23[:, ilat, ilon][nonzero(isnan(sa_23[:, ilat, ilon]) == False)]) |
---|
| 281 | spec_anom_89[imo, ilat, ilon] = mean(sa_89[:, ilat, ilon][nonzero(isnan(sa_89[:, ilat, ilon]) == False)]) |
---|
| 282 | ratio_anom_89[imo, ilat, ilon] = mean(ra_89[:, ilat, ilon][nonzero(isnan(ra_89[:, ilat, ilon]) == False)]) |
---|
[56] | 283 | # calculate monthly cumulation index = (sum(abs(all emissivity parameters))/maximum of this sum)*100 (to have a percentage) |
---|
[52] | 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]) |
---|
| 285 | CP[imo, :, :] = (cumul_params[imo, :, :]/max(cumul_params[imo, :, :][nonzero(isnan(cumul_params[imo, :, :]) == False)])) * 100. |
---|
[51] | 286 | |
---|
[52] | 287 | |
---|
[56] | 288 | ####################### |
---|
| 289 | # map cuulation index # |
---|
| 290 | ####################### |
---|
[52] | 291 | #ion() |
---|
| 292 | x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() |
---|
| 293 | x_coast = x_ind |
---|
| 294 | y_coast = y_ind |
---|
| 295 | z_coast = z_ind |
---|
| 296 | for imo in range (0, M): |
---|
| 297 | print 'map month ' + month[imo] |
---|
| 298 | map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, CP[imo, :, :], 0., 100., 1., cm.s3pcpn_l_r, 'Cumulation index (%)') |
---|
| 299 | title('AMSUA - ' + month[imo] + ' 2009') |
---|
| 300 | savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/cumul_params/cumul_all_parameters/map_cumulation_index_'+ str(MONTH[imo]) + month[imo] + '2009.png') |
---|