"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