UFO
HydrometeorCheckATMS.cc
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2019 UCAR
3  *
4  * This software is licensed under the terms of the Apache Licence Version 2.0
5  * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
6  */
7 
8 #include "ufo/filters/obsfunctions/HydrometeorCheckATMS.h"
9 
10 #include <cmath>
11 
12 #include <algorithm>
13 #include <iomanip>
14 #include <iostream>
15 #include <set>
16 #include <string>
17 #include <vector>
18 
19 #include "ioda/ObsDataVector.h"
20 #include "oops/util/IntSetParser.h"
21 #include "ufo/filters/obsfunctions/CLWRetMW.h"
22 #include "ufo/filters/obsfunctions/ObsErrorModelRamp.h"
23 #include "ufo/filters/Variable.h"
24 #include "ufo/utils/Constants.h"
25 
26 
27 namespace ufo {
28 
29 static ObsFunctionMaker<HydrometeorCheckATMS> makerHydrometeorCheckATMS_("HydrometeorCheckATMS");
30 
31 // -----------------------------------------------------------------------------
32 
33 HydrometeorCheckATMS::HydrometeorCheckATMS(const eckit::LocalConfiguration & conf)
34  : invars_() {
35  // Initialize options
36  options_.deserialize(conf);
37 
38  // Get channels from options
39  std::set<int> channelset = oops::parseIntSet(options_.channelList);
40  std::copy(channelset.begin(), channelset.end(), std::back_inserter(channels_));
41  ASSERT(channels_.size() > 0);
42 
43  // Get test groups from options
44  const std::string &biastermgrp = options_.testBiasTerm.value();
45  const std::string &biasgrp = options_.testBias.value();
46  const std::string &hofxgrp = options_.testHofX.value();
47 
48  // Include required variables from ObsDiag
49  invars_ += Variable("brightness_temperature_jacobian_surface_emissivity@ObsDiag", channels_);
50  invars_ += Variable("brightness_temperature_assuming_clear_sky@ObsDiag", channels_);
51 
52  // Include list of required data from ObsSpace
53  invars_ += Variable("brightness_temperature@ObsValue", channels_);
54  invars_ += Variable("brightness_temperature@"+biasgrp, channels_);
55  invars_ += Variable("brightness_temperature@"+hofxgrp, channels_);
56  invars_ += Variable("constant@"+biastermgrp, channels_);
57  invars_ += Variable("scan_angle_order_4@"+biastermgrp, channels_);
58  invars_ += Variable("scan_angle_order_3@"+biastermgrp, channels_);
59  invars_ += Variable("scan_angle_order_2@"+biastermgrp, channels_);
60  invars_ += Variable("scan_angle@"+biastermgrp, channels_);
61 
62  // Include list of required data from GeoVaLs
63  invars_ += Variable("water_area_fraction@GeoVaLs");
64  invars_ += Variable("land_area_fraction@GeoVaLs");
65  invars_ += Variable("latitude@MetaData");
66  invars_ += Variable("longitude@MetaData");
67 
68  // Include list of required data from ObsFunction
69  const Variable &obserrfunc = options_.obserrFunction.value();
70  invars_ += obserrfunc;
71 
72  const Variable &clwretfunc = options_.clwretFunction.value();
73  invars_ += clwretfunc;
74 }
75 
76 // -----------------------------------------------------------------------------
77 
78 HydrometeorCheckATMS::~HydrometeorCheckATMS() {}
79 
80 // -----------------------------------------------------------------------------
81 
82 void HydrometeorCheckATMS::compute(const ObsFilterData & in,
83  ioda::ObsDataVector<float> & out) const {
84  // Get dimensions
85  size_t nlocs = in.nlocs();
86  size_t nchans = channels_.size();
87 
88  // Set channel index (AMSU-A like channels)
89  int ich238 = 0, ich314 = 1, ich503 = 2, ich528 = 4, ich536 = 5;
90  int ich544 = 6, ich549 = 7, ich890 = 15;
91 
92  // Set channel index (MHS like channels)
93  int ich1650 = 16, ich1830a = 17, ich1830b = 18, ich1830c = 19, ich1830d = 20, ich1830e = 21;
94 
95  // Get test groups from options
96  const std::string &biastermgrp = options_.testBiasTerm.value();
97  const std::string &biasgrp = options_.testBias.value();
98  const std::string &hofxgrp = options_.testHofX.value();
99 
100  // Get clear-sky observation error from options
101  const std::vector<float> &obserr0 = options_.obserrClearSky.value();
102 
103  // Get area fraction of each surface type
104  std::vector<float> water_frac(nlocs);
105  in.get(Variable("water_area_fraction@GeoVaLs"), water_frac);
106 
107  std::vector<float> land_frac(nlocs);
108  in.get(Variable("land_area_fraction@GeoVaLs"), land_frac);
109 
110  std::vector<float> lat(nlocs);
111  in.get(Variable("latitude@MetaData"), lat);
112 
113  std::vector<float> lon(nlocs);
114  in.get(Variable("longitude@MetaData"), lon);
115 
116  // Get surface temperature jacobian
117  std::vector<std::vector<float>> dbtde(nchans, std::vector<float>(nlocs));
118  for (size_t ichan = 0; ichan < nchans; ++ichan) {
119  in.get(Variable("brightness_temperature_jacobian_surface_emissivity@ObsDiag", channels_)[ichan],
120  dbtde[ichan]);
121  }
122 
123  // Get HofX for clear-sky simulation (53.6 GHz)
124  std::vector<float> hofxclr536(nlocs);
125  in.get(Variable("brightness_temperature_assuming_clear_sky@ObsDiag", channels_)[ich536],
126  hofxclr536);
127 
128  // Get HofX for clear-sky simulation (89 GHz)
129  std::vector<float> hofxclr890(nlocs);
130  in.get(Variable("brightness_temperature_assuming_clear_sky@ObsDiag", channels_)[ich890],
131  hofxclr890);
132 
133  // Get HofX for clear-sky simulation (165 GHz)
134  std::vector<float> hofxclr1650(nlocs);
135  in.get(Variable("brightness_temperature_assuming_clear_sky@ObsDiag", channels_)[ich1650],
136  hofxclr1650);
137 
138  // Get ObsBiasTerm: constant term for 23.8GHz channel
139  std::vector<float> bias_const238(nlocs);
140  in.get(Variable("constant@"+biastermgrp, channels_)[ich238], bias_const238);
141 
142  // Get ObsBiasTerm: scan angle terms for 23.8GHz channel
143  size_t nangs = 4;
144  std::vector<float> values(nlocs);
145  std::vector<std::string> scanterms(nangs);
146  std::vector<float> bias_scanang238(nlocs);
147  scanterms[0] = "scan_angle_order_4@"+biastermgrp;
148  scanterms[1] = "scan_angle_order_3@"+biastermgrp;
149  scanterms[2] = "scan_angle_order_2@"+biastermgrp;
150  scanterms[3] = "scan_angle@"+biastermgrp;
151  for (size_t iang = 0; iang < nangs; ++iang) {
152  in.get(Variable(scanterms[iang], channels_)[ich238], values);
153  for (size_t iloc = 0; iloc < nlocs; ++iloc) {
154  bias_scanang238[iloc] = bias_scanang238[iloc] + values[iloc];
155  }
156  }
157 
158  // Calculate bias-corrected innovation: Observation - HofX - bias
159  std::vector<std::vector<float>> btobs(nchans, std::vector<float>(nlocs));
160  std::vector<std::vector<float>> bias(nchans, std::vector<float>(nlocs));
161  std::vector<std::vector<float>> innov(nchans, std::vector<float>(nlocs));
162  std::vector<float> hofx(nlocs);
163  for (size_t ichan = 0; ichan < nchans; ++ichan) {
164  in.get(Variable("brightness_temperature@ObsValue", channels_)[ichan], btobs[ichan]);
165  in.get(Variable("brightness_temperature@"+biasgrp, channels_)[ichan], bias[ichan]);
166  in.get(Variable("brightness_temperature@"+hofxgrp, channels_)[ichan], hofx);
167  for (size_t iloc = 0; iloc < nlocs; ++iloc) {
168  innov[ichan][iloc] = btobs[ichan][iloc] - hofx[iloc] - bias[ichan][iloc];
169  }
170  }
171 
172  // Get all-sky observation error from ObsFunction
173  const Variable &obserrvar = options_.obserrFunction.value();
174  ioda::ObsDataVector<float> obserr(in.obsspace(), obserrvar.toOopsVariables());
175  in.get(obserrvar, obserr);
176 
177  // Get CLW retrieval based on observation from ObsFunction
178  const Variable &clwretvar = options_.clwretFunction.value();
179  ioda::ObsDataVector<float> clwobs(in.obsspace(), clwretvar.toOopsVariables());
180  in.get(clwretvar, clwobs);
181 
182  // Set parameters
183  float w1f6 = 1.0/10.0, w2f6 = 1.0/0.80;
184  float w1f4 = 1.0/0.30, w2f4 = 1.0/1.80;
185 
186  std::vector<std::vector<int>> affected_channels(nchans, std::vector<int>(nlocs));
187  // Loop over locations
188  // Combined cloud-precipitation-surface checks
189  for (size_t iloc = 0; iloc < nlocs; ++iloc) {
190  // Check surface type
191  bool luse = (water_frac[iloc] > 0.99 || land_frac[iloc] > 0.99);
192 
193  // Initialization
194  for (size_t ich = 0; ich < nchans; ++ich) {
195  affected_channels[ich][iloc] = 0;
196  }
197 
198  // Calculate cloud effect from 53.6 GHz (Channel 6)
199  float cldeff_obs536 = 0.0;
200  if (water_frac[iloc] > 0.99) {
201  cldeff_obs536 = btobs[ich536][iloc] - hofxclr536[iloc] - bias[ich536][iloc];
202  }
203 
204  // Calculate cloud effect from 89 GHz (Channel 16)
205  float cldeff_obs890 = 0.0;
206  if (water_frac[iloc] > 0.99) {
207  cldeff_obs890 = btobs[ich890][iloc] - hofxclr890[iloc] - bias[ich890][iloc];
208  }
209 
210  // Calculate cloud effect from 165 GHz (Channel 17)
211  float cldeff_obs1650 = 0.0;
212  if (water_frac[iloc] > 0.99) {
213  cldeff_obs1650 = btobs[ich1650][iloc] - hofxclr1650[iloc] - bias[ich1650][iloc];
214  }
215 
216  // Calculate scattering effect
217  std::vector<float> factch4(nlocs);
218  std::vector<float> factch6(nlocs);
219  float btobsbc238 = btobs[ich238][iloc] - bias_const238[iloc]
220  - bias_scanang238[iloc];
221  float clwx = 0.6;
222  float dsval = 0.8;
223  if (water_frac[iloc] > 0.99) {
224  clwx = 0.0;
225  dsval = ((2.410 - 0.0098 * btobsbc238) * innov[ich238][iloc] +
226  0.454 * innov[ich314][iloc] - innov[ich890][iloc]) * w1f6;
227  dsval = std::max(static_cast<float>(0.0), dsval);
228  }
229  factch4[iloc] = pow(clwx, 2) + pow(innov[ich528][iloc] * w2f4, 2);
230  factch6[iloc] = pow(dsval, 2) + pow(innov[ich544][iloc] * w2f6, 2);
231 
232  // Window channel sanity check
233  // If any of the window channels is bad, skip all window channels
234  // List of surface sensitivity channels
235  std::vector<float> OmFs{std::abs(innov[ich238][iloc]), std::abs(innov[ich314][iloc]),
236  std::abs(innov[ich503][iloc]), std::abs(innov[ich528][iloc]),
237  std::abs(innov[ich536][iloc]), std::abs(innov[ich544][iloc]),
238  std::abs(innov[ich890][iloc])};
239  bool result = false;
240  result = any_of(OmFs.begin(), OmFs.end(), [](float x){return x > 200.0;});
241 
242  if (result) {
243  // remove channels 1-7, 16, 17-22
244  for (size_t ich = ich238; ich <= ich544; ++ich) {
245  affected_channels[ich][iloc] = 1;
246  }
247  affected_channels[ich890][iloc] = 1;
248  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
249  affected_channels[ich][iloc] = 1;
250  }
251  } else {
252  // Hydrometeor check over water surface
253  if (water_frac[iloc] > 0.99) {
254  // Cloud water retrieval sanity check
255  if (clwobs[0][iloc] > 999.0) {
256  for (size_t ich = ich238; ich <= ich544; ++ich) {
257  affected_channels[ich][iloc] = 1;
258  }
259  affected_channels[ich890][iloc] = 1;
260  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
261  affected_channels[ich][iloc] = 1;
262  }
263  }
264  // Precipitation check (factch6: 54.4 GHz)
265  if (factch6[iloc] >= 1.0) {
266  for (size_t ich = ich238; ich <= ich544; ++ich) {
267  affected_channels[ich][iloc] = 1;
268  }
269  affected_channels[ich890][iloc] = 1;
270  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
271  affected_channels[ich][iloc] = 1;
272  }
273  // Scattering check (53.6GHz cloud effect)
274  } else if (cldeff_obs536 < -0.5) {
275  for (size_t ich = ich238; ich <= ich544; ++ich) {
276  affected_channels[ich][iloc] = 1;
277  }
278  affected_channels[ich890][iloc] = 1;
279  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
280  affected_channels[ich][iloc] = 1;
281  }
282  // Scattering checks (89GHz vs 166GHz)
283  } else if (std::abs(cldeff_obs890 - cldeff_obs1650) > 10.0) {
284  affected_channels[ich890][iloc] = 1;
285  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
286  affected_channels[ich][iloc] = 1;
287  }
288  if (std::abs(cldeff_obs890 - cldeff_obs1650) > 15.0) {
289  for (size_t ich = ich238; ich <= ich544; ++ich) {
290  affected_channels[ich][iloc] = 1;
291  }
292  }
293  // Sensitivity of BT to the surface emissivity check
294  } else {
295  float thrd238 = 0.025, thrd314 = 0.015, thrd503 = 0.030, thrd890 = 0.030;
296  float de238 = 0.0, de314 = 0.0, de503 = 0.0, de890 = 0.0;
297  float dbtde238 = dbtde[ich238][iloc];
298  float dbtde314 = dbtde[ich314][iloc];
299  float dbtde503 = dbtde[ich503][iloc];
300  float dbtde890 = dbtde[ich890][iloc];
301  if (dbtde238 != 0.0) de238 = std::abs(innov[ich238][iloc]) / dbtde238 *
302  (obserr0[ich238] / obserr[ich238][iloc]) *
303  (1.0 - std::max(1.0, 10.0*clwobs[0][iloc]));
304  if (dbtde314 != 0.0) de314 = std::abs(innov[ich314][iloc]) / dbtde314 *
305  (obserr0[ich314] / obserr[ich314][iloc]) *
306  (1.0 - std::max(1.0, 10.0*clwobs[0][iloc]));
307  if (dbtde503 != 0.0) de503 = std::abs(innov[ich503][iloc]) / dbtde503 *
308  (obserr0[ich503] / obserr[ich503][iloc]) *
309  (1.0 - std::max(1.0, 10.0*clwobs[0][iloc]));
310  if (dbtde890 != 0.0) de890 = std::abs(innov[ich890][iloc]) / dbtde890 *
311  (obserr0[ich890] / obserr[ich890][iloc]) *
312  (1.0 - std::max(1.0, 10.0*clwobs[0][iloc]));
313  bool qcemiss = false;
314  qcemiss = de238 > thrd238 || de314 > thrd314 || de503 > thrd503 || de890 > thrd890;
315  if (qcemiss) {
316  for (size_t ich = ich238; ich <= ich536; ++ich) {
317  affected_channels[ich][iloc] = 1;
318  }
319  affected_channels[ich890][iloc] = 1;
320  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
321  affected_channels[ich][iloc] = 1;
322  }
323  }
324  }
325  } else {
326  // Hydrometeor check over non-water (land/sea ice/snow) surface
327  // Precipitation check (factch6)
328  if (factch6[iloc] >= 1.0 || luse == false) {
329  for (size_t ich = ich238; ich <= ich544; ++ich) {
330  affected_channels[ich][iloc] = 1;
331  }
332  affected_channels[ich890][iloc] = 1;
333  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
334  affected_channels[ich][iloc] = 1;
335  }
336  // Thick cloud check (factch4)
337  } else if (factch4[iloc] > 0.5) {
338  for (size_t ich = ich238; ich <= ich536; ++ich) {
339  affected_channels[ich][iloc] = 1;
340  }
341  affected_channels[ich890][iloc] = 1;
342  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
343  affected_channels[ich][iloc] = 1;
344  }
345  // Sensitivity of BT to the surface emissivity check
346  } else {
347  float thrd238 = 0.020, thrd314 = 0.015, thrd503 = 0.035, thrd890 = 0.015;
348  float de238 = 0.0, de314 = 0.0, de503 = 0.0, de890 = 0.0;
349  float dbtde238 = dbtde[ich238][iloc];
350  float dbtde314 = dbtde[ich314][iloc];
351  float dbtde503 = dbtde[ich503][iloc];
352  float dbtde890 = dbtde[ich890][iloc];
353  if (dbtde238 != 0.0) de238 = std::abs(innov[ich238][iloc]) / dbtde238;
354  if (dbtde314 != 0.0) de314 = std::abs(innov[ich314][iloc]) / dbtde314;
355  if (dbtde503 != 0.0) de503 = std::abs(innov[ich503][iloc]) / dbtde503;
356  if (dbtde890 != 0.0) de890 = std::abs(innov[ich890][iloc]) / dbtde890;
357  bool qcemiss = false;
358  qcemiss = de238 > thrd238 || de314 > thrd314 || de503 > thrd503 || de890 > thrd890;
359  if (qcemiss) {
360  for (size_t ich = ich238; ich <= ich536; ++ich) {
361  affected_channels[ich][iloc] = 1;
362  }
363  affected_channels[ich890][iloc] = 1;
364  for (size_t ich = ich1650; ich <= ich1830e; ++ich) {
365  affected_channels[ich][iloc] = 1;
366  }
367  }
368  }
369  // surface type
370  }
371  // window channel sanity check
372  }
373  // loop over locations
374  }
375  // Output
376  for (size_t ichan = 0; ichan < nchans; ++ichan) {
377  for (size_t iloc = 0; iloc < nlocs; ++iloc) {
378  out[ichan][iloc] = affected_channels[ichan][iloc];
379  }
380  }
381 }
382 
383 // -----------------------------------------------------------------------------
384 
385 const ufo::Variables & HydrometeorCheckATMS::requiredVariables() const {
386  return invars_;
387 }
388 
389 // -----------------------------------------------------------------------------
390 
391 } // namespace ufo
HydrometeorCheckATMS.h
ufo::ObsFilterData::nlocs
size_t nlocs() const
Returns number of locations.
Definition: ObsFilterData.cc:66
ufo::Variables
Definition: src/ufo/filters/Variables.h:24
ufo::HydrometeorCheckATMSParameters::testBiasTerm
oops::Parameter< std::string > testBiasTerm
Definition: HydrometeorCheckATMS.h:60
CLWRetMW.h
ufo::HydrometeorCheckATMSParameters::channelList
oops::RequiredParameter< std::string > channelList
List of channels available for assimilation.
Definition: HydrometeorCheckATMS.h:36
ufo::ObsFilterData::obsspace
ioda::ObsSpace & obsspace() const
Returns reference to ObsSpace associated with ObsFilterData.
Definition: src/ufo/filters/ObsFilterData.h:84
ufo::Variable::toOopsVariables
oops::Variables toOopsVariables() const
Definition: Variable.cc:129
ufo
Definition: RunCRTM.h:27
ufo::ObsFunctionMaker
Definition: ObsFunctionBase.h:60
ufo::HydrometeorCheckATMS::compute
void compute(const ObsFilterData &, ioda::ObsDataVector< float > &) const
compute the result of the function
Definition: HydrometeorCheckATMS.cc:82
ufo::HydrometeorCheckATMS::invars_
ufo::Variables invars_
Definition: HydrometeorCheckATMS.h:89
ufo::HydrometeorCheckATMS::HydrometeorCheckATMS
HydrometeorCheckATMS(const eckit::LocalConfiguration &)
Definition: HydrometeorCheckATMS.cc:33
ufo::abs
util::Duration abs(const util::Duration &duration)
Definition: TrackCheckUtils.h:31
ioda::ObsDataVector
Definition: BackgroundCheck.h:26
ufo::HydrometeorCheckATMSParameters::obserrFunction
oops::RequiredParameter< Variable > obserrFunction
Definition: HydrometeorCheckATMS.h:43
ufo::HydrometeorCheckATMSParameters::testHofX
oops::Parameter< std::string > testHofX
Name of the HofX group used to replace the default group (default is HofX)
Definition: HydrometeorCheckATMS.h:49
ufo::HydrometeorCheckATMS::channels_
std::vector< int > channels_
Definition: HydrometeorCheckATMS.h:90
ufo::HydrometeorCheckATMS::~HydrometeorCheckATMS
~HydrometeorCheckATMS()
Definition: HydrometeorCheckATMS.cc:78
ufo::HydrometeorCheckATMSParameters::obserrClearSky
oops::RequiredParameter< std::vector< float > > obserrClearSky
Observation error for each channel under clear-sky condition.
Definition: HydrometeorCheckATMS.h:39
ufo::HydrometeorCheckATMSParameters::testBias
oops::Parameter< std::string > testBias
Definition: HydrometeorCheckATMS.h:54
ufo::HydrometeorCheckATMSParameters::clwretFunction
oops::RequiredParameter< Variable > clwretFunction
Function used to retrieve the cloud liquid water from observation (CLWRetMW)
Definition: HydrometeorCheckATMS.h:46
Constants.h
ufo::ObsFilterData::get
void get(const Variable &, std::vector< float > &) const
Gets requested data from ObsFilterData.
Definition: ObsFilterData.cc:130
ufo::HydrometeorCheckATMS::requiredVariables
const ufo::Variables & requiredVariables() const
geovals required to compute the function
Definition: HydrometeorCheckATMS.cc:385
ufo::Variable
Definition: Variable.h:23
ObsErrorModelRamp.h
ufo::makerHydrometeorCheckATMS_
static ObsFunctionMaker< HydrometeorCheckATMS > makerHydrometeorCheckATMS_("HydrometeorCheckATMS")
ufo::ObsFilterData
ObsFilterData provides access to all data related to an ObsFilter.
Definition: src/ufo/filters/ObsFilterData.h:40
ufo::HydrometeorCheckATMS::options_
HydrometeorCheckATMSParameters options_
Definition: HydrometeorCheckATMS.h:91
conf
Definition: conf.py:1
Variable.h