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

Last change on this file since 782 was 782, checked in by mhnguyen, 8 years ago

Reading attributes of curvilinear grid from file

+) Correct some minor bugs detecting type of grid
+) Use constant string for attributes conforming to CF convention
+) Add part of code to read attributes of curvilinear grid

Test
+) On Curie
+) test_remap passes

  • 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: 13.8 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*/)
12    : SuperClass()
13    , SuperClassWriter(filename, &comm_file, multifile)
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    }
34
35    return 1;
36  }
37
38  void CNc4DataInput::readFieldData_(CField* field)
39  {
40    CContext* context = CContext::getCurrent();
41    CContextServer* server = context->server;
42
43    CGrid* grid = field->grid;
44
45    if (!grid->doGridHaveDataToWrite())
46      if (SuperClass::type==MULTI_FILE || !isCollective) return;
47
48    StdString fieldId = field->getFieldOutputName();
49
50    CArray<double,1> fieldData(grid->getWrittenDataSize());
51    if (!field->default_value.isEmpty()) fieldData = field->default_value;
52
53    switch (SuperClass::type)
54    {
55      case MULTI_FILE:
56        SuperClassWriter::getData(fieldData, fieldId, isCollective, field->getNStep() - 1);
57        break;
58      case ONE_FILE:
59      {
60/*
61        std::vector<int> nZoomBeginGlobal = grid->getDistributionServer()->getZoomBeginGlobal();
62        std::vector<int> nZoomBeginServer = grid->getDistributionServer()->getZoomBeginServer();
63        std::vector<int> nZoomSizeServer  = grid->getDistributionServer()->getZoomSizeServer();
64
65        int ssize = nZoomBeginGlobal.size();
66
67        std::vector<StdSize> start(ssize);
68        std::vector<StdSize> count(ssize);
69
70        for (int i = 0; i < ssize; ++i)
71        {
72          start[i] = nZoomBeginServer[ssize - i - 1] - nZoomBeginGlobal[ssize - i - 1];
73          count[i] = nZoomSizeServer[ssize - i - 1];
74        }
75*/
76
77        std::vector<int> nZoomBeginGlobal = grid->getDistributionServer()->getZoomBeginGlobal();
78        std::vector<int> nZoomBeginServer = grid->getDistributionServer()->getZoomBeginServer();
79        std::vector<int> nZoomSizeServer  = grid->getDistributionServer()->getZoomSizeServer();
80
81        std::vector<StdSize> start, count;
82
83        CArray<bool,1> axisDomainOrder = grid->axis_domain_order;
84        std::vector<StdString> domainList = grid->getDomainList();
85        std::vector<StdString> axisList   = grid->getAxisList();
86        int numElement = axisDomainOrder.numElements();
87        int idxDomain = domainList.size() - 1, idxAxis = axisList.size() - 1;
88        int idx = nZoomBeginGlobal.size() - 1;
89
90        start.reserve(nZoomBeginGlobal.size());
91        count.reserve(nZoomBeginGlobal.size());
92
93        for (int i = numElement - 1; i >= 0; --i)
94        {
95          if (axisDomainOrder(i))
96          {
97            CDomain* domain = CDomain::get(domainList[idxDomain]);
98            if ((domain->type) != CDomain::type_attr::unstructured)
99            {
100              start.push_back(nZoomBeginServer[idx] - nZoomBeginGlobal[idx]);
101              count.push_back(nZoomSizeServer[idx]);
102            }
103            --idx ;
104            start.push_back(nZoomBeginServer[idx] - nZoomBeginGlobal[idx]);
105            count.push_back(nZoomSizeServer[idx]);
106            --idx ;
107            --idxDomain;
108          }
109          else
110          {
111            start.push_back(nZoomBeginServer[idx] - nZoomBeginGlobal[idx]);
112            count.push_back(nZoomSizeServer[idx]);
113            --idx;
114           }
115        }
116
117        SuperClassWriter::getData(fieldData, fieldId, isCollective, field->getNStep() - 1, &start, &count);
118        break;
119      }
120    }
121
122    field->inputField(fieldData);
123
124    if (!field->scale_factor.isEmpty() || !field->add_offset.isEmpty())
125    {
126      double scaleFactor = 1.0, addOffset = 0.0;
127      if (!field->scale_factor.isEmpty()) scaleFactor = field->scale_factor;
128      if (!field->add_offset.isEmpty()) addOffset = field->add_offset;
129      field->invertScaleFactorAddOffset(scaleFactor, addOffset);
130    }
131  }
132
133  void CNc4DataInput::readFieldAttributes_(CField* field, bool readAttributeValues)
134  {
135    StdString fieldId = field->getFieldOutputName();
136
137    CGrid* grid = field->grid;
138
139    std::vector<CDomain*> domainP = grid->getDomains();
140    std::vector<CAxis*> axisP = grid->getAxis();
141    int gridDim = domainP.size() * 2 + axisP.size();
142
143    // Verify the compatibility of dimension of declared grid and real grid in file
144    int realGridDim = 1;
145    std::map<StdString, StdSize> dimSizeMap = SuperClassWriter::getDimensions(&fieldId);
146    realGridDim = SuperClassWriter::isTemporal(fieldId) ? dimSizeMap.size() - 1 : dimSizeMap.size();
147
148    if (gridDim != realGridDim)
149       ERROR("CNc4DataInput::readFieldAttributes_(CField* field, bool readAttributeValues)",
150        << "Field '" << fieldId << "' has incorrect dimension " << std::endl
151        << "Verify dimension of grid defined by 'grid_ref' or 'domain_ref'/'axis_ref' and dimension of grid in read file.");
152
153    // Remove unlimited dimension from the map, we dont need it anymore
154    if (SuperClassWriter::isTemporal(fieldId)) dimSizeMap.erase(SuperClassWriter::getUnlimitedDimensionName());
155    int mapSize = dimSizeMap.size() - 1;
156
157    // Now process domain and axis
158    CArray<bool,1> axisDomainOrder = grid->axis_domain_order;
159    int numElement = domainP.size() + axisP.size();
160    int elementPosition = 0;
161    int idxDomain = 0, idxAxis = 0;
162
163    std::pair<std::set<StdString>::iterator,bool> it;
164    for (int i = 0; i < numElement; ++i)
165    {
166      if(axisDomainOrder(i))
167      {
168        if (readAttributeValues)
169        {
170           it = readValueDomains_.insert(domainP[idxDomain]->getId());
171           if (it.second) readDomainAttributeValueFromFile(domainP[idxDomain], dimSizeMap, mapSize - 1 - elementPosition, fieldId);
172        }
173        else
174        {
175          it = readMetaDataDomains_.insert(domainP[idxDomain]->getId());
176          if (it.second) readDomainAttributesFromFile(domainP[idxDomain], dimSizeMap, mapSize - 1 - elementPosition, fieldId);
177        }
178        ++idxDomain;
179        elementPosition += 2;
180      }
181      else
182      {
183        if (readAttributeValues)
184        {
185          it = readValueAxis_.insert(axisP[idxAxis]->getId());
186          if (it.second) readAxisAttributeValueFromFile(axisP[idxAxis], dimSizeMap, mapSize - elementPosition, fieldId);
187        }
188        else
189        {
190          it = readMetaDataAxis_.insert(axisP[idxAxis]->getId());
191          if (it.second) readAxisAttributesFromFile(axisP[idxAxis], dimSizeMap, mapSize - elementPosition, fieldId);
192        }
193        ++idxAxis;
194        ++elementPosition;
195      }
196    }
197  }
198
199  /*!
200    Read attributes of a domain from a file
201    \param [in] domain domain whose attributes are read from the file
202    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
203    \param [in] emelentPosition position of domain in grid
204    \param [in] fieldId id (or name) associated with the grid
205  */
206  void CNc4DataInput::readDomainAttributeValueFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
207                                                       int elementPosition, const StdString& fieldId)
208  {
209    // There are some optional attributes of a domain to retrieve from file    // + lon lat?
210    std::map<StdString, StdSize>::const_iterator itMapNj = dimSizeMap.begin(), itMapNi,
211                                                 iteMap  = dimSizeMap.end();
212
213    for (int i = 0; i < elementPosition; ++i, ++itMapNj) {}
214    itMapNi = itMapNj; ++itMapNi;
215
216    if (this->isRectilinear(fieldId))
217    {
218      // Ok, try to read some f.. attributes such as longitude and latitude
219      domain->latvalue_rectilinear_read_from_file.resize(itMapNj->second);
220      std::vector<StdSize> nBeginLat(1, 0), nSizeLat(1, itMapNj->second);
221      readFieldVariableValue(domain->latvalue_rectilinear_read_from_file, itMapNj->first, nBeginLat, nSizeLat, true);
222
223      domain->lonvalue_rectilinear_read_from_file.resize(itMapNi->second);
224      std::vector<StdSize> nBeginLon(1, 0), nSizeLon(1, itMapNi->second);
225      readFieldVariableValue(domain->lonvalue_rectilinear_read_from_file, itMapNi->first, nBeginLon, nSizeLon, true);
226      domain->fillInRectilinearLonLat();
227    }
228    else if (this->isCurvilinear(fieldId))
229    {
230      int ni = domain->ni;
231      int nj = domain->nj;
232      std::vector<StdSize> nBeginLatLon(2), nSizeLatLon(2);
233      nBeginLatLon[0] = domain->jbegin.getValue(); nBeginLatLon[1] = domain->ibegin.getValue();
234      nSizeLatLon[0]  = nj; nSizeLatLon[1] = ni;
235
236      StdString latName = this->getLatCoordName(fieldId);
237      domain->latvalue_2d.resize(ni,nj);
238      readFieldVariableValue(domain->latvalue_2d, latName, nBeginLatLon, nSizeLatLon);
239      StdString lonName = this->getLonCoordName(fieldId);
240      domain->lonvalue_2d.resize(ni,nj);
241      readFieldVariableValue(domain->lonvalue_2d, lonName, nBeginLatLon, nSizeLatLon);
242
243      StdString boundsLatName = this->getBoundsId(latName);
244      if (0 == boundsLatName.compare(""))
245         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
246              << "Field '" << fieldId << std::endl
247              << "Trying to read attributes from curvilinear grid."
248              << "Latitude variable " << latName << " does not have bounds.");
249      StdString boundsLonName = this->getBoundsId(lonName);
250      if (0 == boundsLonName.compare(""))
251         ERROR("CNc4DataInput::readDomainAttributeValueFromFile(...)",
252              << "Field '" << fieldId << std::endl
253              << "Trying to read attributes from curvilinear grid."
254              << "Longitude variable " << lonName << " does not have bounds.");
255
256      domain->nvertex.setValue(4);
257      std::vector<StdSize> nBeginBndsLatLon(3), nSizeBndsLatLon(3);
258      nBeginBndsLatLon[0] = domain->jbegin.getValue(); nSizeBndsLatLon[0] = nj;
259      nBeginBndsLatLon[1] = domain->ibegin.getValue(); nSizeBndsLatLon[1] = ni;
260      nBeginBndsLatLon[2] = 0; nSizeBndsLatLon[2] = 4;
261
262      domain->bounds_lat_2d.resize(4,ni,nj);
263      readFieldVariableValue(domain->bounds_lat_2d, boundsLatName, nBeginBndsLatLon, nSizeBndsLatLon);
264      domain->bounds_lon_2d.resize(4,ni,nj);
265      readFieldVariableValue(domain->bounds_lon_2d, boundsLonName, nBeginBndsLatLon, nSizeBndsLatLon);
266    }
267    else if (this->isUnstructured(fieldId))
268    {
269
270    }
271  }
272
273  /*!
274    Read attribute value of a domain from a file
275    \param [in] domain domain whose attributes are read from the file
276    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
277    \param [in] emelentPosition position of domain in grid
278    \param [in] fieldId id (or name) associated with the grid
279  */
280  void CNc4DataInput::readDomainAttributesFromFile(CDomain* domain, std::map<StdString, StdSize>& dimSizeMap,
281                                                   int elementPosition, const StdString& fieldId)
282  {
283    // There are some mandatory attributes of a domain to retrieve from file
284    // + ni_glo, nj_glo
285    std::map<StdString, StdSize>::const_iterator itMapNj = dimSizeMap.begin(), itMapNi,
286                                                 iteMap  = dimSizeMap.end();
287    for (int i = 0; i < elementPosition; ++i, ++itMapNj) {}
288    itMapNi = itMapNj; ++itMapNi;
289
290    if (this->isRectilinear(fieldId) || this->isCurvilinear(fieldId))
291    {
292      domain->nj_glo.setValue(itMapNj->second);
293      domain->ni_glo.setValue((itMapNi)->second);
294    }
295    else if (this->isUnstructured(fieldId))
296    {
297
298    }
299  }
300
301  /*!
302    Read attributes of an axis from a file
303    \param [in] axis axis whose attributes are read from the file
304    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
305    \param [in] emelentPosition position of axis in grid
306    \param [in] fieldId id (or name) associated with the grid
307  */
308  void CNc4DataInput::readAxisAttributesFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
309                                                 int elementPosition, const StdString& fieldId)
310  {
311    std::map<StdString, StdSize>::const_iterator itMapN = dimSizeMap.begin(),
312                                                 iteMap  = dimSizeMap.end();
313    for (int i = 0; i < elementPosition; ++i, ++itMapN) {}
314    axis->n_glo.setValue(itMapN->second);
315  }
316
317  /*!
318    Read attribute value of an axis from a file
319    \param [in] axis axis whose attributes are read from the file
320    \param [in] dimSizeMap Dimensions and and their corresponding names and size read from file
321    \param [in] emelentPosition position of axis in grid
322    \param [in] fieldId id (or name) associated with the grid
323  */
324  void CNc4DataInput::readAxisAttributeValueFromFile(CAxis* axis, std::map<StdString, StdSize>& dimSizeMap,
325                                                    int elementPosition, const StdString& fieldId)
326  {
327    std::map<StdString, StdSize>::const_iterator itMapN = dimSizeMap.begin(),
328                                                 iteMap  = dimSizeMap.end();
329    for (int i = 0; i < elementPosition; ++i, ++itMapN) {}
330
331    { // Read axis value
332      std::vector<StdSize> nBegin(1, 0), nSize(1, itMapN->second);
333      CArray<double,1> readAxisValue(itMapN->second);
334      readFieldVariableValue(readAxisValue, itMapN->first, nBegin, nSize, true);
335      int begin = 0, n = itMapN->second;
336      if (!axis->begin.isEmpty()) begin = axis->begin.getValue();
337      if (!axis->n.isEmpty()) n = axis->n.getValue();
338      axis->value.resize(n);
339      for (int i = 0; i < n; ++i) axis->value(i) = readAxisValue(begin + i);
340    }
341  }
342
343  void CNc4DataInput::closeFile_(void)
344  {
345    SuperClassWriter::close();
346  }
347} // namespace xios
Note: See TracBrowser for help on using the repository browser.