-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgsi_ncdiag.py
170 lines (142 loc) · 6.04 KB
/
gsi_ncdiag.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
__all__ = ['read_rad_ncdiag','read_cnv_ncdiag','read_rad_ncdiag0']
import numpy as np
import xarray as xa
def read_rad_ncdiag0(infile,**kwargs):
# Read in the old converted NetCDF diagnostic files
select_wavenumber = kwargs.get('chkwvn',None)
cal_ae = kwargs.get('cal_ae',False)
get_water_frac = kwargs.get('get_water_frac',False)
ds=xa.open_dataset(infile)
npts=ds.obsloc.size
nchs=ds.channel.size
rlat=ds.locinfo[0,:].values
rlon=ds.locinfo[1,:].values
rlon=(rlon+180)%360-180
qcflags=ds.qcflag.values.T
obs=ds.tb_obs.values.T
varinv=ds.errinv.values.T
omb_bc=ds.tbc.values.T
omb_nbc=ds.tbcnob.values.T
sim=obs-omb_nbc
clr=sim
if cal_ae:
aereff_fg=sim-clr
aereff_obs=obs-clr
aereff=0.5*abs(aereff_fg)+0.5*abs(aereff_obs)
if get_water_frac: water_frac=ds.locinfo[10,:].values
data_dict={'rlon':(['obsloc'],rlon),
'rlat':(['obsloc'],rlat),
'qcflag':(['obsloc','wavenumber'],qcflags),
'tb_obs':(['obsloc','wavenumber'],obs),
'tb_sim':(['obsloc','wavenumber'],sim),
'tb_clr':(['obsloc','wavenumber'],clr),
'varinv':(['obsloc','wavenumber'],varinv),
'omb_bc':(['obsloc','wavenumber'],omb_bc),
'omb_nbc':(['obsloc','wavenumber'],omb_nbc),
}
if cal_ae: data_dict['Ae']=(['obsloc','wavenumber'],aereff)
if get_water_frac: data_dict['water_frac']=(['obsloc'],water_frac)
coords_dict={'obsloc':np.arange(npts),'wavenumber':ds.wavenumber.values}
tmpds=xa.Dataset(data_dict,coords=coords_dict)
if type(select_wavenumber)==list or type(select_wavenumber)==float:
tmpds=tmpds.sel(wavenumber=select_wavenumber)
return tmpds
def read_rad_ncdiag(infile,**kwargs):
# Read in the new NetCDF diagnostic files
area = kwargs.get('area',None)
cornerll = kwargs.get('cornerll',None)
select_wavenumber = kwargs.get('chkwvn',None)
#select_qcflag = kwargs.get('sel_qc',-1)
cal_ae = kwargs.get('cal_ae',False)
get_water_frac = kwargs.get('get_water_frac',False)
if area != 'Glb':
minlat, maxlat, minlon, maxlon = cornerll
ds=xa.open_dataset(infile)
npts=int(ds.nobs.size/ds.nchans.size)
nchs=ds.nchans.size
rlat=np.reshape(ds.Latitude.values,(npts,nchs))[:,0]
rlon=np.reshape(ds.Longitude.values,(npts,nchs))[:,0]
rlon=(rlon+180)%360-180
qcflags=np.reshape(ds.QC_Flag.values,(npts,nchs))
obs=np.reshape(ds.Observation.values,(npts,nchs))
sim=np.reshape(ds.Simulated_Tb.values,(npts,nchs))
clr=np.reshape(ds.Clearsky_Tb.values,(npts,nchs))
varinv=np.reshape(ds.Inverse_Observation_Error.values,(npts,nchs))
omb_bc=np.reshape(ds.Obs_Minus_Forecast_adjusted.values,(npts,nchs))
omb_nbc=np.reshape(ds.Obs_Minus_Forecast_unadjusted.values,(npts,nchs))
mask = (~np.isnan(rlat))
if area != 'Glb':
mask = mask&((rlon<maxlon)&(rlon>minlon)&(rlat>minlat)&(rlat<maxlat))
if cal_ae:
aereff_fg=sim-clr
aereff_obs=obs-clr
aereff=0.5*abs(aereff_fg)+0.5*abs(aereff_obs)
if get_water_frac: water_frac=np.reshape(ds.Water_Fraction.values,(npts,nchs))
data_dict={'rlon':(['obsloc'],rlon[mask]),
'rlat':(['obsloc'],rlat[mask]),
'qcflag':(['obsloc','wavenumber'],qcflags[mask,:]),
'tb_obs':(['obsloc','wavenumber'],obs[mask,:]),
'tb_sim':(['obsloc','wavenumber'],sim[mask,:]),
'tb_clr':(['obsloc','wavenumber'],clr[mask,:]),
'varinv':(['obsloc','wavenumber'],varinv[mask,:]),
'omb_bc':(['obsloc','wavenumber'],omb_bc[mask,:]),
'omb_nbc':(['obsloc','wavenumber'],omb_nbc[mask,:]),
}
if cal_ae: data_dict['Ae']=(['obsloc','wavenumber'],aereff)
if get_water_frac: data_dict['water_frac']=(['obsloc','wavenumber'],water_frac)
out_npts=np.count_nonzero(mask)
coords_dict={'obsloc':np.arange(out_npts),'wavenumber':ds.wavenumber.values}
tmpds=xa.Dataset(data_dict,coords=coords_dict)
if select_wavenumber is not None:
if (type(select_wavenumber)==list or
type(select_wavenumber)==float or
select_wavenumber.size>0):
tmpds=tmpds.sel(wavenumber=select_wavenumber)
#if select_qcflag != -1:
# qcmask = (tmpds.qcflag==select_qcflag)
# tmpds = tmpds
return tmpds
def read_cnv_ncdiag(infile,**kwargs):
useqc = kwargs.get('useqc',None)
sel_bufr = kwargs.get('sel_bufr',None)
is_sfc = kwargs.get('is_sfc',None)
is_uv = kwargs.get('is_uv',None)
area = kwargs.get('area',None)
cornerll = kwargs.get('cornerll',None)
if area != 'Glb':
minlat, maxlat, minlon, maxlon = cornerll
ds = xa.open_dataset(infile)
rlat = ds.Latitude.data
rlon = ds.Longitude.data
rlon = (rlon+180)%360-180
sta_id = ds.Station_ID.data
obstype = ds.Observation_Type.data
sta_elev = ds.Station_Elevation.data
qcflag = ds.Analysis_Use_Flag.data
errinv = ds.Errinv_Final.data
obs = ds.Observation.data
omb_bc = ds.Obs_Minus_Forecast_adjusted.data
omb_nbc = ds.Obs_Minus_Forecast_unadjusted.data
mask = (~np.isnan(rlat))
if area != 'Glb':
mask = mask&((rlon<maxlon)&(rlon>minlon)&(rlat>minlat)&(rlat<maxlat))
if useqc:
mask = mask&(qcflag==1)
if sel_bufr != 'all':
mask = mask&(obstype==int(bufrtype))
npts = rlat[mask].size
data_dict = {'rlon':(['obsloc'],rlon[mask]),
'rlat':(['obsloc'],rlat[mask]),
'qcflag':(['obsloc'],qcflag[mask]),
'obs':(['obsloc'],obs[mask]),
'omb_bc':(['obsloc'],omb_bc[mask]),
'omb_nbc':(['obsloc'],omb_nbc[mask]),
'sta_id':(['obsloc'],sta_id[mask]),
'obstype':(['obsloc'],obstype[mask]),
}
if not is_sfc:
pres = ds.Pressure.data
data_dict['pres'] = (['obsloc'],pres[mask])
coords_dict = {'obsloc':np.arange(npts)}
tmpds = xa.Dataset(data_dict,coords=coords_dict)
return tmpds