UFO
MWCLWCheck.cc
Go to the documentation of this file.
1 /*
2  * (C) Copyright 2018-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 
9 
10 #include <math.h>
11 
12 #include <algorithm>
13 #include <iostream>
14 #include <vector>
15 
16 #include "eckit/config/Configuration.h"
17 
18 #include "ioda/ObsDataVector.h"
19 #include "ioda/ObsSpace.h"
20 #include "oops/util/Logger.h"
21 #include "oops/util/missingValues.h"
22 #include "ufo/filters/QCflags.h"
23 
24 namespace ufo {
25 
26 // -----------------------------------------------------------------------------
27 
28 MWCLWCheck::MWCLWCheck(ioda::ObsSpace & obsdb, const Parameters_ & parameters,
29  std::shared_ptr<ioda::ObsDataVector<int> > flags,
30  std::shared_ptr<ioda::ObsDataVector<float> > obserr)
31  : FilterBase(obsdb, parameters, flags, obserr), parameters_(parameters) {
32  oops::Log::debug() << "MWCLWCheck: config = " << parameters_ << std::endl;
33 
34  const oops::Variables &invars = parameters_.clwVariables.value();
35  ASSERT(invars.size() == 2);
36  const Variable var0(invars[0] + "@HofX");
37  const Variable var1(invars[1] + "@HofX");
38  allvars_ += var0;
39  allvars_ += var1;
40 }
41 
42 // -----------------------------------------------------------------------------
43 
45 
46 // -----------------------------------------------------------------------------
47 
48 void MWCLWCheck::applyFilter(const std::vector<bool> & apply,
49  const Variables & filtervars,
50  std::vector<std::vector<bool>> & flagged) const {
51  oops::Log::trace() << "MWCLWCheck postFilter" << std::endl;
52 
53  const oops::Variables observed = obsdb_.obsvariables();
54  const float missing = util::missingValue(missing);
55  float amsua_clw(float, float, float);
56 
57 // Get config
58  const std::vector<float> &clw_thresholds = parameters_.clwThresholds;
59 // clw_option controls how the clw is calculated:
60 // 1) Use observed BTs
61 // 2) Use calculated BTs
62 // 3) Symmetric calculation
63  const int clw_option = parameters_.clwOption;
64  const oops::Variables &invars = parameters_.clwVariables.value();
65  oops::Log::debug() << "MWCLWCheck: config = " << parameters_ << std::endl;
66 
67 // Number of channels to be tested and number of thresholds needs to be the same
68  ASSERT(clw_thresholds.size() == filtervars.nvars());
69 // Thresholds need to have non-missing values
70  for (size_t jv = 0; jv < filtervars.nvars(); ++jv) {
71  ASSERT(clw_thresholds[jv] != missing);
72  }
73 
74  ioda::ObsDataVector<float> obs(obsdb_, filtervars.toOopsVariables(), "ObsValue");
75  ioda::ObsDataVector<float> obs_for_calc(obsdb_, invars, "ObsValue");
76  ioda::ObsDataVector<float> sza(obsdb_, "sensor_zenith_angle", "MetaData");
77  ioda::ObsDataVector<float> clw(obsdb_, "cloud_liquid_water", "Diagnostic", false);
78  ioda::ObsDataVector<float> clw_guess_out(obsdb_, "clws_guess", "Diagnostic", false);
79  ioda::ObsDataVector<float> clw_obs_out(obsdb_, "clw_obs", "Diagnostic", false);
80 
81 // H(x)
82  const Variable var0(invars[0] + "@HofX");
83  const Variable var1(invars[1] + "@HofX");
84  std::vector<float> hofx0;
85  data_.get(var0, hofx0);
86  std::vector<float> hofx1;
87  data_.get(var1, hofx1);
88 
89  auto clw_obs = [&](size_t jobs) { return amsua_clw(
90  obs_for_calc[0][jobs], obs_for_calc[1][jobs], sza[0][jobs]); };
91 
92  auto clw_guess = [&](size_t jobs) { return amsua_clw(
93  hofx0[jobs], hofx1[jobs], sza[0][jobs]); };
94 
95 // Loop over obs locations calculating CLW from observations
96  for (size_t jobs = 0; jobs < obs.nlocs(); ++jobs) {
97  switch (clw_option) {
98  case 1 :
99  clw[0][jobs] = clw_obs(jobs);
100  break;
101  case 2 :
102  clw[0][jobs] = clw_guess(jobs);
103  break;
104  case 3 :
105  clw_obs_out[0][jobs] = clw_obs(jobs);
106  clw_guess_out[0][jobs] = clw_guess(jobs);
107  if (clw_obs_out[0][jobs] == missing || clw_guess_out[0][jobs] == missing) {
108  clw[0][jobs] = missing;
109  } else {
110  clw[0][jobs] = (clw_obs_out[0][jobs] + clw_guess_out[0][jobs])/2.0;
111  }
112  break;
113  default:
114  oops::Log::error() << "Invalid value for clw_option:" << clw_option << std::endl;
115  ABORT("Invalid value for clw_option");
116  }
117 
118 // Apply CLW threshold to observations
119  for (size_t jv = 0; jv < filtervars.nvars(); ++jv) {
120  if (apply[jobs]) {
121  if (clw[0][jobs] == missing ||
122  clw[0][jobs] > clw_thresholds[jv]) flagged[jv][jobs] = true;
123  }
124  }
125  }
126  clw.save("Derived");
127  clw_obs_out.save("Derived");
128  clw_guess_out.save("Derived");
129 }
130 
131 // -----------------------------------------------------------------------------
132 
133 float amsua_clw(float tobs1, float tobs2, float sza) {
134  const float d1 = 0.754;
135  const float d2 = -2.265;
136  const float missing = util::missingValue(missing);
137 
138  float clw;
139 
140  if (tobs1 != missing && tobs2 != missing && sza != missing &&
141  tobs1 <= 284.0 && tobs2 <= 284.0 && tobs1 > 0.0 && tobs2 > 0.0) {
142  float cossza = cos(M_PI * sza/180.0);
143  float d0 = 8.240 - (2.622 - 1.846*cossza)*cossza;
144  clw = cossza*(d0 + d1*std::log(285.0-tobs1)) + d2*std::log(285.0-tobs2);
145  clw = std::max(0.0f, clw);
146  } else {
147  clw = missing;
148  }
149 
150  return clw;
151 }
152 
153 // -----------------------------------------------------------------------------
154 
155 void MWCLWCheck::print(std::ostream & os) const {
156  os << "MWCLWCheck: config = " << parameters_ << std::endl;
157 }
158 
159 // -----------------------------------------------------------------------------
160 
161 } // namespace ufo
Base class for UFO QC filters.
Definition: FilterBase.h:45
MWCLWCheck(ioda::ObsSpace &, const Parameters_ &, std::shared_ptr< ioda::ObsDataVector< int > >, std::shared_ptr< ioda::ObsDataVector< float > >)
Definition: MWCLWCheck.cc:28
void print(std::ostream &) const override
Definition: MWCLWCheck.cc:155
Parameters_ parameters_
Definition: MWCLWCheck.h:73
void applyFilter(const std::vector< bool > &, const Variables &, std::vector< std::vector< bool >> &) const override
Definition: MWCLWCheck.cc:48
Parameters controlling the operation of the MWCLWCheck filter.
Definition: MWCLWCheck.h:36
oops::RequiredParameter< std::vector< float > > clwThresholds
Definition: MWCLWCheck.h:42
oops::RequiredParameter< oops::Variables > clwVariables
Definition: MWCLWCheck.h:40
oops::RequiredParameter< int > clwOption
Definition: MWCLWCheck.h:49
void get(const Variable &varname, std::vector< float > &values) const
Fills a std::vector with values of the specified variable.
ufo::Variables allvars_
ioda::ObsSpace & obsdb_
size_t nvars() const
Return the number of constituent "primitive" (single-channel) variables.
Definition: Variables.cc:104
oops::Variables toOopsVariables() const
Definition: Variables.cc:156
constexpr int clw
Definition: QCflags.h:28
constexpr int missing
Definition: QCflags.h:20
Definition: RunCRTM.h:27
float amsua_clw(float tobs1, float tobs2, float sza)
Definition: MWCLWCheck.cc:133