-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathread_trmm.py
171 lines (116 loc) · 4.46 KB
/
read_trmm.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
'''
Read TRMM hdf4 format
Raul Valenzuela
December, 2015
'''
from pyhdf.SD import SD, SDC
from datetime import datetime
from brightness_temperature import BT
import numpy as np
def retrieve_1B01(*arg):
hdf = SD(arg[0], SDC.READ)
Lon, Lat, date_beg, date_end, st, en = retrieve_ancillary(hdf)
' calculate brightness temperature from IR channel'
'**************************************************'
chIR = hdf.select('channels')[st:en, :, 4] # [mW cm^-2 micrometer^-1 sr^-1]
I = chIR * 1e-3 * (100 ** 2) * 1e6 # [J s^-1 m^-3]
wavelength = 12e-6 # [m]
T = BT(I, wavelength)
return Lon[st:en, :], Lat[st:en, :], [date_beg, date_end], T
def retrieve_1C21(*arg):
hdf = SD(arg[0], SDC.READ)
Lon, Lat, date_beg, date_end, st, en = retrieve_ancillary(hdf)
' osRain includes rays #20 to #30'
# data=hdf.select('osRain')[st:en,:,0] # [dBZ]
' normalSample correspond to dBZ in 140 levels'
data = hdf.select('normalSample')[st:en, :, 100] # [dBZ*100]
data = data.astype(float)
data[data <= -32700.] = np.nan
data = data / 100. # [dBZ]
# foo = hdf.select('normalSample')[:, :, :]
# print foo.shape
# return Lon[:,20:31], Lat[:,20:31], [date_beg,date_end], osRain
data = np.ma.masked_array(data, np.isnan(data))
return Lon[st:en, :], Lat[st:en, :], [date_beg, date_end], data
def retrieve_2A12(*arg):
hdf = SD(arg[0], SDC.READ)
Lon, Lat, date_beg, date_end, st, en = retrieve_ancillary(hdf)
# data=hdf.select('windSpeed')[st:en,:] # [m/s]
# data[data==-9999.9]=np.nan
# data=hdf.select('surfaceRain')[st:en,:] # [mm/hr]
# data[data==-9999.9]=np.nan
# data=hdf.select('surfacePrecipitation')[st:en,:] # [mm/hr]
# data[data==-9999.9]=np.nan
# data=hdf.select('freezingHeight')[st:en,:] # [m]
# data[data==-9999.9]=np.nan
# data=hdf.select('cloudWaterPath')[st:en,:] # [kg/m2]
# data[data==-9999.9]=np.nan
# data=hdf.select('rainWaterPath')[st:en,:] # [kg/m2]
# data[data==-9999.9]=np.nan
# data=hdf.select('iceWaterPath')[st:en,:] # [kg/m2]
# data[data==-9999.9]=np.nan
# data=hdf.select('totalPrecipitableWater')[st:en,:] # [mm]
# data[data==-9999.9]=np.nan
return Lon, Lat, [date_beg, date_end], data
def retrieve_2A25(*arg):
hdf = SD(arg[0], SDC.READ)
if arg[1] == 'rainrate':
field = 'nearSurfRain' # [mm/hr]
elif arg[1] == 'dBZnearSurf':
field = 'nearSurfZ' # [dBZ]
elif arg[1] == 'dBZlevel':
field = 'correctZFactor' # [dBZ]
Lon, Lat, date_beg, date_end, st, en = retrieve_ancillary(hdf)
data=hdf.select(field)[st:en,:] # [mm/hr]
data[data==-99.99]=np.nan
# data=hdf.select('zeta')[st:en,:,1] # [mm/hr]
# data[data==-99.99]=np.nan
# data=hdf.select('correctZFactor')[st:en,:,10] # [dBZ]
# data=hdf.select('rainFlag')[st:en]
# data=hdf.select('zmmax')[st:en,:] # [dBZ]
data = np.ma.masked_array(data, np.isnan(data))
return Lon[st:en, :], Lat[st:en, :], [date_beg, date_end], data
def print_dataset_1B01(*arg):
# FILE_NAME = arg[0] + '1B01.' + arg[1]
hdf = SD(arg[0], SDC.READ)
'List available SDS datasets'
for ds in hdf.datasets():
print ds
def print_dataset_1C21(*arg):
# FILE_NAME = arg[0] + '1C21.' + arg[1]
# print FILE_NAME
hdf = SD(arg[0], SDC.READ)
'List available SDS datasets'
for ds in hdf.datasets():
print ds
def print_dataset_2A12(*arg):
FILE_NAME = arg[0] + '1B01.' + arg[1]
hdf = SD(FILE_NAME, SDC.READ)
'List available SDS datasets'
for ds in hdf.datasets():
print ds
def print_dataset_2A25(*arg):
FILE_NAME = arg[0] + '1C21.' + arg[1]
print FILE_NAME
hdf = SD(FILE_NAME, SDC.READ)
'List available SDS datasets'
for ds in hdf.datasets():
print ds
def retrieve_ancillary(hdf):
Latitude = hdf.select('Latitude')[:]
Longitude = hdf.select('Longitude')[:]
' find indices of lon corresponding to Chile'
lon1D = np.amax(Longitude, axis=1)
idxsub = np.where((lon1D >= -90) & (lon1D <= -60))
st = idxsub[0][0]
en = idxsub[0][-1]
' get datetime of swath'
Yr = hdf.select('Year')[0]
Mo = hdf.select('Month')[0]
Dy = hdf.select('DayOfMonth')[0]
Hr = hdf.select('Hour')[st:en]
Mn = hdf.select('Minute')[st:en]
date_beg = datetime(Yr, Mo, Dy, Hr[0], Mn[0], 0)
date_end = datetime(Yr, Mo, Dy, Hr[-1], Mn[-1], 0)
return Longitude, Latitude, date_beg, date_end, st, en