source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/calcul_spec_lamb_nadir_bis.py @ 55

Last change on this file since 55 was 44, checked in by lahlod, 10 years ago

modifs
File size: 7.8 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3import string
4import numpy as np
5import matplotlib.pyplot as plt
6from pylab import *
7from mpl_toolkits.basemap import Basemap
8from mpl_toolkits.basemap import shiftgrid, cm
9from 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
12import scipy.special
13
14
15
16#####################
17# read fichier .DAT #
18#####################
19month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER'])
20month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]), 
21M = len(month)
22
23
24
25
26for imo in range (0, M):
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_AMSUA50.dat', 'a')
29    ############################################
30    print 'read file' + month[imo]
31    fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/GLACE/AMSUA/GLACE_AMSUA_EMIS_' + month[imo] + '2009.DAT','r')
32    numlines = 0
33    for line in fichier: numlines += 1
34    fichier.close
35    fichier = open('/net/argos/data/parvati/lahlod/ARCTIC/GLACE/AMSUA/GLACE_AMSUA_EMIS_' + month[imo] + '2009.DAT','r')
36    nbtotal = numlines - 1
37    iligne = 0
38    lat = np.zeros([nbtotal],float)
39    lon = np.zeros([nbtotal],float)
40    jjr = np.zeros([nbtotal],float)
41    zen = np.zeros([nbtotal],float)
42    fov = np.zeros([nbtotal],float)
43    ts = np.zeros([nbtotal],float)
44    emis= np.zeros([4, nbtotal],float)
45    tb = np.zeros([4, nbtotal],float)
46    tup = np.zeros([4, nbtotal],float)
47    tdn = np.zeros([4, nbtotal],float)
48    trans = np.zeros([4, nbtotal],float)
49    # format:
50    # 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))' 
51    while (iligne < nbtotal) :
52         line=fichier.readline()
53         # exemple : line = "0.22 2.3 5.0 6"
54         liste = line.split()
55         # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es)
56         lat[iligne] = float(liste[1])
57         lon[iligne] = float(liste[0])
58         jjr[iligne] = float(liste[4])
59         ts[iligne] = float(liste[8])
60         fov[iligne] = float(liste[10])
61         zen[iligne] = float(liste[9])
62         for ki in range(0,4):
63             emis[ki,iligne] = float(liste[13+ki])
64             tb[ki,iligne] = float(liste[17+ki])
65             tdn[ki,iligne] = float(liste[21+ki])
66             tup[ki,iligne] = float(liste[25+ki])
67             trans[ki,iligne] = float(liste[29+ki])
68         iligne=iligne+1
69    fichier.close()
70    for i in range (0, nbtotal):
71        if (emis[2, i] < 0.):
72            emis[2, i] = NaN
73        if (emis[2, i] > 1.):
74            emis[2, i] = 0.99 
75    ############################
76    # definition of parameters #
77    ############################
78    bblat = nonzero(lat >= 65.) # only high latitudes
79    bbzen = nonzero(zen[bblat] < 20.) # only high latitudes
80    jour = jjr[bblat][bbzen]
81    llon = lon[bblat][bbzen] # longitude
82    llat = lat[bblat][bbzen] # latitude
83    zen_angle = zen[bblat][bbzen] # zenith angle
84    field_view = fov[bblat][bbzen] # field of view
85    e = emis[2, :][bblat][bbzen] # emissivity
86    t_surf = ts[bblat][bbzen] # surface temperature
87    t_up = tup[2, :][bblat][bbzen] # tup
88    t_dn = tdn[2, :][bblat][bbzen] # tdown
89    transm = trans[2, :][bblat][bbzen] # transmissivity
90    t_b = tb[2, :][bblat][bbzen] # brightness temperature
91    L = len(llon)
92    #######################################
93    # calculation of lamb_spec parameters #
94    #######################################
95    r = pi / 180.
96    opac = - cos(zen_angle * r) * np.log(transm) 
97    E3 = scipy.special.expn(3, opac)
98    theta_eff = np.zeros([L], float)
99    tdn_lamb = np.zeros([L], float)
100    tb_spec = np.zeros([L], float)
101    tb_lamb = np.zeros([L], float)
102    e_spec = np.zeros([L], float)
103    e_lamb = np.zeros([L], float)
104    e_spec_lamb = np.zeros([L], float)
105    print 'calculate lamb_spec'
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    # mixed spec and lamb assumptions #
116    ###################################
117    print 'calculate mixed emis'
118    s = np.array([0., 0.25, 0.5, 0.75, 1.]) # specularity coefficient (s=1 <=> specular)
119    S = len(s)
120    new_tdn = np.zeros([L, S], float)
121    new_emis = np.zeros([L, S], float)
122    for s_ind in range (0, S):
123        for ii in range (0, L):
124            new_tdn[ii, s_ind] = s[s_ind] * t_dn[ii] + (1. - s[s_ind]) * tdn_lamb[ii]
125            new_emis[ii, s_ind] = ((t_b[ii] - t_up[ii]) * exp(opac[ii]) - new_tdn[ii, s_ind]) / (t_surf[ii] - new_tdn[ii, s_ind])
126    ##############################################
127    # stack of mixed emis lamb_spec in .dat file #
128    ##############################################
129    #data_spec_lamb = open ('/mma/hermozol/Documents/Data_tests/monthly_GLACE/lamb_spec_param_near_nadir_SEPTEMBER2009.dat', 'a')
130    print 'stack data'
131    for ii in range (0, L):
132        data_spec_lamb.write(('%(jour)2.0f    %(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    %(e_mixed_00)5.4f    %(e_mixed_25)5.4f    %(e_mixed_50)5.4f    %(e_mixed_75)5.4f   %(e_mixed_100)5.4f \n' % {
133                                            'jour':jour[ii],
134                                            'lon':llon[ii],
135                                            'lat':llat[ii],
136                                            'za':zen_angle[ii],
137                                            'fov':field_view[ii],
138                                            'emis':e[ii],
139                                            'tsurf':t_surf[ii],
140                                            'tup':t_up[ii],
141                                            'tdn':t_dn[ii],
142                                            'trans':transm[ii],
143                                            'tb':t_b[ii],
144                                            'opac':opac[ii],
145                                            'theta_eff':theta_eff[ii],
146                                            'tdn_lamb':tdn_lamb[ii],
147                                            'tb_spec':tb_spec[ii],
148                                            'tb_lamb':tb_lamb[ii],
149                                            'e_spec':e_spec[ii],
150                                            'e_lamb':e_lamb[ii],
151                                            'e_diff':e_spec_lamb[ii],
152                                            'e_mixed_00':new_emis[ii, 0],
153                                            'e_mixed_25':new_emis[ii, 1],
154                                            'e_mixed_50':new_emis[ii, 2],
155                                            'e_mixed_75':new_emis[ii, 3],
156                                            'e_mixed_100':new_emis[ii, 4],
157                                            }))
158    data_spec_lamb.close()
159
Note: See TracBrowser for help on using the repository browser.