source: XIOS/trunk/src/io/nc4_data_input.cpp @ 811

Last change on this file since 811 was 811, checked in by ymipsl, 7 years ago

Bug fix when time dimension is not the unlimited dimension

YM

  • Property copyright set to
    Software name : XIOS (Xml I/O Server)
    http://forge.ipsl.jussieu.fr/ioserver
    Creation date : January 2009
    Licence : CeCCIL version2
    see license file in root directory : Licence_CeCILL_V2-en.txt
    or http://www.cecill.info/licences/Licence_CeCILL_V2-en.html
    Holder : CEA/LSCE (Laboratoire des Sciences du CLimat et de l'Environnement)
    CNRS/IPSL (Institut Pierre Simon Laplace)
    Project Manager : Yann Meurdesoif
    yann.meurdesoif@cea.fr
File size: 18.0 KB
Line 
1#include "nc4_data_input.hpp"
2
3#include "context.hpp"
4#include "context_server.hpp"
5#include "context_client.hpp"
6#include "domain.hpp"
7#include "axis.hpp"
8
9namespace xios
10{
11  CNc4DataInput::CNc4DataInput(const StdString& filename, MPI_Comm comm_file, bool multifile, bool isCollective /*= true*/, const StdString& timeCounterName /*= "time_counter"*/)
12    : SuperClass()
13    , SuperClassWriter(filename, &comm_file, multifile, timeCounterName)
14    , comm_file(comm_file)
15    , filename(filename)
16    , isCollective(isCollective)
17    , readMetaDataDomains_(), readValueDomains_()
18    , readMetaDataAxis_(), readValueAxis_()
19  {
20    SuperClass::type = multifile ? MULTI_FILE : ONE_FILE;
21  }
22
23  CNc4DataInput::~CNc4DataInput(void)
24  { /* Nothing more to do */ }
25
26  StdSize CNc4DataInput::getFieldNbRecords_(CField* field)
27  {
28    StdString fieldId = field->getFieldOutputName();
29
30    if (SuperClassWriter::isTemporal(fieldId))
31    {
32//      return SuperClassWriter::getDimensions(&fieldId)[SuperClassWriter::getUnlimitedDimensionName()];
33      return SuperClassWriter::getDimensions(&fieldId)[SuperClassWriter::getTimeCounterName()];
34    }
35
36    return 1;
37  }
38
39  void CNc4DataInput::readFieldData_(CField* field)
40  {
41    CContext* context = CContext::getCurrent();
42    CContextServer* server = context->server;
43
44    CGrid* grid = field->grid;
45
46    if (!grid->doGridHaveDataToWrite())
47      if (SuperClass::type==MULTI_FILE || !isCollective) return;
48
49    StdString fieldId = field->getFieldOutputName();
50
51    CArray<double,1> fieldData(grid->getWrittenDataSize());
52    if (!field->default_value.isEmpty()) fieldData = field->default_value;
53
54    switch (SuperClass::type)
55    {
56      case MULTI_FILE:
57        SuperClassWriter::getData(fieldData, fieldId, isCollective, field->getNStep() - 1);
58        break;
59      case ONE_FILE:
60      {
61/*
62        std::vector<int> nZoomBeginGlobal = grid->getDistributionServer()->getZoomBeginGlobal();
63        std::vector<int> nZoomBeginServer = grid->getDistributionServer()->getZoomBeginServer();
64        std::vector<int> nZoomSizeServer  = grid->getDistributionServer()->getZoomSizeServer();
65
66        int ssize = nZoomBeginGlobal.size();
67
68        std::vector<StdSize> start(ssize);
69        std::vector<StdSize> count(ssize);
70
71        for (int i = 0; i < ssize; ++i)
72        {
73          start[i] = nZoomBeginServer[ssize - i - 1] - nZoomBeginGlobal[ssize - i - 1];
74          count[i] = nZoomSizeServer[ssize - i - 1];
75        }
76*/
77
78        std::vector<int> nZoomBeginGlobal = grid->getDistributionServer()->getZoomBeginGlobal();
79        std::vector<int> nZoomBeginServer = grid->getDistributionServer()->getZoomBeginServer();
80        std::vector<int> nZoomSizeServer  = grid->getDistributionServer()->getZoomSizeServer();
81
82        std::vector<StdSize> start, count;
83
84        CArray<bool,1> axisDomainOrder = grid->axis_domain_order;
85        std::vector<StdString> domainList = grid->getDomainList();
86        std::vector<StdString> axisList   = grid->getAxisList();
87        int numElement = axisDomainOrder.numElements();
88        int idxDomain = domainList.size() - 1, idxAxis = axisList.size() - 1;
89        int idx = nZoomBeginGlobal.size() - 1;
90
91        start.reserve(nZoomBeginGlobal.size());
92        count.reserve(nZoomBeginGlobal.size());
93
94        for (int i = numElement - 1; i >= 0; --i)
95        {
96          if (axisDomainOrder(i))
97          {
98            CDomain* domain = CDomain::get(domainList[idxDomain]);
99            if ((domain->type) != CDomain::type_attr::unstructured)
100            {
101              start.push_back(nZoomBeginServer[idx] - nZoomBeginGlobal[idx]);
102              count.push_back(nZoomSizeServer[idx]);
103            }
104            --idx ;
105            start.push_back(nZoomBeginServer[idx] - nZoomBeginGlobal[idx]);
106            count.push_back(nZoomSizeServer[idx]);
107            --idx ;
108            --idxDomain;
109          }
110          else
111          {
112            start.push_back(nZoomBeginServer[idx] - nZoomBeginGlobal[idx]);
113            count.push_back(nZoomSizeServer[idx]);
114            --idx;
115           }
116        }
117
118        SuperClassWriter::getData(fieldData, fieldId, isCollective, field->getNStep() - 1, &start, &count);
119        break;
120      }
121    }
122
123    field->inputField(fieldData);
124
125    if (!field->scale_factor.isEmpty() || !field->add_offset.isEmpty())
126    {
127      double scaleFactor = 1.0, addOffset = 0.0;
128      if (!field->scale_factor.isEmpty()) scaleFactor = field->scale_factor;
129      if (!field->add_offset.isEmpty()) addOffset = field->add_offset;
130      field->invertScaleFactorAddOffset(scaleFactor, addOffset);
131    }
132  }
133
134  void CNc4DataInput::readFieldAttributes_(CField* field, bool readAttributeValues)
135  {
136    StdString fieldId = field->getFieldOutputName();
137
138    CGrid* grid = field->grid;
139
140    std::vector<CDomain*> domainP = grid->getDomains();
141    std::vector<CAxis*> axisP = grid->getAxis();
142    int gridDim = domainP.size() * 2 + axisP.size();
143
144    // Verify the compatibility of dimension of declared grid and real grid in file
145    int realGridDim = 1;
146    bool isUnstructuredGrid = SuperClassWriter::isUnstructured(fieldId);
147    std::map<StdString, StdSize> dimSizeMap = SuperClassWriter::getDimensions(&fieldId);
148    std::list<StdString> dimList = SuperClassWriter::getDimensionsList(&fieldId);
149   
150    realGridDim = SuperClassWriter::isTemporal(fieldId) ? dimSizeMap.size() - 1 : dimSizeMap.size();
151    if (isUnstructuredGrid) ++realGridDim;
152
153    if (gridDim != realGridDim)
154       ERROR("CNc4DataInput::readFieldAttributes_(CField* field, bool readAttributeValues)",
155        << "Field '" << fieldId << "' has incorrect dimension " << std::endl
156        << "Verify dimension of grid defined by 'grid_ref' or 'domain_ref'/'axis_ref' and dimension of grid in read file.");
157
158    // Remove unlimited dimension from the map, we dont need it anymore
159    if (SuperClassWriter::isTemporal(fieldId)) 
160    {
161      dimSizeMap.erase(SuperClassWriter::getUnlimitedDimensionName());
162      dimList.pop_front() ;  // assume time dimension is first
163    }
164   
165    int mapSize = dimSizeMap.size() - 1;
166    std::list<std::pair<StdString, StdSize> > listDimSize;
167/*
168    for (std::map<StdString, StdSize>::const_iterator itMap = dimSizeMap.begin(); itMap != dimSizeMap.end(); ++itMap)
169      listDimSize.push_front(*itMap);
170*/
171    for (std::list<StdString>::const_iterator it = dimList.begin(); it != dimList.end(); ++it)
172      listDimSize.push_front(*dimSizeMap.find(*it));
173
174    // Now process domain and axis
175    CArray<bool,1> axisDomainOrder = grid->axis_domain_order;
176    int numElement = domainP.size() + axisP.size();
177    int elementPosition = 0;
178    int idxDomain = 0, idxAxis = 0;
179
180    std::pair<std::set<StdString>::iterator,bool> it;
181    for (int i = 0; i < numElement; ++i)
182    {
183      if(axisDomainOrder(i))
184      {
185        if (readAttributeValues)
186        {
187           it = readValueDomains_.insert(domainP[idxDomain]->getId());
188           if (it.second) readDomainAttributeValueFromFile(domainP[idxDomain], listDimSize, elementPosition, fieldId);
189        }
190        else
191        {
192          it = readMetaDataDomains_.insert(domainP[idxDomain]->getId());
193          if (it.second) readDomainAttributesFromFile(domainP[idxDomain], listDimSize, elementPosition, fieldId);
194        }
195        ++idxDomain;
196        if (isUnstructuredGrid) ++elementPosition;
197        else elementPosition += 2;
198      }
199      else
200      {
201        if (readAttributeValues)
202        {
203          it = readValueAxis_.insert(axisP[idxAxis]->getId());
204          if (it.second) readAxisAttributeValueFromFile(axisP[idxAxis], listDimSize, elementPosition, fieldId);
205        }
206        else
207        {
208          it = readMetaDataAxis_.insert(axisP[idxAxis]->getId());
209          if (it.second) readAxisAttributesFromFile(axisP[idxAxis], listDimSize, elementPosition, fieldId);
210        }
211        ++idxAxis;
212        ++elementPosition;
213      }
214    }
215  }
216
217  /*!
218    Read attributes of a domain from a file
219    \param [in] domain domain whose attributes are read from the file
220    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
221    \param [in] emelentPosition position of domain in grid
222    \param [in] fieldId id (or name) associated with the grid
223  */
224  void CNc4DataInput::readDomainAttributeValueFromFile(CDomain* domain, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
225                                                       int elementPosition, const StdString& fieldId)
226  {
227    // There are some optional attributes of a domain to retrieve from file    // + lon lat?
228    std::list<std::pair<StdString, StdSize> >::const_iterator itMapNi = dimSizeMap.begin(), itMapNj,
229                                                              iteMap  = dimSizeMap.end();
230
231    for (int i = 0; i < elementPosition; ++i, ++itMapNi) {}
232    itMapNj = itMapNi; ++itMapNj;
233
234    if (this->isRectilinear(fieldId))
235    {
236      // Ok, try to read some f.. attributes such as longitude and latitude
237      domain->latvalue_rectilinear_read_from_file.resize(itMapNj->second);
238      std::vector<StdSize> nBeginLat(1, 0), nSizeLat(1, itMapNj->second);
239      readFieldVariableValue(domain->latvalue_rectilinear_read_from_file, itMapNj->first, nBeginLat, nSizeLat, true);
240
241      domain->lonvalue_rectilinear_read_from_file.resize(itMapNi->second);
242      std::vector<StdSize> nBeginLon(1, 0), nSizeLon(1, itMapNi->second);
243      readFieldVariableValue(domain->lonvalue_rectilinear_read_from_file, itMapNi->first, nBeginLon, nSizeLon, true);
244      domain->fillInRectilinearLonLat();
245    }
246    else if (this->isCurvilinear(fieldId))
247    {
248      int ni = domain->ni;
249      int nj = domain->nj;
250      std::vector<StdSize> nBeginLatLon(2), nSizeLatLon(2);
251      nBeginLatLon[0] = domain->jbegin.getValue(); nBeginLatLon[1] = domain->ibegin.getValue();
252      nSizeLatLon[0]  = nj; nSizeLatLon[1] = ni;
253
254      StdString latName = this->getLatCoordName(fieldId);
255      domain->latvalue_2d.resize(ni,nj);
256      readFieldVariableValue(domain->latvalue_2d, latName, nBeginLatLon, nSizeLatLon);
257      StdString lonName = this->getLonCoordName(fieldId);
258      domain->lonvalue_2d.resize(ni,nj);
259      readFieldVariableValue(domain->lonvalue_2d, lonName, nBeginLatLon, nSizeLatLon);
260
261      StdString boundsLatName = this->getBoundsId(latName);
262      if (0 == boundsLatName.compare(""))
263         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
264              << "Field '" << fieldId << std::endl
265              << "Trying to read attributes from curvilinear grid."
266              << "Latitude variable " << latName << " does not have bounds.");
267      StdString boundsLonName = this->getBoundsId(lonName);
268      if (0 == boundsLonName.compare(""))
269         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
270              << "Field '" << fieldId << std::endl
271              << "Trying to read attributes from curvilinear grid."
272              << "Longitude variable " << lonName << " does not have bounds.");
273
274      int nbVertex = this->getNbVertex(fieldId);
275      domain->nvertex.setValue(nbVertex);
276      std::vector<StdSize> nBeginBndsLatLon(3), nSizeBndsLatLon(3);
277      nBeginBndsLatLon[0] = domain->jbegin.getValue(); nSizeBndsLatLon[0] = nj;
278      nBeginBndsLatLon[1] = domain->ibegin.getValue(); nSizeBndsLatLon[1] = ni;
279      nBeginBndsLatLon[2] = 0; nSizeBndsLatLon[2] = nbVertex;
280
281      domain->bounds_lat_2d.resize(nbVertex,ni,nj);
282      readFieldVariableValue(domain->bounds_lat_2d, boundsLatName, nBeginBndsLatLon, nSizeBndsLatLon);
283      domain->bounds_lon_2d.resize(nbVertex,ni,nj);
284      readFieldVariableValue(domain->bounds_lon_2d, boundsLonName, nBeginBndsLatLon, nSizeBndsLatLon);
285    }
286    else if (this->isUnstructured(fieldId))
287    {
288      /*
289      if (domain->i_index.isEmpty())
290         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
291              << "Field '" << fieldId << std::endl
292              << "Trying to read attributes from unstructured grid."
293              << "i_index of domain" << domain->getId() << " is mandatory");
294     
295      int ni = domain->i_index.numElements();
296*/
297      int ni     = domain->ni;
298      int ibegin = domain->ibegin;
299      if (domain->i_index.isEmpty())
300      {
301        domain->i_index.resize(ni) ;
302        for(int idx = 0; idx < ni; ++idx) domain->i_index(idx)=ibegin+idx ;
303      }
304     
305      std::vector<StdSize> nBeginLatLon(1,0), nSizeLatLon(1,0);
306      nSizeLatLon[0]  = domain->ni_glo.getValue();
307      CArray<double,1> globalLonLat(domain->ni_glo.getValue());
308
309      StdString latName = this->getLatCoordName(fieldId);
310      readFieldVariableValue(globalLonLat, latName, nBeginLatLon, nSizeLatLon);
311      domain->latvalue_1d.resize(ni);
312      for (int idx = 0; idx < ni; ++idx)
313        domain->latvalue_1d(idx) =  globalLonLat(domain->i_index(idx));
314
315      StdString lonName = this->getLonCoordName(fieldId);
316      readFieldVariableValue(globalLonLat, lonName, nBeginLatLon, nSizeLatLon);
317      domain->lonvalue_1d.resize(ni);
318      for (int idx = 0; idx < ni; ++idx)
319        domain->lonvalue_1d(idx) = globalLonLat(domain->i_index(idx));
320
321      StdString boundsLatName = this->getBoundsId(latName);
322      if (0 == boundsLatName.compare(""))
323         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
324              << "Field '" << fieldId << std::endl
325              << "Trying to read attributes from unstructured grid."
326              << "Latitude variable " << latName << " does not have bounds.");
327      StdString boundsLonName = this->getBoundsId(lonName);
328      if (0 == boundsLonName.compare(""))
329         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
330              << "Field '" << fieldId << std::endl
331              << "Trying to read attributes from unstructured grid."
332              << "Longitude variable " << lonName << " does not have bounds.");
333
334      int nbVertex = this->getNbVertex(fieldId);
335      domain->nvertex.setValue(nbVertex);
336      std::vector<StdSize> nBeginBndsLatLon(2), nSizeBndsLatLon(2);
337      nBeginBndsLatLon[0] = 0; nSizeBndsLatLon[0] = domain->ni_glo.getValue();
338      nBeginBndsLatLon[1] = 0; nSizeBndsLatLon[1] = nbVertex;
339
340      CArray<double,2> globalBndsLonLat(nSizeBndsLatLon[1], nSizeBndsLatLon[0]);
341      readFieldVariableValue(globalBndsLonLat, boundsLatName, nBeginBndsLatLon, nSizeBndsLatLon);
342      domain->bounds_lat_1d.resize(nbVertex,ni);
343      for (int idx = 0; idx < ni; ++idx)
344        for (int jdx = 0; jdx < nbVertex; ++jdx)
345          domain->bounds_lat_1d(jdx,idx) = globalBndsLonLat(jdx, domain->i_index(idx));
346
347      readFieldVariableValue(globalBndsLonLat, boundsLonName, nBeginBndsLatLon, nSizeBndsLatLon);
348      domain->bounds_lon_1d.resize(nbVertex,ni);
349      for (int idx = 0; idx < ni; ++idx)
350        for (int jdx = 0; jdx < nbVertex; ++jdx)
351          domain->bounds_lon_1d(jdx,idx) = globalBndsLonLat(jdx, domain->i_index(idx));
352    }
353  }
354
355  /*!
356    Read attribute value of a domain from a file
357    \param [in] domain domain whose attributes are read from the file
358    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
359    \param [in] emelentPosition position of domain in grid
360    \param [in] fieldId id (or name) associated with the grid
361  */
362  void CNc4DataInput::readDomainAttributesFromFile(CDomain* domain, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
363                                                   int elementPosition, const StdString& fieldId)
364  {
365    // There are some mandatory attributes of a domain to retrieve from file
366    // + ni_glo, nj_glo
367    std::list<std::pair<StdString, StdSize> >::const_iterator itMapNi = dimSizeMap.begin(), itMapNj,
368                                                              iteMap  = dimSizeMap.end();
369    for (int i = 0; i < elementPosition; ++i, ++itMapNi) {}
370    itMapNj = itMapNi; ++itMapNj;
371
372    if (this->isRectilinear(fieldId) || this->isCurvilinear(fieldId))
373    {
374      domain->nj_glo.setValue(itMapNj->second);
375      domain->ni_glo.setValue(itMapNi->second);
376    }
377    else if (this->isUnstructured(fieldId))
378    {
379      domain->nj_glo.setValue(1);
380      domain->ni_glo.setValue(itMapNi->second);
381    }
382  }
383
384  /*!
385    Read attributes of an axis from a file
386    \param [in] axis axis whose attributes are read from the file
387    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
388    \param [in] emelentPosition position of axis in grid
389    \param [in] fieldId id (or name) associated with the grid
390  */
391  void CNc4DataInput::readAxisAttributesFromFile(CAxis* axis, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
392                                                 int elementPosition, const StdString& fieldId)
393  {
394    std::list<std::pair<StdString, StdSize> >::const_iterator itMapN = dimSizeMap.begin(),
395                                                              iteMap = dimSizeMap.end();
396    for (int i = 0; i < elementPosition; ++i, ++itMapN) {}
397    axis->n_glo.setValue(itMapN->second);
398  }
399
400  /*!
401    Read attribute value of an axis from a file
402    \param [in] axis axis whose attributes are read from the file
403    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
404    \param [in] emelentPosition position of axis in grid
405    \param [in] fieldId id (or name) associated with the grid
406  */
407  void CNc4DataInput::readAxisAttributeValueFromFile(CAxis* axis, std::list<std::pair<StdString, StdSize> >& dimSizeMap,
408                                                    int elementPosition, const StdString& fieldId)
409  {
410    std::list<std::pair<StdString, StdSize> >::const_iterator itMapN = dimSizeMap.begin(),
411                                                              iteMap = dimSizeMap.end();
412    for (int i = 0; i < elementPosition; ++i, ++itMapN) {}
413
414    { // Read axis value
415      std::vector<StdSize> nBegin(1, 0), nSize(1, itMapN->second);
416      CArray<double,1> readAxisValue(itMapN->second);
417      readFieldVariableValue(readAxisValue, itMapN->first, nBegin, nSize, true);
418      int begin = 0, n = itMapN->second;
419      if (!axis->begin.isEmpty()) begin = axis->begin.getValue();
420      if (!axis->n.isEmpty()) n = axis->n.getValue();
421      axis->value.resize(n);
422      for (int i = 0; i < n; ++i) axis->value(i) = readAxisValue(begin + i);
423    }
424  }
425
426  void CNc4DataInput::closeFile_(void)
427  {
428    SuperClassWriter::close();
429  }
430} // namespace xios
Note: See TracBrowser for help on using the repository browser.