"Fossies" - the Fresh Open Source Software Archive

Member "fimex-1.4.2/src/metgm/MetGmGroup5Ptr.cc" (6 Jan 2020, 25381 Bytes) of package /linux/privat/fimex-1.4.2.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 "MetGmGroup5Ptr.cc" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.4.1_vs_1.4.2.

    1 /*
    2  * Fimex
    3  *
    4  * (C) Copyright 2011-2019, 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 
   24 // internals
   25 //
   26 #include "MetGmUtils.h"
   27 #include "MetGmGroup5Ptr.h"
   28 #include "MetGmFileHandlePtr.h"
   29 #include "MetGmHandlePtr.h"
   30 #include "MetGmDimensionsTag.h"
   31 #include "MetGmVersion.h"
   32 
   33 // fimex
   34 //
   35 #include "fimex/Type2String.h"
   36 #include "fimex/Units.h"
   37 
   38 // standard
   39 #include <cassert>
   40 #include <cmath>
   41 
   42 namespace MetNoFimex {
   43 
   44 MetGmGroup5Ptr::MetGmGroup5Ptr(const std::shared_ptr<MetGmGroup3Ptr> gp3, const std::shared_ptr<MetGmHDTag> hdTag, const shared_array<float> data,
   45                                const std::string fillValue)
   46     : pGp3_(gp3)
   47     , hdTag_(hdTag)
   48     , data_(data)
   49     , fillValue_(fillValue)
   50 {
   51     }
   52 
   53 
   54 /*
   55  * transforming ND index to 1D index:
   56  *
   57  * (x, y)         -> x + (nx * y)
   58  * (y, x)         -> y + (ny * x)
   59  *
   60  * (x, y, z)      -> x + (nx * y) + (nx * ny) * z
   61  * (z, x, y)      -> z + (nz * x) + (nz * nx) * y
   62  *
   63  * (x, y, z, t)   -> x + (nx * y) + (nx * ny) * z + (nx * ny * nz) * t
   64  * (z, x, y, t)   -> z + (nz * x) + (nz * nx) * y + (nx * ny * nz) * t
   65  */
   66 
   67 
   68 /*
   69  * (z, x, y, t) = (x, y, z, t)
   70  * (z, x, y, t)   -> z + nz * x + Nzx * y + Nxyz * t -> z + zLeap + zxLeap + xyzLeap
   71  * (x, y, z, t)   -> x + nx * y + Nxy * z + Nxyz * t -> x + xLeap + xyLeap + xyzLeap
   72  *
   73  *
   74  *     reindex data - all slices
   75  *     Fimex               METGM
   76  * (x, y, z, slice) -> (z, x, y, slice)
   77  */
   78 void MetGmGroup5Ptr::toMetGmLayout()
   79 {
   80     if(hdTag_->asShort() !=  MetGmHDTag::HD_3D_T)
   81         return;
   82 
   83     shared_array<float> dataT(new float[hdTag_->totalSize()]);
   84 
   85     float* pos = data_.get();
   86     float* posT = dataT.get();
   87 
   88     size_t nz = hdTag_->zSize();
   89     size_t ny = hdTag_->ySize();
   90     size_t nx = hdTag_->xSize();
   91     size_t nt = hdTag_->tSize();
   92 
   93     size_t Nxy  = nx * ny;
   94     size_t Nzx  = nz * nx;
   95     size_t Nxyz = Nxy * nz;
   96 
   97 
   98     size_t xyzLeap = -Nxyz;
   99 
  100     for(size_t t = 0; t < nt; ++t) {
  101 
  102         xyzLeap += Nxyz;
  103 
  104         size_t xLeap  = -nx;
  105         size_t zxLeap = -Nzx;
  106 
  107         for(size_t y = 0; y < ny; ++y) {
  108 
  109             xLeap += nx;
  110             zxLeap += Nzx;
  111 
  112             size_t zLeapF = -nz;
  113 
  114             for(size_t x = 0; x < nx; ++x) {
  115 
  116                 zLeapF += nz;
  117 
  118                 size_t xyLeapF = -Nxy;
  119                 size_t xyLeapB = nz * Nxy;
  120                 // from front and back
  121                 //   at the same time
  122                 // to meet in the middle
  123                 size_t midz = nz / 2 + nz % 2;
  124                 for(size_t z = 0; z < midz; ++z) {
  125 
  126                     xyLeapF += Nxy;
  127                     xyLeapB -= Nxy;
  128 
  129                     float valueFz = *(pos + x + xLeap + xyLeapF + xyzLeap);
  130                     *(posT + z + zLeapF + zxLeap + xyzLeap) = valueFz;
  131                     float valueBz = *(pos + x + xLeap + xyLeapB + xyzLeap);
  132                     *(posT + (nz - 1 - z) + zLeapF + zxLeap + xyzLeap) = valueBz;
  133 
  134                 } // z
  135             } // x
  136         } // y
  137     } // slice
  138 
  139     data_ = dataT;
  140 }
  141 
  142 /*
  143  *     reindex data - all slices
  144  *     METGM               Fimex
  145  * (z, x, y, slice) -> (x, y, z, slice)
  146  */
  147 void MetGmGroup5Ptr::toFimexLayout()
  148 {
  149     if(hdTag_->asShort() !=  MetGmHDTag::HD_3D_T)
  150         return;
  151 
  152     shared_array<float> dataT(new float[hdTag_->totalSize()]);
  153 
  154     float* pos = data_.get();
  155     float* posT = dataT.get();
  156 
  157     size_t nz = hdTag_->zSize();
  158     size_t ny = hdTag_->ySize();
  159     size_t nx = hdTag_->xSize();
  160     size_t nt = hdTag_->tSize();
  161 
  162     size_t Nxy  = nx * ny;
  163     size_t Nzx  = nz * nx;
  164     size_t Nxyz = Nxy * nz;
  165 
  166 
  167     size_t xyzLeap = -Nxyz;
  168 
  169     for(size_t t = 0; t < nt; ++t) {
  170 
  171         xyzLeap += Nxyz;
  172 
  173         size_t xLeap  = -nx;
  174         size_t zxLeap = -Nzx;
  175 
  176         for(size_t y = 0; y < ny; ++y) {
  177 
  178             xLeap += nx;
  179             zxLeap += Nzx;
  180 
  181             size_t xyLeap = -Nxy;
  182 
  183             for(size_t z = 0; z < nz; ++z) {
  184 
  185                 xyLeap += Nxy;
  186 
  187                 size_t zLeapF = -nz;
  188                 size_t zLeapB = nx * nz;
  189 
  190                 // from front and back
  191                 //   at the same time
  192                 // to meet in the middle
  193                 size_t midx = nx / 2 + nx % 2;
  194                 for(size_t x = 0; x < midx; ++x) {
  195 
  196                     zLeapF += nz;
  197                     zLeapB -= nz;
  198 
  199                     float valueFx = *(pos + z + zLeapF + zxLeap + xyzLeap);
  200                     *(posT + x + xLeap + xyLeap + xyzLeap) = valueFx;
  201                     float valueBx = *(pos + z + zLeapB + zxLeap + xyzLeap);
  202                     *(posT + (nx - 1 - x) + xLeap + xyLeap + xyzLeap) = valueBx;
  203                 } //x
  204             } // z
  205         } // y
  206     } // slice
  207 
  208     data_ = dataT;
  209 }
  210 
  211 std::shared_ptr<MetGmGroup5Ptr> MetGmGroup5Ptr::createMetGmGroup5PtrForWriting(const CDMReader_p pCdmReader, const CDMVariable* pVariable,
  212                                                                                const std::shared_ptr<MetGmGroup3Ptr> pg3)
  213 {
  214     std::shared_ptr<MetGmHDTag> hdtag = MetGmHDTag::createMetGmDimensionsTagForWriting(pCdmReader, pVariable);
  215 
  216     switch(hdtag->asShort()) {
  217         case MetGmHDTag::HD_2D:
  218         case MetGmHDTag::HD_3D_T:
  219             {
  220                 std::string mgmUnits;
  221 
  222                 if(pg3->p_id() == 7) {
  223                     /**
  224                       * METGM is tricky here as it can accomodate both m and hPa units
  225                       * depending if it is geopotential_height or pressure but the
  226                       * mgm_get_param_unit is, for some reason, always returning 'm'
  227                       *
  228                       * it should honor the pr settings
  229                       */
  230                     if(hdtag->zTag()->pr() == 2) {
  231                         /**
  232                           * If the value specified for pr is 2,
  233                           * and the values for pid=7 (if reported)
  234                           * are heights given in meters above MSL
  235                           */
  236                         mgmUnits = "m";
  237                     } else {
  238                         mgmUnits = "hPa";
  239                     }
  240                 } else {
  241                     mgmUnits = std::string(mgm_get_param_unit(pg3->p_id(), *(pg3->mgmHandle())));
  242                 }
  243 
  244                 const CDM& cdmRef = pCdmReader->getCDM();
  245                 const std::string varName = pVariable->getName();
  246 
  247                 Units unitsConvertor;
  248                 if(! unitsConvertor.areConvertible(cdmRef.getUnits(varName), mgmUnits)) {
  249                     std::string msg(" can't convert from ");
  250                     msg.append(cdmRef.getUnits(varName)).append(" to ").append(mgmUnits).append(" for variable ").append(varName).append(" -- excluding");
  251                     MGM_MESSAGE_POINT(msg)
  252                     shared_array<float> empty;
  253                     return std::shared_ptr<MetGmGroup5Ptr>(new MetGmGroup5Ptr(pg3, hdtag, empty));
  254                 }
  255 
  256                 DataPtr raw_data  = pCdmReader->getScaledDataInUnit(varName, mgmUnits);
  257 
  258                 std::shared_ptr<MetGmGroup5Ptr> gp5(new MetGmGroup5Ptr(pg3, hdtag, raw_data->asFloat()));
  259 
  260                 /*
  261                  *     reindex data - all slices
  262                  *     Fimex               METGM
  263                  * (x, y, z, slice) -> (z, x, y, slice)
  264                  */
  265                 gp5->toMetGmLayout();
  266 
  267                 return gp5;
  268             }
  269             break;
  270         case MetGmHDTag::HD_0D:
  271         case MetGmHDTag::HD_0D_T:
  272         case MetGmHDTag::HD_1D:
  273         case MetGmHDTag::HD_1D_T:
  274         case MetGmHDTag::HD_2D_T:
  275         case MetGmHDTag::HD_3D:
  276         default:
  277         throw CDMException(std::string(__FUNCTION__) + std::string(": dimensionality not supported yet :") + hdtag->asString() + " for " + pVariable->getName());
  278     }
  279 }
  280 
  281 std::shared_ptr<MetGmGroup5Ptr> MetGmGroup5Ptr::createMetGmGroup5PtrForReading(const std::shared_ptr<MetGmGroup3Ptr> gp3,
  282                                                                                const std::shared_ptr<MetGmHDTag> hdTag)
  283 {
  284     switch (hdTag->asShort()) {
  285     case MetGmHDTag::HD_2D:
  286     case MetGmHDTag::HD_2D_T:
  287     case MetGmHDTag::HD_3D_T: {
  288         shared_array<float> data(new float[hdTag->totalSize()]);
  289 
  290         std::shared_ptr<MetGmGroup5Ptr> gp5(new MetGmGroup5Ptr(gp3, hdTag, data));
  291 
  292         gp5->sOffset_ = ftell(gp3->mgmHandle()->fileHandle()->handle());
  293         MGM_THROW_ON_ERROR(mgm_read_group5(*gp3->mgmHandle()->fileHandle(), *gp3->mgmHandle(), data.get()))
  294         gp5->eOffset_ = ftell(gp3->mgmHandle()->fileHandle()->handle());
  295 
  296         MGM_THROW_ON_ERROR(mgm_param_is_convertible(gp3->p_id(), *gp3->mgmHandle()->version()))
  297 
  298         /*
  299          *           reindex data
  300          *     METGM               Fimex
  301          * (z, x, y, slice) -> (x, y, z, slice)
  302          */
  303         gp5->toFimexLayout();
  304         return gp5;
  305     } break;
  306     case MetGmHDTag::HD_0D:
  307     case MetGmHDTag::HD_0D_T:
  308     case MetGmHDTag::HD_1D:
  309     case MetGmHDTag::HD_1D_T:
  310     case MetGmHDTag::HD_3D:
  311     default: {
  312         std::ostringstream msg;
  313         msg << __FUNCTION__ << ": dimensionality not supported yet :" << hdTag->asString() << " for p_id =" << gp3->p_id();
  314         throw CDMException(msg.str());
  315     }
  316     }
  317 
  318     return std::shared_ptr<MetGmGroup5Ptr>(new MetGmGroup5Ptr(gp3, hdTag, shared_array<float>()));
  319     }
  320 
  321     std::shared_ptr<MetGmGroup5Ptr> MetGmGroup5Ptr::createMetGmGroup5PtrForSlicedReading(const std::shared_ptr<MetGmGroup3Ptr> gp3,
  322                                                                                          const std::shared_ptr<MetGmHDTag> hdTag)
  323     {
  324         switch(hdTag->asShort()) {
  325             case MetGmHDTag::HD_2D:
  326             case MetGmHDTag::HD_2D_T:
  327             case MetGmHDTag::HD_3D_T:
  328                 {
  329                 std::shared_ptr<MetGmGroup5Ptr> gp5(new MetGmGroup5Ptr(gp3, hdTag, shared_array<float>()));
  330 
  331                 /**
  332                  * must skip data to keep METGM C API lib satisifed
  333                  */
  334 
  335                 gp5->sOffset_ = ftell(gp3->mgmHandle()->fileHandle()->handle());
  336                 MGM_THROW_ON_ERROR(mgm_skip_group5(*gp3->mgmHandle()->fileHandle(), *gp3->mgmHandle()))
  337                 gp5->eOffset_ = ftell(gp3->mgmHandle()->fileHandle()->handle());
  338 
  339                 return gp5;
  340                 }
  341                 break;
  342             case MetGmHDTag::HD_0D:
  343             case MetGmHDTag::HD_0D_T:
  344             case MetGmHDTag::HD_1D:
  345             case MetGmHDTag::HD_1D_T:
  346             case MetGmHDTag::HD_3D:
  347             default: {
  348                 std::ostringstream msg;
  349                 msg << __FUNCTION__ << ": dimensionality not supported yet :" << hdTag->asString() << " for p_id =" << gp3->p_id();
  350                 throw CDMException(msg.str());
  351             }
  352         }
  353 
  354         return std::shared_ptr<MetGmGroup5Ptr>(new MetGmGroup5Ptr(gp3, hdTag, shared_array<float>()));
  355     }
  356 
  357     /*
  358      *  reindex data - several slices
  359      *     METGM               Fimex
  360      * (z, x, y, slice) -> (x, y, z, slice)
  361      */
  362     void MetGmGroup5Ptr::slicesToFimexLayout(shared_array<float>& slices, size_t numberOfSlices)
  363     {
  364         if(hdTag_->asShort() !=  MetGmHDTag::HD_3D_T)
  365             return;
  366 
  367         shared_array<float> dataT(new float[hdTag_->sliceSize() * numberOfSlices]);
  368 
  369         float* pos = slices.get();
  370         float* posT = dataT.get();
  371 
  372         size_t nz = hdTag_->zSize();
  373         size_t ny = hdTag_->ySize();
  374         size_t nx = hdTag_->xSize();
  375 
  376         size_t Nxy  = nx * ny;
  377         size_t Nzx  = nz * nx;
  378         size_t Nxyz = Nxy * nz;
  379 
  380         size_t xyzLeap = -Nxyz;
  381 
  382         for(size_t slice_index = 0; slice_index < numberOfSlices; ++slice_index) {
  383 
  384             xyzLeap += Nxyz;
  385 
  386             size_t xLeap  = -nx;
  387             size_t zxLeap = -Nzx;
  388 
  389             for(size_t y = 0; y < ny; ++y) {
  390 
  391                 xLeap += nx;
  392                 zxLeap += Nzx;
  393 
  394                 size_t xyLeap = -Nxy;
  395 
  396                 for(size_t z = 0; z < nz; ++z) {
  397 
  398                     xyLeap += Nxy;
  399 
  400                     size_t zLeapF = -nz;
  401                     size_t zLeapB = nx * nz;
  402 
  403                     size_t midx = nx / 2 + nx % 2;
  404                     for(size_t x = 0; x < midx; ++x) {
  405 
  406                         zLeapF += nz;
  407                         zLeapB -= nz;
  408 
  409                         float valueFx = *(pos + z + zLeapF + zxLeap + xyzLeap);
  410                         *(posT + x + xLeap + xyLeap + xyzLeap) = valueFx;
  411 
  412                         float valueBx = *(pos + z + zLeapB + zxLeap + xyzLeap);
  413                         *(posT + (nx - 1 - x) + xLeap + xyLeap + xyzLeap) = valueBx;
  414                     } // x
  415                 } // z
  416             } // y
  417         } // slice
  418 
  419         slices = dataT;
  420     }
  421 
  422     shared_array<float> MetGmGroup5Ptr::readDataSlices(size_t pos, size_t numberOfSlices)
  423     {
  424         assert( (pos + numberOfSlices - 1) >= 1 && (pos + numberOfSlices - 1) <= hdTag_->tSize());
  425 
  426         switch(hdTag_->asShort()) {
  427             case MetGmHDTag::HD_2D:
  428             case MetGmHDTag::HD_2D_T:
  429             case MetGmHDTag::HD_3D_T:
  430                 {
  431                     FILE* fh = fopen(pGp3_->mgmHandle()->fileHandle()->fileName().c_str(), "rb");;
  432                     if(!fh) {
  433                         return shared_array<float>();
  434                     }
  435 
  436                     mgm_handle* mh = mgm_new_handle();
  437                     if(!mh) {
  438                         fclose(fh);
  439                         return shared_array<float>();
  440                     }
  441 
  442                     mgm_group3* gp3 = mgm_new_group3();
  443                     if(!gp3) {
  444                         mgm_free_handle(mh);
  445                         fclose(fh);
  446                         return shared_array<float>();
  447                     }
  448 
  449                     int call_result = MGM_OK;
  450 
  451                     int n = 0;
  452                     int np = 0;
  453                     int ndp = 0;
  454 
  455                     call_result = mgm_read_header(fh, mh);
  456 
  457                     if(call_result != MGM_OK) {
  458                         mgm_free_group3(gp3);
  459                         mgm_free_handle(mh);
  460                         fclose(fh);
  461                         return shared_array<float>();
  462                     }
  463 
  464                     np = mgm_get_number_of_params(mh);
  465                     ndp = mgm_get_number_of_dist_params(mh);
  466 
  467                     for (n = 0; n < np; n++)
  468                     {
  469                         call_result = mgm_read_group3(fh, mh, gp3);
  470 
  471                         if(call_result != MGM_OK) {
  472                             mgm_free_group3(gp3);
  473                             mgm_free_handle(mh);
  474                             fclose(fh);
  475                             return shared_array<float>();
  476                         }
  477 
  478                         if (mgm_get_pz(gp3) > 0) {
  479                             call_result = mgm_skip_group4(fh, mh);
  480                             if(call_result != MGM_OK) {
  481                                 mgm_free_group3(gp3);
  482                                 mgm_free_handle(mh);
  483                                 fclose(fh);
  484                                 return shared_array<float>();
  485                             }
  486 
  487                         }
  488                         long cOffset = ftell(fh);
  489                         if(cOffset < sOffset_) {
  490                             // skip group5 data
  491                             int nt = mgm_get_nt(gp3);
  492                             for(int slice_index = 1; slice_index <= nt; ++slice_index)
  493                             {
  494                                 size_t cSlicePos = 0;
  495                                 call_result = mgm_skip_group5_slice(fh, mh, &cSlicePos);
  496                                 if(call_result != MGM_OK) {
  497                                     mgm_free_group3(gp3);
  498                                     mgm_free_handle(mh);
  499                                     fclose(fh);
  500                                     return shared_array<float>();
  501                                 }
  502                             }
  503                         } else if(cOffset > sOffset_) {
  504                             // something wrong
  505                             mgm_free_group3(gp3);
  506                             mgm_free_handle(mh);
  507                             fclose(fh);
  508                             return shared_array<float>();
  509                         } else {
  510                             shared_array<float> data(new float[hdTag_->sliceSize() * numberOfSlices]);
  511 
  512                             size_t slices_read = 0;
  513 
  514                             for(size_t slice_index = 1; slice_index <= hdTag_->tSize(); ++slice_index)
  515                             {
  516                                 size_t cSlicePos = -1;
  517                                 call_result = mgm_read_group5_slice(fh, mh, data.get() + (slices_read * hdTag_->sliceSize()), &cSlicePos);
  518                                 if(call_result != MGM_OK) {
  519                                     mgm_free_group3(gp3);
  520                                     mgm_free_handle(mh);
  521                                     fclose(fh);
  522                                     return shared_array<float>();
  523                                 }
  524                                 if(cSlicePos >= pos && cSlicePos < pos + numberOfSlices -1) {
  525                                     ++slices_read;
  526                                     continue;
  527                                 } else if(cSlicePos == pos + numberOfSlices -1) {
  528                                     mgm_free_group3(gp3);
  529                                     mgm_free_handle(mh);
  530                                     fclose(fh);
  531 
  532                                     MGM_THROW_ON_ERROR(mgm_param_is_convertible(pGp3_->p_id(), *pGp3_->mgmHandle()->version()))
  533 
  534                                     /*
  535                                      * reindex data - slice
  536                                      *   METGM      Fimex
  537                                      * (z, x, y -> (x, y, z)
  538                                      */
  539                                     slicesToFimexLayout(data, numberOfSlices);
  540 
  541                                     return data;
  542                                 }
  543                             }
  544                         }
  545                     }
  546 
  547                     mgm_free_group3(gp3);
  548                     mgm_free_handle(mh);
  549                     fclose(fh);
  550                     return shared_array<float>();
  551                 }
  552                 break;
  553             case MetGmHDTag::HD_0D:
  554             case MetGmHDTag::HD_0D_T:
  555             case MetGmHDTag::HD_1D:
  556             case MetGmHDTag::HD_1D_T:
  557             case MetGmHDTag::HD_3D:
  558             default:
  559                 throw CDMException(std::string(__FUNCTION__) + std::string(": dimensionality not supported yet :") + hdTag_->asString() +
  560                                    " for p_id =" + type2string<int>(pGp3_->p_id()));
  561         }
  562 
  563         return shared_array<float>();
  564     }
  565 
  566     std::shared_ptr<MetGmGroup5Ptr> MetGmGroup5Ptr::createMetGmGroup5PtrForSlicedWriting(const CDMReader_p pCdmReader, const CDMVariable* pVariable,
  567                                                                                          const std::shared_ptr<MetGmGroup3Ptr> pg3)
  568     {
  569         std::shared_ptr<MetGmHDTag> hdtag = MetGmHDTag::createMetGmDimensionsTagForWriting(pCdmReader, pVariable);
  570 
  571         std::string mgmUnits;
  572 
  573         switch(hdtag->asShort()) {
  574             case MetGmHDTag::HD_2D:
  575             case MetGmHDTag::HD_3D:
  576             case MetGmHDTag::HD_3D_T:
  577                 {
  578                     if(pg3->p_id() == 7) {
  579                         /**
  580                           * METGM is tricky here as it can accomodate both m and hPa units
  581                           * depending if it is geopotential_height or pressure but the
  582                           * mgm_get_param_unit is, for some reason, always returning 'm'
  583                           *
  584                           * it should honor the pr settings
  585                           */
  586                         if(hdtag->zTag()->pr() == 2) {
  587                             /**
  588                               * If the value specified for pr is 2,
  589                               * and the values for pid=7 (if reported)
  590                               * are heights given in meters above MSL
  591                               */
  592                             mgmUnits = "m";
  593                         } else {
  594                             mgmUnits = "hPa";
  595                         }
  596                     } else {
  597                         mgmUnits = std::string(mgm_get_param_unit(pg3->p_id(), *(pg3->mgmHandle())));
  598                     }
  599 
  600                     const CDM& cdmRef = pCdmReader->getCDM();
  601                     const std::string varName = pVariable->getName();
  602 
  603                     Units unitsConvertor;
  604                     if(! unitsConvertor.areConvertible(cdmRef.getUnits(varName), mgmUnits)) {
  605                         std::string msg(" can't convert from ");
  606                         msg.append(cdmRef.getUnits(varName)).append(" to ").append(mgmUnits).append(" for variable ").append(varName).append(" -- excluding");
  607                         MGM_MESSAGE_POINT(msg)
  608                         return std::shared_ptr<MetGmGroup5Ptr>();
  609                     }
  610 
  611                     std::shared_ptr<MetGmGroup5Ptr> gp5(new MetGmGroup5Ptr(pg3, hdtag, shared_array<float>()));
  612                     gp5->units_ = mgmUnits;
  613 
  614                     return gp5;
  615                 }
  616                 break;
  617             case MetGmHDTag::HD_0D:
  618             case MetGmHDTag::HD_0D_T:
  619             case MetGmHDTag::HD_1D:
  620             case MetGmHDTag::HD_1D_T:
  621             case MetGmHDTag::HD_2D_T:
  622             default:
  623             throw CDMException(std::string(__FUNCTION__) + std::string(": dimensionality not supported yet :") + hdtag->asString() + " for " + pVariable->getName());
  624         }
  625     }
  626 
  627     /*
  628      * reindex data - one slice
  629      *   Fimex      METGM
  630      * (x, y, z) -> (z, x, y)
  631      */
  632     void MetGmGroup5Ptr::sliceToMetGmLayout(shared_array<float>& slice)
  633     {
  634         if(hdTag_->asShort() !=  MetGmHDTag::HD_3D_T)
  635             return;
  636 
  637         shared_array<float> dataT(new float[hdTag_->sliceSize()]);
  638 
  639         float* pos = slice.get();
  640         float* posT = dataT.get();
  641 
  642         size_t nz = hdTag_->zSize();
  643         size_t ny = hdTag_->ySize();
  644         size_t nx = hdTag_->xSize();
  645 
  646         size_t Nxy  = nx * ny;
  647         size_t Nzx  = nz * nx;
  648 
  649         size_t xLeap  = -nx;
  650         size_t zxLeap = -Nzx;
  651 
  652         for(size_t y = 0; y < ny; ++y) {
  653 
  654             xLeap += nx;
  655             zxLeap += Nzx;
  656 
  657             size_t zLeapF = -nz;
  658 
  659             for(size_t x = 0; x < nx; ++x) {
  660 
  661                 zLeapF += nz;
  662 
  663                 size_t xyLeapF = -Nxy;
  664                 size_t xyLeapB = nz * Nxy;
  665 
  666                 size_t midz = nz / 2 + nz % 2;
  667                 for(size_t z = 0; z < midz; ++z) {
  668 
  669                     xyLeapF += Nxy;
  670                     xyLeapB -= Nxy;
  671 
  672                     float valueFz = *(pos + x + xLeap + xyLeapF);
  673                     *(posT + z + zLeapF + zxLeap) = valueFz;
  674                     float valueBz = *(pos + x + xLeap + xyLeapB);
  675                     *(posT + (nz - 1 - z) + zLeapF + zxLeap) = valueBz;
  676                 } // z
  677 
  678             } // x
  679 
  680         } // y
  681 
  682         slice = dataT;
  683     }
  684 
  685     void MetGmGroup5Ptr::dumpFimexLayout()
  686     {
  687         std::cerr << "dumping group5 in Fimex layout [START]" << std::endl;
  688         for(size_t sIndex = 0; sIndex < hdTag_->tSize(); ++sIndex) {
  689 
  690             float* slice = data_.get() + sIndex * hdTag_->sliceSize();
  691 
  692             for(size_t z_index = 0; z_index < hdTag_->zSize(); ++z_index) {
  693 
  694                 for(size_t y_index = 0; y_index < hdTag_->ySize(); ++y_index) {
  695 
  696                     for(size_t x_index = 0; x_index < hdTag_->xSize(); ++x_index) {
  697 
  698                         std::cerr << "[pid=" << pGp3_->p_id()  << "]"
  699                                   << "[T=" << sIndex  << "]"
  700                                   << "[X=" << x_index << "]"
  701                                   << "[Y=" << y_index << "]"
  702                                   << "[Z=" << z_index << "]"
  703                                   << "{" << slice[x_index + y_index * hdTag_->xSize() + z_index * hdTag_->xSize() * hdTag_->ySize()] << "}"
  704                                   << "    ";
  705                     } // x_index
  706 
  707                     std::cerr << std::endl << "end of Y" << std::endl;
  708 
  709                 } // y_index
  710 
  711                 std::cerr << std::endl << "end of Z" << std::endl;
  712 
  713             } // z_index
  714 
  715             std::cerr << std::endl << "end of T" << std::endl;
  716 
  717         } // sliceIndex
  718         std::cerr << "dumping group5 in Fimex layout"   << std::endl;
  719     }
  720 
  721 }