| 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 | |
|---|
| 13 | |
|---|
| 14 | le_mois=np.array(['NOVEMBER2009', 'DECEMBER2009', 'JANUARY2010', 'FEBRUARY2010']) |
|---|
| 15 | month_days = np.array([30, 31, 31, 28]) |
|---|
| 16 | |
|---|
| 17 | |
|---|
| 18 | |
|---|
| 19 | for imo in range (0, 1): |
|---|
| 20 | fichier=open("/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_EMIS_GLACE_"+str(le_mois[imo])+".DAT",'r') |
|---|
| 21 | numlines = 0 |
|---|
| 22 | for line in fichier: numlines += 1 |
|---|
| 23 | fichier.close() |
|---|
| 24 | fichier=open("/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_EMIS_GLACE_"+str(le_mois[imo])+".DAT",'r') |
|---|
| 25 | nbtotal=numlines-1 |
|---|
| 26 | iligne=0 |
|---|
| 27 | lat=np.zeros([nbtotal],float) |
|---|
| 28 | lon=np.zeros([nbtotal],float) |
|---|
| 29 | jjr=np.zeros([nbtotal],float) |
|---|
| 30 | zen=np.zeros([nbtotal],float) |
|---|
| 31 | fov=np.zeros([nbtotal],float) |
|---|
| 32 | ts=np.zeros([nbtotal],float) |
|---|
| 33 | emis=np.zeros([4,nbtotal],float) |
|---|
| 34 | tb=np.zeros([4,nbtotal],float) |
|---|
| 35 | tup=np.zeros([4,nbtotal],float) |
|---|
| 36 | tdn=np.zeros([4,nbtotal],float) |
|---|
| 37 | #trans=np.zeros([4,nbtotal],float) |
|---|
| 38 | while (iligne < nbtotal) : |
|---|
| 39 | line = fichier.readline() |
|---|
| 40 | # exemple : line = "0.22 2.3 5.0 6" |
|---|
| 41 | liste = line.split() |
|---|
| 42 | # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es) |
|---|
| 43 | lat[iligne] = float(liste[1]) |
|---|
| 44 | lon[iligne] = float(liste[0]) |
|---|
| 45 | jjr[iligne] = float(liste[4]) |
|---|
| 46 | ts[iligne] = float(liste[8]) |
|---|
| 47 | fov[iligne] = float(liste[10]) |
|---|
| 48 | #trans[iligne] = float(liste[17]) |
|---|
| 49 | for ki in range(0,4): |
|---|
| 50 | emis[ki,iligne] = float(liste[12+ki]) |
|---|
| 51 | tb[ki,iligne] = float(liste[16+ki]) |
|---|
| 52 | tdn[ki,iligne] = float(liste[20+ki]) |
|---|
| 53 | tup[ki,iligne] = float(liste[24+ki]) |
|---|
| 54 | #trans[ki,iligne] = float(liste[28+ki]) |
|---|
| 55 | #tup_ssmis[ki,iligne] = float(liste[28+ki]) |
|---|
| 56 | #tdn_ssmis[ki,iligne] = float(liste[37+ki]) |
|---|
| 57 | #trans_ssmis[ki,iligne] = float(liste[44+ki]) |
|---|
| 58 | iligne=iligne+1 |
|---|
| 59 | fichier.close() |
|---|
| 60 | jour = jjr |
|---|
| 61 | month = month_days[imo] #number of days in month |
|---|
| 62 | longi = lon |
|---|
| 63 | lati = lat |
|---|
| 64 | z_ini = tb[0,:] |
|---|
| 65 | z0 = min(z_ini) |
|---|
| 66 | z1 = max(z_ini) |
|---|
| 67 | dx = 50. |
|---|
| 68 | dy = 50. |
|---|
| 69 | zgrid_output, ngrid_output, z2grid_output, sigmagrid_output, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(jour, month, longi, lati, z_ini, z0, z1, dx, dy) |
|---|
| 70 | tbch1_nov2009 = zgrid_output |
|---|
| 71 | |
|---|
| 72 | |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | ## study zone 1 (ice after summer melt) ## |
|---|
| 76 | # x around -587.198 // y around 1187.5 |
|---|
| 77 | index_x = np.where(xvec == -550.)[0][0] |
|---|
| 78 | index_y = np.where(yvec == 1150.)[0][0] |
|---|
| 79 | a = tbch1_nov2009[index_y, index_x, :] |
|---|
| 80 | b = ngrid_output[index_y, index_x, :] |
|---|
| 81 | c = sigmagrid_output[index_y, index_x, :] |
|---|
| 82 | |
|---|
| 83 | |
|---|
| 84 | plt.ion() |
|---|
| 85 | plt.figure() |
|---|
| 86 | plt.plot(a, 'r', label = 'Tb ch1') |
|---|
| 87 | legend() |
|---|
| 88 | plt.figure() |
|---|
| 89 | plt.plot(b, 'g', label = 'n values in cell') |
|---|
| 90 | legend() |
|---|
| 91 | plt.figure() |
|---|
| 92 | plt.plot(c, 'b', label = 'variance Tb values in cell' ) |
|---|
| 93 | legend() |
|---|
| 94 | |
|---|
| 95 | ## interessant de comparer les Tb avec une resolution plus ou moins fine ## |
|---|
| 96 | |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | tbmean = np.zeros([len(yvec), len(xvec)], float) |
|---|
| 100 | nmean = np.zeros([len(yvec), len(xvec)], float) |
|---|
| 101 | sigmamean = np.zeros([len(yvec), len(xvec)], float) |
|---|
| 102 | for ix in range (0, len(xvec)): |
|---|
| 103 | for iy in range (0, len(yvec)): |
|---|
| 104 | tbmean[iy, ix] = tbch1_nov2009[iy, ix, :][nonzero(isnan(tbch1_nov2009[iy, ix, :]) == False)].mean() |
|---|
| 105 | nmean[iy, ix] = ngrid_output[iy, ix, :][nonzero(isnan(ngrid_output[iy, ix, :]) == False)].mean() |
|---|
| 106 | sigmamean[iy, ix] = sigmagrid_output[iy, ix, :][nonzero(isnan(sigmagrid_output[iy, ix, :]) == False)].mean() |
|---|
| 107 | |
|---|
| 108 | |
|---|
| 109 | |
|---|
| 110 | |
|---|
| 111 | x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() |
|---|
| 112 | plt.figure() |
|---|
| 113 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 114 | clevs = np.arange(137., 269., 1.) |
|---|
| 115 | cs = plt.contourf(xgrid_cart,ygrid_cart, tbmean, clevs, cmap=cm.s3pcpn_l_r) |
|---|
| 116 | cbar = colorbar(cs) |
|---|
| 117 | xlim(-2500., 3000.) |
|---|
| 118 | ylim(-3000., 3000.) |
|---|
| 119 | #txt = 'variance - spatial resolution = ' + str(delta_xy[idel]) + 'km' |
|---|
| 120 | #cbar.set_label(txt) |
|---|
| 121 | plt.figure() |
|---|
| 122 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 123 | clevs = np.arange(0, 30, 1) |
|---|
| 124 | cs = plt.contourf(xgrid_cart,ygrid_cart, nmean, clevs, cmap=cm.s3pcpn_l_r) |
|---|
| 125 | cbar = colorbar(cs) |
|---|
| 126 | xlim(-2500., 3000.) |
|---|
| 127 | ylim(-3000., 3000.) |
|---|
| 128 | plt.figure() |
|---|
| 129 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 130 | clevs = np.arange(0., 16., 0.1) |
|---|
| 131 | cs = plt.contourf(xgrid_cart,ygrid_cart, sigmamean, clevs, cmap=cm.s3pcpn_l_r) |
|---|
| 132 | cbar = colorbar(cs) |
|---|
| 133 | xlim(-2500., 3000.) |
|---|
| 134 | ylim(-3000., 3000.) |
|---|
| 135 | |
|---|
| 136 | |
|---|
| 137 | |
|---|
| 138 | |
|---|
| 139 | |
|---|
| 140 | |
|---|
| 141 | #for canal in range(0,4): |
|---|
| 142 | # for ii in range (0,size(emis[canal,:])): |
|---|
| 143 | # if (emis[canal,ii] > 1): emis[canal,ii] = NaN |
|---|
| 144 | bblat_nov2009 = nonzero(lat >= 65.) |
|---|
| 145 | #bbfov = nonzero(fov[bblat] == 0.) |
|---|
| 146 | #ijr = 0 |
|---|
| 147 | bbjr_nov2009 = nonzero(jjr[bblat_nov2009] == 15.) |
|---|
| 148 | lon_zon_nov2009 = lon[bblat_nov2009][bbjr_nov2009] |
|---|
| 149 | lat_zon_nov2009 = lat[bblat_nov2009][bbjr_nov2009] |
|---|
| 150 | #emisch4 = emis[3,:][bblat][bbjr] |
|---|
| 151 | #emisch3 = emis[2,:][bblat][bbjr] |
|---|
| 152 | tbch1_nov2009 = tb[0,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 153 | tbch2_nov2009 = tb[1,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 154 | tbch3_nov2009 = tb[2,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 155 | tbch4_nov2009 = tb[3,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 156 | #tsurf = ts[bblat][bbjr] |
|---|
| 157 | #tdnch4 = tdn[3,:][bblat][bbjr] |
|---|
| 158 | #tdnch3 = tdn[2,:][bblat][bbjr] |
|---|
| 159 | #tupch4 = tup[3,:][bblat][bbjr] |
|---|
| 160 | #tupch3 = tup[2,:][bblat][bbjr] |
|---|
| 161 | |
|---|
| 162 | |
|---|
| 163 | |
|---|
| 164 | # DECEMBER2009 |
|---|
| 165 | imo = 1 |
|---|
| 166 | fichier=open("/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_EMIS_GLACE_"+str(le_mois[imo])+".DAT",'r') |
|---|
| 167 | numlines = 0 |
|---|
| 168 | for line in fichier: numlines += 1 |
|---|
| 169 | |
|---|
| 170 | fichier.close() |
|---|
| 171 | fichier=open("/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_EMIS_GLACE_"+str(le_mois[imo])+".DAT",'r') |
|---|
| 172 | nbtotal=numlines-1 |
|---|
| 173 | iligne=0 |
|---|
| 174 | lat=np.zeros([nbtotal],float) |
|---|
| 175 | lon=np.zeros([nbtotal],float) |
|---|
| 176 | jjr=np.zeros([nbtotal],float) |
|---|
| 177 | zen=np.zeros([nbtotal],float) |
|---|
| 178 | fov=np.zeros([nbtotal],float) |
|---|
| 179 | ts=np.zeros([nbtotal],float) |
|---|
| 180 | emis=np.zeros([4,nbtotal],float) |
|---|
| 181 | tb=np.zeros([4,nbtotal],float) |
|---|
| 182 | tup=np.zeros([4,nbtotal],float) |
|---|
| 183 | tdn=np.zeros([4,nbtotal],float) |
|---|
| 184 | #trans=np.zeros([4,nbtotal],float) |
|---|
| 185 | while (iligne < nbtotal) : |
|---|
| 186 | line = fichier.readline() |
|---|
| 187 | # exemple : line = "0.22 2.3 5.0 6" |
|---|
| 188 | liste = line.split() |
|---|
| 189 | # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es) |
|---|
| 190 | lat[iligne] = float(liste[1]) |
|---|
| 191 | lon[iligne] = float(liste[0]) |
|---|
| 192 | jjr[iligne] = float(liste[4]) |
|---|
| 193 | ts[iligne] = float(liste[8]) |
|---|
| 194 | fov[iligne] = float(liste[10]) |
|---|
| 195 | #trans[iligne] = float(liste[17]) |
|---|
| 196 | for ki in range(0,4): |
|---|
| 197 | emis[ki,iligne] = float(liste[12+ki]) |
|---|
| 198 | tb[ki,iligne] = float(liste[16+ki]) |
|---|
| 199 | tdn[ki,iligne] = float(liste[20+ki]) |
|---|
| 200 | tup[ki,iligne] = float(liste[24+ki]) |
|---|
| 201 | #trans[ki,iligne] = float(liste[28+ki]) |
|---|
| 202 | #tup_ssmis[ki,iligne] = float(liste[28+ki]) |
|---|
| 203 | #tdn_ssmis[ki,iligne] = float(liste[37+ki]) |
|---|
| 204 | #trans_ssmis[ki,iligne] = float(liste[44+ki]) |
|---|
| 205 | iligne=iligne+1 |
|---|
| 206 | |
|---|
| 207 | fichier.close() |
|---|
| 208 | #for canal in range(0,4): |
|---|
| 209 | # for ii in range (0,size(emis[canal,:])): |
|---|
| 210 | # if (emis[canal,ii] > 1): emis[canal,ii] = NaN |
|---|
| 211 | bblat_dec2009 = nonzero(lat >= 65.) |
|---|
| 212 | #bbfov = nonzero(fov[bblat] == 0.) |
|---|
| 213 | #ijr = 0 |
|---|
| 214 | bbjr_dec2009 = nonzero(jjr[bblat_dec2009] == 15.) |
|---|
| 215 | lon_zon_dec2009 = lon[bblat_dec2009][bbjr_dec2009] |
|---|
| 216 | lat_zon_dec2009 = lat[bblat_dec2009][bbjr_dec2009] |
|---|
| 217 | #emisch4 = emis[3,:][bblat][bbjr] |
|---|
| 218 | #emisch3 = emis[2,:][bblat][bbjr] |
|---|
| 219 | tbch1_dec2009 = tb[0,:][bblat_dec2009][bbjr_dec2009] |
|---|
| 220 | tbch2_dec2009 = tb[1,:][bblat_dec2009][bbjr_dec2009] |
|---|
| 221 | tbch3_dec2009 = tb[2,:][bblat_dec2009][bbjr_dec2009] |
|---|
| 222 | tbch4_dec2009 = tb[3,:][bblat_dec2009][bbjr_dec2009] |
|---|
| 223 | #tsurf = ts[bblat][bbjr] |
|---|
| 224 | #tdnch4 = tdn[3,:][bblat][bbjr] |
|---|
| 225 | #tdnch3 = tdn[2,:][bblat][bbjr] |
|---|
| 226 | #tupch4 = tup[3,:][bblat][bbjr] |
|---|
| 227 | #tupch3 = tup[2,:][bblat][bbjr] |
|---|
| 228 | |
|---|
| 229 | |
|---|
| 230 | # JANUARY2009 |
|---|
| 231 | imo = 2 |
|---|
| 232 | fichier=open("/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_EMIS_GLACE_"+str(le_mois[imo])+".DAT",'r') |
|---|
| 233 | numlines = 0 |
|---|
| 234 | for line in fichier: numlines += 1 |
|---|
| 235 | |
|---|
| 236 | fichier.close() |
|---|
| 237 | fichier=open("/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_EMIS_GLACE_"+str(le_mois[imo])+".DAT",'r') |
|---|
| 238 | nbtotal=numlines-1 |
|---|
| 239 | iligne=0 |
|---|
| 240 | lat=np.zeros([nbtotal],float) |
|---|
| 241 | lon=np.zeros([nbtotal],float) |
|---|
| 242 | jjr=np.zeros([nbtotal],float) |
|---|
| 243 | zen=np.zeros([nbtotal],float) |
|---|
| 244 | fov=np.zeros([nbtotal],float) |
|---|
| 245 | ts=np.zeros([nbtotal],float) |
|---|
| 246 | emis=np.zeros([4,nbtotal],float) |
|---|
| 247 | tb=np.zeros([4,nbtotal],float) |
|---|
| 248 | tup=np.zeros([4,nbtotal],float) |
|---|
| 249 | tdn=np.zeros([4,nbtotal],float) |
|---|
| 250 | #trans=np.zeros([4,nbtotal],float) |
|---|
| 251 | while (iligne < nbtotal) : |
|---|
| 252 | line = fichier.readline() |
|---|
| 253 | # exemple : line = "0.22 2.3 5.0 6" |
|---|
| 254 | liste = line.split() |
|---|
| 255 | # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es) |
|---|
| 256 | lat[iligne] = float(liste[1]) |
|---|
| 257 | lon[iligne] = float(liste[0]) |
|---|
| 258 | jjr[iligne] = float(liste[4]) |
|---|
| 259 | ts[iligne] = float(liste[8]) |
|---|
| 260 | fov[iligne] = float(liste[10]) |
|---|
| 261 | #trans[iligne] = float(liste[17]) |
|---|
| 262 | for ki in range(0,4): |
|---|
| 263 | emis[ki,iligne] = float(liste[12+ki]) |
|---|
| 264 | tb[ki,iligne] = float(liste[16+ki]) |
|---|
| 265 | tdn[ki,iligne] = float(liste[20+ki]) |
|---|
| 266 | tup[ki,iligne] = float(liste[24+ki]) |
|---|
| 267 | #trans[ki,iligne] = float(liste[28+ki]) |
|---|
| 268 | #tup_ssmis[ki,iligne] = float(liste[28+ki]) |
|---|
| 269 | #tdn_ssmis[ki,iligne] = float(liste[37+ki]) |
|---|
| 270 | #trans_ssmis[ki,iligne] = float(liste[44+ki]) |
|---|
| 271 | iligne=iligne+1 |
|---|
| 272 | |
|---|
| 273 | fichier.close() |
|---|
| 274 | #for canal in range(0,4): |
|---|
| 275 | # for ii in range (0,size(emis[canal,:])): |
|---|
| 276 | # if (emis[canal,ii] > 1): emis[canal,ii] = NaN |
|---|
| 277 | bblat_jan2009 = nonzero(lat >= 65.) |
|---|
| 278 | #bbfov = nonzero(fov[bblat] == 0.) |
|---|
| 279 | #ijr = 0 |
|---|
| 280 | bbjr_jan2009 = nonzero(jjr[bblat_jan2009] == 15.) |
|---|
| 281 | lon_zon_jan2009 = lon[bblat_jan2009][bbjr_jan2009] |
|---|
| 282 | lat_zon_jan2009 = lat[bblat_jan2009][bbjr_jan2009] |
|---|
| 283 | #emisch4 = emis[3,:][bblat][bbjr] |
|---|
| 284 | #emisch3 = emis[2,:][bblat][bbjr] |
|---|
| 285 | tbch1_jan2009 = tb[0,:][bblat_jan2009][bbjr_jan2009] |
|---|
| 286 | tbch2_jan2009 = tb[1,:][bblat_jan2009][bbjr_jan2009] |
|---|
| 287 | tbch3_jan2009 = tb[2,:][bblat_jan2009][bbjr_jan2009] |
|---|
| 288 | tbch4_jan2009 = tb[3,:][bblat_jan2009][bbjr_jan2009] |
|---|
| 289 | #tsurf = ts[bblat][bbjr] |
|---|
| 290 | #tdnch4 = tdn[3,:][bblat][bbjr] |
|---|
| 291 | #tdnch3 = tdn[2,:][bblat][bbjr] |
|---|
| 292 | #tupch4 = tup[3,:][bblat][bbjr] |
|---|
| 293 | #tupch3 = tup[2,:][bblat][bbjr] |
|---|
| 294 | |
|---|
| 295 | |
|---|
| 296 | |
|---|
| 297 | |
|---|
| 298 | |
|---|
| 299 | ## test of variance distribution with different spatial resolution ##################################### |
|---|
| 300 | imo = 0 |
|---|
| 301 | longitude = lon_zon_nov2009 |
|---|
| 302 | latitude = lat_zon_nov2009 |
|---|
| 303 | z = tbch1_nov2009 |
|---|
| 304 | z0 = min(z) |
|---|
| 305 | z1 = max(z) |
|---|
| 306 | delta_xy = np.array([25., 30., 35., 40., 45., 50.]) |
|---|
| 307 | Lxy = len(delta_xy) |
|---|
| 308 | plt.ion() |
|---|
| 309 | nber_val = np.zeros([6], float) |
|---|
| 310 | for idel in range (0, Lxy): |
|---|
| 311 | dx = delta_xy[idel] |
|---|
| 312 | dy = delta_xy[idel] |
|---|
| 313 | x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() |
|---|
| 314 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 315 | plt.figure() |
|---|
| 316 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 317 | clevs = np.arange(137., 269, 1.) |
|---|
| 318 | cs = plt.contourf(xgrid_cart,ygrid_cart, zgrid_output[:, :, 0], clevs, cmap=cm.s3pcpn_l_r) |
|---|
| 319 | cbar = colorbar(cs) |
|---|
| 320 | txt = 'variance - spatial resolution = ' + str(delta_xy[idel]) + 'km' |
|---|
| 321 | cbar.set_label(txt) |
|---|
| 322 | xlim(-4000., 3000.) |
|---|
| 323 | ylim(-4000., 4000.) |
|---|
| 324 | plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/variance_new_grid_spatial_res_'+ str(delta_xy[idel]) + 'km_' + le_mois[imo] + '.png') |
|---|
| 325 | ######################################################################################################## |
|---|
| 326 | |
|---|
| 327 | ## new grid parameters for mapping ## |
|---|
| 328 | dx=50. |
|---|
| 329 | dy=50. |
|---|
| 330 | longitude = lon_zon_nov2009 |
|---|
| 331 | latitude = lat_zon_nov2009 |
|---|
| 332 | z = tbch1_nov2009 |
|---|
| 333 | z0 = min(z) |
|---|
| 334 | z1 = max(z) |
|---|
| 335 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 336 | tb_mesh_nov2009 = ZGRID |
|---|
| 337 | x_mesh = xvec |
|---|
| 338 | y_mesh = yvec |
|---|
| 339 | ## map parameters for arctic coast limits mapping ## |
|---|
| 340 | x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() |
|---|
| 341 | ## mapping ## |
|---|
| 342 | plt.ion() |
|---|
| 343 | plt.figure() |
|---|
| 344 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 345 | clevs1 = np.arange(0, 30, 1.) |
|---|
| 346 | cs = plt.contourf(xgrid_cart,ygrid_cart, ngrid, clevs1, cmap=cm.s3pcpn_l_r) |
|---|
| 347 | cbar = colorbar(cs) |
|---|
| 348 | txt = 'number of values per grid - Tb ch1 (16th, november 2009) - AMSU-A' |
|---|
| 349 | cbar.set_label(txt) |
|---|
| 350 | xlim(-3000., 2500.) |
|---|
| 351 | ylim(-3000., 3000.) |
|---|
| 352 | |
|---|
| 353 | |
|---|
| 354 | |
|---|
| 355 | |
|---|
| 356 | ## calculation of standard deviation for november in each new_grid cell ## |
|---|
| 357 | ## 1 ## matrix of new gridded Tb per day #### |
|---|
| 358 | x0 = -3000. |
|---|
| 359 | x1 = 2500. |
|---|
| 360 | dx=25. |
|---|
| 361 | xvec = np.arange(x0, x1+dx, dx) |
|---|
| 362 | nx = len(xvec) |
|---|
| 363 | y0 = -3000. |
|---|
| 364 | y1 = 3000. |
|---|
| 365 | dy=25. |
|---|
| 366 | yvec = np.arange(y0, y1+dy, dy) |
|---|
| 367 | ny = len(yvec) |
|---|
| 368 | tbch1_mesh_nov2jan009 = np.zeros([ny, nx, 30], float) |
|---|
| 369 | tbch2_mesh_jan2009 = np.zeros([ny, nx, 30], float) |
|---|
| 370 | tsurf_mesh_jan2009 = np.zeros([ny, nx, 30], float) |
|---|
| 371 | for ijr in range (0, 30): |
|---|
| 372 | bbjr_jan2009 = nonzero(jjr[bblat_jan2009] == ijr + 1.) |
|---|
| 373 | print 'date: ' + str(ijr + 1) + ' nov 2009' |
|---|
| 374 | lon_zon_nov2009 = lon[bblat_nov2009][bbjr_nov2009] |
|---|
| 375 | lat_zon_nov2009 = lat[bblat_nov2009][bbjr_nov2009] |
|---|
| 376 | tbch1_nov2009 = tb[0,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 377 | #tbch2_nov2009 = tb[1,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 378 | #tbch3_nov2009 = tb[2,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 379 | #tbch4_nov2009 = tb[3,:][bblat_nov2009][bbjr_nov2009] |
|---|
| 380 | tsurf_nov2009 = ts[bblat_nov2009][bbjr_nov2009] |
|---|
| 381 | longitude = lon_zon_nov2009 |
|---|
| 382 | latitude = lat_zon_nov2009 |
|---|
| 383 | z = tbch1_nov2009 # channel 1 |
|---|
| 384 | z0 = min(z) |
|---|
| 385 | z1 = max(z) |
|---|
| 386 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 387 | tbch1_mesh_nov2009[:, :, ijr] = ZGRID |
|---|
| 388 | print 'channel 1' |
|---|
| 389 | x_mesh = xvec |
|---|
| 390 | y_mesh = yvec |
|---|
| 391 | z = tbch2_nov2009 # channel 2 |
|---|
| 392 | z0 = min(z) |
|---|
| 393 | z1 = max(z) |
|---|
| 394 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 395 | tbch2_mesh_nov2009[:, :, ijr] = ZGRID |
|---|
| 396 | print 'channel 2' |
|---|
| 397 | z = tsurf_nov2009 # surface temperature |
|---|
| 398 | z0 = min(z) |
|---|
| 399 | z1 = max(z) |
|---|
| 400 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 401 | tsurf_mesh_nov2009[:, :, ijr] = ZGRID |
|---|
| 402 | print 'surf temp' |
|---|
| 403 | |
|---|
| 404 | |
|---|
| 405 | |
|---|
| 406 | ## 2 ## std in each grid cell #### |
|---|
| 407 | tbsum_tbch1 = np.zeros([241, 261], float) |
|---|
| 408 | tbsum_tbch2 = np.zeros([241, 261], float) |
|---|
| 409 | tbsum_tsurf = np.zeros([241, 261], float) |
|---|
| 410 | for ix in range (0, 261): |
|---|
| 411 | for iy in range (0, 241): |
|---|
| 412 | bb_value_tbch1 = nonzero(isnan(tbch1_mesh_nov2009[iy, ix, :]) == False) |
|---|
| 413 | bb_value_tbch2 = nonzero(isnan(tbch2_mesh_nov2009[iy, ix, :]) == False) |
|---|
| 414 | bb_value_tsurf = nonzero(isnan(tsurf_mesh_nov2009[iy, ix, :]) == False) |
|---|
| 415 | # channel 1 |
|---|
| 416 | if (bb_value_tbch1[0] == np.array([])): |
|---|
| 417 | tbsum_tbch1[iy, ix] = NaN # take out grid cells where no obs in the month |
|---|
| 418 | else: |
|---|
| 419 | tbsum_tbch1[iy, ix] = sum((tbch1_mesh_nov2009[iy, ix, :][bb_value_tbch1] - mean(tbch1_mesh_nov2009[iy, ix, :][bb_value_tbch1])) ** 2) |
|---|
| 420 | # channel 2 |
|---|
| 421 | if (bb_value_tbch2[0] == np.array([])): |
|---|
| 422 | tbsum_tbch2[iy, ix] = NaN # take out grid cells where no obs in the month |
|---|
| 423 | else: |
|---|
| 424 | tbsum_tbch2[iy, ix] = sum((tbch2_mesh_nov2009[iy, ix, :][bb_value_tbch2] - mean(tbch2_mesh_nov2009[iy, ix, :][bb_value_tbch2])) ** 2) |
|---|
| 425 | # surface temperature |
|---|
| 426 | if (bb_value_tsurf[0] == np.array([])): |
|---|
| 427 | tbsum_tsurf[iy, ix] = NaN # take out grid cells where no obs in the month |
|---|
| 428 | else: |
|---|
| 429 | tbsum_tsurf[iy, ix] = sum((tsurf_mesh_nov2009[iy, ix, :][bb_value_tsurf] - mean(tsurf_mesh_nov2009[iy, ix, :][bb_value_tsurf])) ** 2) |
|---|
| 430 | |
|---|
| 431 | |
|---|
| 432 | |
|---|
| 433 | std_tbch1_nov2009 = sqrt((1. / (30. - 1.)) * tbsum_tbch1) |
|---|
| 434 | std_tbch2_nov2009 = sqrt((1. / (30. - 1.)) * tbsum_tbch2) |
|---|
| 435 | std_tsurf_nov2009 = sqrt((1. / (30. - 1.)) * tbsum_tsurf) |
|---|
| 436 | |
|---|
| 437 | for ix in range (0, 261): |
|---|
| 438 | for iy in range (0, 241): |
|---|
| 439 | if (std_tbch1_nov2009[iy, ix] < 0.): |
|---|
| 440 | std_tbch1_nov2009[iy, ix] == NaN # take out the negative values of std |
|---|
| 441 | if (std_tbch2_nov2009[iy, ix] < 0.): |
|---|
| 442 | std_tbch2_nov2009[iy, ix] == NaN # take out the negative values of std |
|---|
| 443 | if (std_tsurf_nov2009[iy, ix] < 0.): |
|---|
| 444 | std_tsurf_nov2009[iy, ix] == NaN # take out the negative values of std |
|---|
| 445 | |
|---|
| 446 | |
|---|
| 447 | tbch1_mean_nov2009 = np.zeros([241, 261], float) |
|---|
| 448 | tbch2_mean_nov2009 = np.zeros([241, 261], float) |
|---|
| 449 | tsurf_mean_nov2009 = np.zeros([241, 261], float) |
|---|
| 450 | for ix in range (0, 261): |
|---|
| 451 | for iy in range (0, 241): |
|---|
| 452 | bb_val_tbch1 = nonzero(isnan(tbch1_mesh_nov2009[iy, ix, :]) == False) |
|---|
| 453 | bb_val_tbch2 = nonzero(isnan(tbch2_mesh_nov2009[iy, ix, :]) == False) |
|---|
| 454 | bb_val_tsurf = nonzero(isnan(tsurf_mesh_nov2009[iy, ix, :]) == False) |
|---|
| 455 | tbch1_mean_nov2009[iy, ix] = tbch1_mesh_nov2009[iy, ix, :][bb_val_tbch1].mean() |
|---|
| 456 | tbch2_mean_nov2009[iy, ix] = tbch2_mesh_nov2009[iy, ix, :][bb_val_tbch2].mean() |
|---|
| 457 | tsurf_mean_nov2009[iy, ix] = tsurf_mesh_nov2009[iy, ix, :][bb_val_tsurf].mean() |
|---|
| 458 | |
|---|
| 459 | ## study zone ## |
|---|
| 460 | ## 1 ## monthly mean (temporal) values #### |
|---|
| 461 | study_x = np.where((xvec >= -800.) & (xvec <= -350.))[0] |
|---|
| 462 | study_y = np.where((yvec >= 950.) & (yvec <= 1300.))[0] |
|---|
| 463 | Sx = len(study_x) |
|---|
| 464 | Sy = len(study_y) |
|---|
| 465 | tbch1_mean_zone_nov2009 = tbch1_mean_nov2009[study_y[0]: study_y[Sy - 1] + 1, study_x[0]: study_x[Sx - 1] + 1] |
|---|
| 466 | tbch2_mean_zone_nov2009 = tbch2_mean_nov2009[study_y[0]: study_y[Sy - 1] + 1, study_x[0]: study_x[Sx - 1] + 1] |
|---|
| 467 | tsurf_mean_zone_nov2009 = tsurf_mean_nov2009[study_y[0]: study_y[Sy - 1] + 1, study_x[0]: study_x[Sx - 1] + 1] |
|---|
| 468 | ## 2 ## daily mean (spatial) values #### |
|---|
| 469 | tbch1_mean_day_nov2009 = np.zeros([30], float) |
|---|
| 470 | tbch2_mean_day_nov2009 = np.zeros([30], float) |
|---|
| 471 | tsurf_mean_day_nov2009 = np.zeros([30], float) |
|---|
| 472 | for ijr in range (0, 30): |
|---|
| 473 | a = np.reshape(tbch1_mesh_nov2009[:, :, ijr], size(tbch1_mesh_nov2009[:, :, ijr])) |
|---|
| 474 | bba = nonzero(isnan(a) == False) |
|---|
| 475 | moy_a = mean(a[bba]) |
|---|
| 476 | tbch1_mean_day_nov2009[ijr] = moy_a |
|---|
| 477 | b = np.reshape(tbch2_mesh_nov2009[:, :, ijr], size(tbch2_mesh_nov2009[:, :, ijr])) |
|---|
| 478 | bbb = nonzero(isnan(b) == False) |
|---|
| 479 | moy_b = mean(b[bbb]) |
|---|
| 480 | tbch2_mean_day_nov2009[ijr] = moy_b |
|---|
| 481 | c = np.reshape(tsurf_mesh_nov2009[:, :, ijr], size(tsurf_mesh_nov2009[:, :, ijr])) |
|---|
| 482 | bbc = nonzero(isnan(c) == False) |
|---|
| 483 | moy_c = mean(c[bbc]) |
|---|
| 484 | tsurf_mean_day_nov2009[ijr] = moy_c |
|---|
| 485 | anom_tbch1_nov2009[ijr] |
|---|
| 486 | |
|---|
| 487 | |
|---|
| 488 | plt.figure() |
|---|
| 489 | plt.plot(tbch1_mean_day_nov2009, 'b', label = 'Tb ch1') |
|---|
| 490 | plt.plot(tbch2_mean_day_nov2009, 'g', label = 'Tb ch2') |
|---|
| 491 | legend() |
|---|
| 492 | twinx() |
|---|
| 493 | plt.plot(tsurf_mean_day_nov2009, 'r', label = 'Tsurf') |
|---|
| 494 | legend() |
|---|
| 495 | |
|---|
| 496 | |
|---|
| 497 | # channel 1 |
|---|
| 498 | plt.ion() |
|---|
| 499 | plt.figure() |
|---|
| 500 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 501 | clevs1 = np.arange(140., 255., 1.) |
|---|
| 502 | cs = plt.contourf(xgrid_cart,ygrid_cart, tbch1_mesh_nov2009[:, :, ijr], clevs1, cmap=cm.s3pcpn_l_r) |
|---|
| 503 | cbar = colorbar(cs) |
|---|
| 504 | txt = 'std Tb ch1 (november 2009) - AMSU-A' |
|---|
| 505 | cbar.set_label(txt) |
|---|
| 506 | xlim(-3000., 2500.) |
|---|
| 507 | ylim(-2500., 3500.) |
|---|
| 508 | |
|---|
| 509 | # channel 2 |
|---|
| 510 | plt.figure() |
|---|
| 511 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 512 | clevs1 = np.arange(0.8, 39, 1.) |
|---|
| 513 | cs = plt.contourf(xgrid_cart,ygrid_cart,std_tbch2_nov2009, clevs1, cmap=cm.s3pcpn_l_r) |
|---|
| 514 | cbar = colorbar(cs) |
|---|
| 515 | txt = 'std Tb ch2 (november 2009) - AMSU-A' |
|---|
| 516 | cbar.set_label(txt) |
|---|
| 517 | xlim(-4000., 3000.) |
|---|
| 518 | ylim(-4000., 4000.) |
|---|
| 519 | |
|---|
| 520 | # surface temperature |
|---|
| 521 | plt.figure() |
|---|
| 522 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 523 | clevs1 = np.arange(0., 10., 0.1) |
|---|
| 524 | cs = plt.contourf(xgrid_cart,ygrid_cart,std_tsurf_nov2009, clevs1, cmap=cm.s3pcpn_l_r) |
|---|
| 525 | cbar = colorbar(cs) |
|---|
| 526 | txt = 'std Tsurf (november 2009) - AMSU-A' |
|---|
| 527 | cbar.set_label(txt) |
|---|
| 528 | xlim(-4000., 3000.) |
|---|
| 529 | ylim(-4000., 4000.) |
|---|
| 530 | |
|---|
| 531 | |
|---|
| 532 | |
|---|
| 533 | |
|---|
| 534 | |
|---|
| 535 | |
|---|
| 536 | longitude = lon_zon_nov2009 |
|---|
| 537 | latitude = lat_zon_nov2009 |
|---|
| 538 | |
|---|
| 539 | delta_xy = np.array([25., 30., 35., 40., 45., 50.]) |
|---|
| 540 | idel = 0 |
|---|
| 541 | x0 = -3000. |
|---|
| 542 | x1 = 2500. |
|---|
| 543 | dx = delta_xy[idel] |
|---|
| 544 | xvec = np.arange(x0, x1+dx, dx) |
|---|
| 545 | nx = len(xvec) |
|---|
| 546 | y0 = -3000. |
|---|
| 547 | y1 = 3000. |
|---|
| 548 | dy= delta_xy[idel] |
|---|
| 549 | yvec = np.arange(y0, y1+dy, dy) |
|---|
| 550 | ny = len(yvec) |
|---|
| 551 | for ix in range (0, nx): |
|---|
| 552 | for iy in range (0, ny): |
|---|
| 553 | tbch1_moy_nov2009[iy, ix] = tbch1_mesh_nov2009[iy, ix, :].mean |
|---|
| 554 | z = tbch1_moy_nov2009 |
|---|
| 555 | z0 = min(z) |
|---|
| 556 | z1 = max(z) |
|---|
| 557 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 558 | moy_ZGRID = ZGRID[nonzero(isnan(ZGRID) == False)].mean() |
|---|
| 559 | good_points = np.zeros([ny, nx], float) |
|---|
| 560 | for ix in range (0, nx): |
|---|
| 561 | for iy in range (0, ny): |
|---|
| 562 | if (sigmagrid[iy, ix] > ((1. * moy_ZGRID) / 100.)): |
|---|
| 563 | good_points[iy, ix] = NaN |
|---|
| 564 | else: |
|---|
| 565 | good_points[iy, ix] = sigmagrid[iy, ix] |
|---|
| 566 | |
|---|
| 567 | |
|---|
| 568 | plt.figure() |
|---|
| 569 | plt.scatter(x_ind, y_ind, c = z_ind, s = volume) |
|---|
| 570 | clevs1 = np.arange(143, 255, 1) |
|---|
| 571 | cs = plt.contourf(xgrid_cart,ygrid_cart,ZGRID, clevs1, cmap=cm.s3pcpn_l_r) |
|---|
| 572 | cbar = colorbar(cs) |
|---|
| 573 | txt = 'best_grid cells: std < 1% of Tbch1 - res' + str(delta_xy[idel]) + 'km (november 2009, AMSU-A)' |
|---|
| 574 | cbar.set_label(txt) |
|---|
| 575 | xlim(-3000., 2500.) |
|---|
| 576 | ylim(-3000., 3000.) |
|---|
| 577 | |
|---|
| 578 | |
|---|
| 579 | |
|---|
| 580 | |
|---|
| 581 | |
|---|
| 582 | |
|---|
| 583 | |
|---|
| 584 | |
|---|
| 585 | |
|---|
| 586 | |
|---|
| 587 | nber_val_res25 = ngrid |
|---|
| 588 | idel = 1 |
|---|
| 589 | dx = delta_xy[idel] |
|---|
| 590 | dy = delta_xy[idel] |
|---|
| 591 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 592 | nber_val_res30 = ngrid |
|---|
| 593 | idel = 2 |
|---|
| 594 | dx = delta_xy[idel] |
|---|
| 595 | dy = delta_xy[idel] |
|---|
| 596 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 597 | nber_val_res35 = ngrid |
|---|
| 598 | idel = 3 |
|---|
| 599 | dx = delta_xy[idel] |
|---|
| 600 | dy = delta_xy[idel] |
|---|
| 601 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 602 | nber_val_res40 = ngrid |
|---|
| 603 | idel = 4 |
|---|
| 604 | dx = delta_xy[idel] |
|---|
| 605 | dy = delta_xy[idel] |
|---|
| 606 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 607 | nber_val_res45 = ngrid |
|---|
| 608 | idel = 5 |
|---|
| 609 | dx = delta_xy[idel] |
|---|
| 610 | dy = delta_xy[idel] |
|---|
| 611 | ZGRID, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_test.new_cartesian_grid(longitude, latitude, z, z0, z1, dx, dy) |
|---|
| 612 | nber_val_res50 = ngrid |
|---|
| 613 | |
|---|
| 614 | |
|---|
| 615 | ngrid_25km_vec = reshape(nber_val_res25, size(nber_val_res25)) |
|---|
| 616 | ngrid_30km_vec = reshape(nber_val_res30, size(nber_val_res30)) |
|---|
| 617 | ngrid_35km_vec = reshape(nber_val_res35, size(nber_val_res35)) |
|---|
| 618 | ngrid_40km_vec = reshape(nber_val_res40, size(nber_val_res40)) |
|---|
| 619 | ngrid_45km_vec = reshape(nber_val_res45, size(nber_val_res45)) |
|---|
| 620 | ngrid_50km_vec = reshape(nber_val_res50, size(nber_val_res50)) |
|---|
| 621 | |
|---|
| 622 | plt.figure() |
|---|
| 623 | hist([ngrid_25km_vec, ngrid_30km_vec, ngrid_35km_vec, ngrid_40km_vec, ngrid_45km_vec, ngrid_50km_vec], histtype = 'bar', label = ['25km', '30km', '35km', '40km', '45km', '50km']) |
|---|
| 624 | |
|---|
| 625 | |
|---|
![(please configure the [header_logo] section in trac.ini)](/twice/chrome/site/your_project_logo.png){kind=logo}
