"Fossies" - the Fresh Open Source Software Archive 
Member "fimex-1.4.1/src/CDMMerger.cc" (30 Oct 2019, 11285 Bytes) of package /linux/privat/fimex-1.4.1.tar.gz:
As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style:
standard) with prefixed line numbers and
code folding option.
Alternatively you can here
view or
download the uninterpreted source code file.
For more information about "CDMMerger.cc" see the
Fossies "Dox" file reference documentation.
1 /*
2 * Fimex, CDMMerger.cc
3 *
4 * (C) Copyright 2012-2013, met.no
5 *
6 * Project Info: https://wiki.met.no/fimex/start
7 *
8 * This library is free software; you can redistribute it and/or modify it
9 * under the terms of the GNU Lesser General Public License as published by
10 * the Free Software Foundation; either version 2.1 of the License, or
11 * (at your option) any later version.
12 *
13 * This library is distributed in the hope that it will be useful, but
14 * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with this library; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301,
21 * USA.
22 *
23 * Created on: Aug 28, 2012
24 * Author: Alexander Bürger
25 */
26
27 #include "fimex/CDMMerger.h"
28
29 #include "fimex/CDM.h"
30 #include "fimex/CDMBorderSmoothing_Linear.h"
31 #include "fimex/CDMException.h"
32 #include "fimex/CDMInterpolator.h"
33 #include "fimex/CDMOverlay.h"
34 #include "fimex/CDMconstants.h"
35 #include "fimex/Data.h"
36 #include "fimex/DataIndex.h"
37 #include "fimex/Logger.h"
38 #include "fimex/SpatialAxisSpec.h"
39 #include "fimex/coordSys/CoordinateAxis.h"
40
41 #include "CDMMergeUtils.h"
42
43 // clang-format off
44 #define THROW(x) do { std::ostringstream t; t << x; throw CDMException(t.str()); } while(false)
45 // clang-format on
46
47 using namespace MetNoFimex;
48 using namespace std;
49
50 // ========================================================================
51
52 namespace MetNoFimex {
53
54 static Logger_p logger(getLogger("fimex.CDMMerger"));
55
56 // ========================================================================
57
58 struct CDMMergerPrivate {
59 CDMReader_p readerI;
60 CDMReader_p readerO;
61
62 CDMInterpolator_p interpolatedOT; //! outer interpolated to target grid
63 CDMBorderSmoothing_p readerSmooth; //! smoothed (inner+outer) on inner grid
64 CDMInterpolator_p interpolatedST; //! smoothed (inner+outer), interpolated to target grid
65 CDMOverlay_p readerOverlay; //! final result, we forward getDataSlice there
66
67 int gridInterpolationMethod;
68 bool useOuterIfInnerUndefined;
69 bool keepOuterVariables;
70 CDMBorderSmoothing::SmoothingFactory_p smoothingFactory;
71
72 CDM makeCDM(const string& tproj, const values_v& tx, const values_v& ty,
73 const std::string& tx_unit, const std::string& ty_unit,
74 const CDMDataType& tx_type, const CDMDataType& ty_type);
75 values_v extendInnerAxis(CoordinateAxis_cp axisI, CoordinateAxis_cp axisO, const char* unit);
76 };
77
78 // ========================================================================
79
80 CDMMerger::CDMMerger(CDMReader_p inner, CDMReader_p outer)
81 : p(new CDMMergerPrivate)
82 {
83 p->readerI = inner;
84 p->readerO = outer;
85 p->gridInterpolationMethod = MIFI_INTERPOL_BILINEAR;
86 p->smoothingFactory = CDMBorderSmoothing::SmoothingFactory_p(new CDMBorderSmoothing_LinearFactory());
87 p->useOuterIfInnerUndefined = true;
88 p->keepOuterVariables = false;
89 }
90
91 CDMMerger::~CDMMerger() {}
92
93 // ------------------------------------------------------------------------
94
95 void CDMMerger::setSmoothing(CDMBorderSmoothing::SmoothingFactory_p smoothingFactory)
96 {
97 p->smoothingFactory = smoothingFactory;
98 if (p->readerSmooth.get() != 0)
99 p->readerSmooth->setSmoothing(smoothingFactory);
100 }
101
102 // ------------------------------------------------------------------------
103
104 void CDMMerger::setUseOuterIfInnerUndefined(bool useOuter)
105 {
106 p->useOuterIfInnerUndefined = useOuter;
107 if (p->readerSmooth.get() != 0)
108 p->readerSmooth->setUseOuterIfInnerUndefined(useOuter);
109 }
110 // ------------------------------------------------------------------------
111
112 void CDMMerger::setKeepOuterVariables(bool keepOuterVariables)
113 {
114 if (p->readerOverlay.get() != 0)
115 LOG4FIMEX(logger, Logger::WARN, "setting keepOuterVariables after setting target grid has"
116 " no effect on the already defined target grid");
117 p->keepOuterVariables = keepOuterVariables;
118 }
119 // ------------------------------------------------------------------------
120
121 void CDMMerger::setGridInterpolationMethod(int method)
122 {
123 if (p->readerOverlay.get() != 0)
124 LOG4FIMEX(logger, Logger::WARN, "setting grid interpolation method after setting target grid has"
125 " no effect on the already defined target grid");
126 p->gridInterpolationMethod = method;
127 }
128
129 // ------------------------------------------------------------------------
130
131 void CDMMerger::setTargetGrid(const string& proj, const string& tx_axis, const string& ty_axis,
132 const string& tx_unit, const string& ty_unit, const string& tx_type, const string& ty_type)
133 {
134 SpatialAxisSpec xAxisSpec(tx_axis);
135 SpatialAxisSpec yAxisSpec(ty_axis);
136 if (xAxisSpec.requireStartEnd() or yAxisSpec.requireStartEnd())
137 throw CDMException("setTargetGrid without axis start/end not implemented");
138
139 const CDMDataType xType = string2datatype(tx_type),
140 yType = string2datatype(ty_type);
141 const vector<double> tx_values = xAxisSpec.getAxisSteps(),
142 ty_values = yAxisSpec.getAxisSteps();
143 setTargetGrid(proj, tx_values, ty_values, tx_unit, ty_unit, xType, yType);
144 }
145
146 // ------------------------------------------------------------------------
147
148 void CDMMerger::setTargetGrid(const string& proj, const values_v& tx, const values_v& ty,
149 const std::string& tx_unit, const std::string& ty_unit,
150 const CDMDataType& tx_type, const CDMDataType& ty_type)
151 {
152 *cdm_ = p->makeCDM(proj, tx, ty, tx_unit, ty_unit, tx_type, ty_type);
153 }
154
155 // ------------------------------------------------------------------------
156
157 void CDMMerger::setTargetGridFromInner()
158 {
159 const CDM &cdmI = p->readerI->getCDM(), &cdmO = p->readerO->getCDM();
160
161 const CoordinateSystem_cp_v allCsI = listCoordinateSystems(p->readerI), allCsO = listCoordinateSystems(p->readerO);
162
163 const CDM::VarVec& varsI = cdmI.getVariables();
164 for (const CDMVariable& varI : varsI) {
165 const string& varName = varI.getName();
166 if (not cdmO.hasVariable(varName))
167 continue;
168 if (cdmI.hasDimension(varName) or cdmO.hasDimension(varName))
169 continue;
170
171 const CoordinateSystem_cp csI = findCompleteCoordinateSystemFor(allCsI, varName), csO = findCompleteCoordinateSystemFor(allCsO, varName);
172 if (!csI.get() || !csO.get())
173 continue;
174
175 if (not (csI->isSimpleSpatialGridded() and csO->isSimpleSpatialGridded()))
176 continue;
177 if (not (csI->hasProjection() and csO->hasProjection()))
178 continue;
179
180 Projection_cp projI = csI->getProjection(), projO = csO->getProjection();
181 if (projI->isDegree() != projO->isDegree())
182 continue;
183
184 CoordinateAxis_cp xAxisI = csI->getGeoXAxis(), yAxisI = csI->getGeoYAxis();
185
186 const char* unit = projI->isDegree() ? "degree" : "m";
187 const values_v vx = p->extendInnerAxis(xAxisI, csO->getGeoXAxis(), unit);
188 const values_v vy = p->extendInnerAxis(yAxisI, csO->getGeoYAxis(), unit);
189
190 LOG4FIMEX(logger, Logger::INFO, "extending grid for inner variable '" << varName << "'");
191 setTargetGrid(projI->getProj4String(), vx, vy,
192 cdmI.getUnits(xAxisI->getName()), cdmI.getUnits(yAxisI->getName()),
193 xAxisI->getDataType(), yAxisI->getDataType());
194 return;
195 }
196
197 LOG4FIMEX(logger, Logger::WARN, "extending grid failed, no inner variable with CS found");
198 }
199
200 // ------------------------------------------------------------------------
201
202 DataPtr CDMMerger::getDataSlice(const std::string &varName, size_t unLimDimPos)
203 {
204 if (not p->readerOverlay)
205 THROW("must call setTargetGrid or setTargetGridFromInner before getDataSlice");
206
207 return p->readerOverlay->getDataSlice(varName, unLimDimPos);
208 }
209
210 // ========================================================================
211
212 CDM CDMMergerPrivate::makeCDM(const string& proj, const values_v& tx, const values_v& ty,
213 const std::string& tx_unit, const std::string& ty_unit,
214 const CDMDataType& tx_type, const CDMDataType& ty_type)
215 {
216 readerSmooth = std::make_shared<CDMBorderSmoothing>(readerI, readerO);
217 readerSmooth->setSmoothing(smoothingFactory);
218 readerSmooth->setUseOuterIfInnerUndefined(useOuterIfInnerUndefined);
219
220 interpolatedST = std::make_shared<CDMInterpolator>(readerSmooth);
221 interpolatedST->changeProjection(gridInterpolationMethod, proj,
222 tx, ty, tx_unit, ty_unit, tx_type, ty_type);
223
224 readerOverlay = std::make_shared<CDMOverlay>(readerO, interpolatedST, gridInterpolationMethod, keepOuterVariables);
225
226 return readerOverlay->getCDM();
227 }
228
229 // ------------------------------------------------------------------------
230
231 values_v CDMMergerPrivate::extendInnerAxis(CoordinateAxis_cp axisI, CoordinateAxis_cp axisO, const char* unit)
232 {
233 const string &nameI = axisI->getName(), &nameO = axisO->getName();
234
235 const CDM &cdmI = readerI->getCDM(), &cdmO = readerO->getCDM();
236 if (not (cdmI.hasDimension(nameI) and cdmO.hasDimension(nameO)))
237 THROW("found an axis that is not a dimension, cannot merge");
238
239 const CDMDimension& dimI = cdmI.getDimension(nameI), &dimO = cdmO.getDimension(nameO);
240 if( dimI.isUnlimited() != dimO.isUnlimited() )
241 THROW("axis '" << nameI << "' / '" << nameO << "' is not unlimited in both inner and outer, cannot merge");
242
243 DataPtr axisDataI = readerI->getScaledDataInUnit(nameI, unit), axisDataO = readerO->getScaledDataInUnit(nameO, unit);
244
245 if( !axisDataI || axisDataI->size() < 2 )
246 THROW("no data for axis '" << nameI << "' in inner");
247 if( !axisDataO || axisDataO->size() < 2 )
248 THROW("no data for axis '" << nameO << "' in outer");
249
250 shared_array<double> valuesI = axisDataI->asDouble(), valuesO = axisDataO->asDouble();
251 const double stepI = valuesI[1] - valuesI[0], stepO = valuesO[1] - valuesO[0];
252
253 for(size_t i=2; i<axisDataI->size(); ++i) {
254 if( !equal(valuesI[i] - valuesI[i-1], stepI) )
255 THROW("axis '" << nameI << "' in inner does not have constant step size, cannot merge");
256 }
257 for(size_t i=2; i<axisDataO->size(); ++i) {
258 if( !equal(valuesO[i] - valuesO[i-1], stepO) )
259 THROW("axis '" << nameO << "' in outer does not have constant step size, cannot merge");
260 }
261
262 const double minI = stepI > 0 ? valuesI[0] : valuesI[axisDataI->size()-1], minO = stepO > 0 ? valuesO[0] : valuesO[axisDataO->size()-1];
263 const double maxI = stepI < 0 ? valuesI[0] : valuesI[axisDataI->size()-1], maxO = stepO < 0 ? valuesO[0] : valuesO[axisDataO->size()-1];
264 if( minI < minO || maxI > maxO )
265 THROW("top not inside bottom");
266
267 vector<double> reverse;
268 for(double nO = valuesI[0] - stepI; nO >= minO && nO <= maxO; nO -= stepI)
269 reverse.push_back(nO);
270 vector<double> extended(reverse.rbegin(), reverse.rend());
271 extended.insert(extended.end(), &valuesI[0], &valuesI[axisDataI->size()]);
272 for(double nO = valuesI[axisDataI->size()-1] + stepI; nO >= minO && nO <= maxO; nO += stepI)
273 extended.push_back(nO);
274 return extended;
275 }
276
277 } // namespace MetNoFimex