| 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 | |
|---|
| 16 | |
|---|
| 17 | |
|---|
| 18 | ##################### |
|---|
| 19 | # read fichier .DAT # |
|---|
| 20 | ##################### |
|---|
| 21 | fichier=open('/mma/hermozol/Documents/Data_tests/GLACE/GLACE_AMSUB89_EMIS_SEPTEMBER2009.DAT','r') |
|---|
| 22 | numlines = 0 |
|---|
| 23 | for line in fichier: numlines += 1 |
|---|
| 24 | |
|---|
| 25 | fichier.close |
|---|
| 26 | |
|---|
| 27 | fichier=open('/mma/hermozol/Documents/Data_tests/GLACE/GLACE_AMSUB89_EMIS_SEPTEMBER2009.DAT','r') |
|---|
| 28 | |
|---|
| 29 | nbtotal=numlines-1 |
|---|
| 30 | |
|---|
| 31 | iligne=0 |
|---|
| 32 | lat=np.zeros([nbtotal],float) |
|---|
| 33 | lon=np.zeros([nbtotal],float) |
|---|
| 34 | jjr=np.zeros([nbtotal],float) |
|---|
| 35 | zen=np.zeros([nbtotal],float) |
|---|
| 36 | fov=np.zeros([nbtotal],float) |
|---|
| 37 | ts=np.zeros([nbtotal],float) |
|---|
| 38 | emis=np.zeros([nbtotal],float) |
|---|
| 39 | tb=np.zeros([nbtotal],float) |
|---|
| 40 | tup=np.zeros([nbtotal],float) |
|---|
| 41 | tdn=np.zeros([nbtotal],float) |
|---|
| 42 | trans=np.zeros([nbtotal],float) |
|---|
| 43 | |
|---|
| 44 | # format: |
|---|
| 45 | # printf, luna2, xlon, ylat, annee, mois, jour, heure, jour_obs, heure_obs, tsy, zeny, fovy, saty, masque, ee89, tb89, tup89, tdn89, trans89, format='(18(f12.3,1X))' |
|---|
| 46 | |
|---|
| 47 | while (iligne < nbtotal) : |
|---|
| 48 | line=fichier.readline() |
|---|
| 49 | # exemple : line = "0.22 2.3 5.0 6" |
|---|
| 50 | liste = line.split() |
|---|
| 51 | # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es) |
|---|
| 52 | lat[iligne] = float(liste[1]) |
|---|
| 53 | lon[iligne] = float(liste[0]) |
|---|
| 54 | jjr[iligne] = float(liste[4]) |
|---|
| 55 | ts[iligne] = float(liste[8]) |
|---|
| 56 | fov[iligne] = float(liste[10]) |
|---|
| 57 | zen[iligne] = float(liste[9]) |
|---|
| 58 | emis[iligne] = float(liste[13]) |
|---|
| 59 | tb[iligne] = float(liste[14]) |
|---|
| 60 | tup[iligne] = float(liste[15]) |
|---|
| 61 | tdn[iligne] = float(liste[16]) |
|---|
| 62 | trans[iligne] = float(liste[17]) |
|---|
| 63 | iligne=iligne+1 |
|---|
| 64 | |
|---|
| 65 | |
|---|
| 66 | fichier.close() |
|---|
| 67 | |
|---|
| 68 | for i in range (0, nbtotal): |
|---|
| 69 | if (emis[i] < 0.): |
|---|
| 70 | emis[i] = NaN |
|---|
| 71 | if (emis[i] > 1.): |
|---|
| 72 | emis[i] = 0.99 |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | |
|---|
| 76 | ############################ |
|---|
| 77 | # definition of parameters # |
|---|
| 78 | ############################ |
|---|
| 79 | bblat = nonzero(lat >= 65.) # only high latitudes |
|---|
| 80 | llon = lon[bblat] # longitude |
|---|
| 81 | llat = lat[bblat] # latitude |
|---|
| 82 | zen_angle = zen[bblat]# zenith angle |
|---|
| 83 | field_view = fov[bblat] # field of view |
|---|
| 84 | e = emis[bblat] # emissivity |
|---|
| 85 | t_surf = ts[bblat] # surface temperature |
|---|
| 86 | t_up = tup[bblat] # tup |
|---|
| 87 | t_dn = tdn[bblat] # tdown |
|---|
| 88 | transm = trans[bblat] # transmissivity |
|---|
| 89 | t_b = tb[bblat] # brightness temperature |
|---|
| 90 | L = shape(bblat)[1] |
|---|
| 91 | |
|---|
| 92 | |
|---|
| 93 | ####################################### |
|---|
| 94 | # calculation of lamb_spec parameters # |
|---|
| 95 | ####################################### |
|---|
| 96 | r = pi / 180. |
|---|
| 97 | opac = - cos(zen_angle * r) * np.log(transm) |
|---|
| 98 | E3 = scipy.special.expn(3, opac) |
|---|
| 99 | theta_eff = np.zeros([L], float) |
|---|
| 100 | tdn_lamb = np.zeros([L], float) |
|---|
| 101 | tb_spec = np.zeros([L], float) |
|---|
| 102 | tb_lamb = np.zeros([L], float) |
|---|
| 103 | e_spec = np.zeros([L], float) |
|---|
| 104 | e_lamb = np.zeros([L], float) |
|---|
| 105 | e_spec_lamb = np.zeros([L], float) |
|---|
| 106 | for ii in range (0, L): |
|---|
| 107 | theta_eff[ii] = math.acos(- opac[ii] / (np.log(2. * E3[ii]))) * (1. / r) |
|---|
| 108 | tdn_lamb[ii] = t_surf[ii] * (1 - exp(- opac[ii] / cos (theta_eff[ii] * r))) # Tdown to calculate with theta_eff |
|---|
| 109 | tb_spec[ii] = (t_surf[ii] * e[ii] * exp(- opac[ii])) + (1. - e[ii]) * exp(- opac[ii]) * t_dn[ii] + t_up[ii] # calculation of Tb spec with theta |
|---|
| 110 | tb_lamb[ii] = (t_surf[ii] * e[ii] * exp(- opac[ii])) + (1. - e[ii]) * exp(- opac[ii]) * tdn_lamb[ii] + t_up[ii] # calculation of Tb spec with theta_eff (with Tdown_lamb) |
|---|
| 111 | e_spec[ii] = ((tb_spec[ii] - t_up[ii]) * exp(opac[ii]) - t_dn[ii]) / (t_surf[ii] - t_dn[ii]) |
|---|
| 112 | e_lamb[ii] = ((tb_lamb[ii] - t_up[ii]) * exp(opac[ii]) - t_dn[ii]) / (t_surf[ii] - t_dn[ii]) |
|---|
| 113 | e_spec_lamb[ii] = e_lamb[ii] - e_spec[ii] |
|---|
| 114 | |
|---|
| 115 | |
|---|
| 116 | |
|---|
| 117 | plt.ion() |
|---|
| 118 | plt.figure() |
|---|
| 119 | plt.scatter(opac, theta_eff) |
|---|
| 120 | xlim(0.14, 0.42) |
|---|
| 121 | ylim(52.9, 56.2) |
|---|
| 122 | xticks(np.arange(0.14, 0.44, 0.02)) |
|---|
| 123 | grid(True) |
|---|
| 124 | xlabel('zenith opacity') |
|---|
| 125 | ylabel('effective incidence angle (deg)') |
|---|
| 126 | title('SEPTEMBER 2009 - AMSUB') |
|---|
| 127 | plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/scatter_theta_eff-opacity_sept2009_AMSUB.png') |
|---|
| 128 | |
|---|
| 129 | |
|---|
| 130 | plt.figure() |
|---|
| 131 | plt.scatter(field_view, e_spec_lamb) |
|---|
| 132 | xlim(0., 90.) |
|---|
| 133 | xticks(np.arange(0, 90, 5)) |
|---|
| 134 | grid(True) |
|---|
| 135 | xlabel('scan position') |
|---|
| 136 | ylabel('emissivity SPEC-LAMB') |
|---|
| 137 | title('SEPTEMBER 2009 - AMSUB') |
|---|
| 138 | plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/scatter_e-spec_e-lamb_fov_sept2009_AMSUB.png') |
|---|
| 139 | |
|---|
| 140 | |
|---|
| 141 | |
|---|
| 142 | e_diff_mean = np.zeros([89], float) |
|---|
| 143 | e_diff_ect = np.zeros([89], float) |
|---|
| 144 | for ifov in range (0, 89): |
|---|
| 145 | ind = np.where(field_view == float(ifov))[0] |
|---|
| 146 | if (ind == []): |
|---|
| 147 | e_diff_mean[ifov] = NaN |
|---|
| 148 | e_diff_ect[ifov] = NaN |
|---|
| 149 | else: |
|---|
| 150 | e_diff_mean[ifov] = mean(e_spec_lamb[ind]) |
|---|
| 151 | e_diff_ect[ifov] = sqrt(sum((e_spec_lamb[ind] - e_diff_mean[ifov]) ** 2.)) |
|---|
| 152 | if (e_diff_ect[ifov] == 0.): |
|---|
| 153 | e_diff_ect[ifov] = NaN |
|---|
| 154 | |
|---|
| 155 | |
|---|
| 156 | |
|---|
| 157 | |
|---|
| 158 | plt.figure() |
|---|
| 159 | plt.plot(e_diff_mean, 'bo', label = 'mean value') |
|---|
| 160 | xlabel('field of view') |
|---|
| 161 | ylabel('emissivity SPEC-LAMB') |
|---|
| 162 | legend(loc = 2) |
|---|
| 163 | twinx() |
|---|
| 164 | plot(e_diff_ect , 'go', label = 'std value') |
|---|
| 165 | xlabel('standard deviation') |
|---|
| 166 | xticks(np.arange(0, 90, 10)) |
|---|
| 167 | legend(loc = 1) |
|---|
| 168 | plt.title('SEPTEMBER 2009 - AMSUB') |
|---|
| 169 | plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/evo_emis_lamb-spec_fov_sept2009_AMSUB.png') |
|---|
| 170 | |
|---|
| 171 | |
|---|
| 172 | |
|---|
| 173 | |
|---|
| 174 | |
|---|
| 175 | ############################################## |
|---|
| 176 | # stack of lamb_spec parameters in .txt file # |
|---|
| 177 | ############################################## |
|---|
| 178 | data_spec_lamb = open ('/mma/hermozol/Documents/Data_tests/monthly_GLACE/lamb_spec_param_all_angles_SEPTEMBER2009.dat', 'a') |
|---|
| 179 | for ii in range (0, L): |
|---|
| 180 | data_spec_lamb.write(('%(lon)4.2f %(lat)4.2f %(za)4.2f %(fov)2.0f %(emis)5.4f %(tsurf)5.2f %(tup)5.2f %(tdn)5.2f %(trans)5.4f %(tb)5.2f %(opac)5.4f %(theta_eff)5.3f %(tdn_lamb)5.2f %(tb_spec)5.2f %(tb_lamb)5.2f %(e_spec)5.4f %(e_lamb)5.4f %(e_diff)5.4f \n' % { |
|---|
| 181 | 'lon':llon[ii], |
|---|
| 182 | 'lat':llat[ii], |
|---|
| 183 | 'za':zen_angle[ii], |
|---|
| 184 | 'fov':field_view[ii], |
|---|
| 185 | 'emis':e[ii], |
|---|
| 186 | 'tsurf':t_surf[ii], |
|---|
| 187 | 'tup':t_up[ii], |
|---|
| 188 | 'tdn':t_dn[ii], |
|---|
| 189 | 'trans':transm[ii], |
|---|
| 190 | 'tb':t_b[ii], |
|---|
| 191 | 'opac':opac[ii], |
|---|
| 192 | 'theta_eff':theta_eff[ii], |
|---|
| 193 | 'tdn_lamb':tdn_lamb[ii], |
|---|
| 194 | 'tb_spec':tb_spec[ii], |
|---|
| 195 | 'tb_lamb':tb_lamb[ii], |
|---|
| 196 | 'e_spec':e_spec[ii], |
|---|
| 197 | 'e_lamb':e_lamb[ii], |
|---|
| 198 | 'e_diff':e_spec_lamb[ii], |
|---|
| 199 | })) |
|---|
| 200 | |
|---|
| 201 | |
|---|
| 202 | data_spec_lamb.close() |
|---|
| 203 | |
|---|
| 204 | |
|---|
| 205 | |
|---|
| 206 | |
|---|
![(please configure the [header_logo] section in trac.ini)](/twice/chrome/site/your_project_logo.png){kind=logo}
