source: XMLIO_V2/dev/dev_rv/src/XMLIO/NetCDF4_data_output.hpp @ 137

Last change on this file since 137 was 137, checked in by hozdoba, 14 years ago

mise à jour
File size: 17.2 KB
Line 
1#ifndef __NETCDF4_DATA_OUTPUT__
2#define __NETCDF4_DATA_OUTPUT__
3
4// Entête NetCDF pour le C++.
5#include <netcdfcpp.h>
6
7namespace XMLIOSERVER
8{
9
10   static const char* TimeName = "time";
11
12   class NetCDF4DataOutput : public AbstractDataOutput
13   {
14      public :
15
16         NetCDF4DataOutput(CFile* const _file)
17            : AbstractDataOutput(_file), dataFile(NULL)
18         { /* Ne rien faire de plus */ }
19
20         const NcFile* getDataFile(void) const { return (dataFile); }
21
22         void writeVarData(const string& id, const Array<float, 1>& vdata)
23         {
24
25         }
26
27         virtual ~NetCDF4DataOutput( void)
28         { if (dataFile != NULL) delete dataFile; }
29
30      protected :
31
32         Poco::HashMap<string, int> getAllLimitedDim(void) const
33         {
34            Poco::HashMap<string, int> dims;
35            const std::set<const CDomain*> sdom = this->getRelFile()->getEnabledDomains();
36            const std::set<const CAxis*>  saxis = this->getRelFile()->getEnabledAxis();
37
38            // std::cout << "Nombre de domaines disponibles pour le fichier : " << sdom.size()  << std::endl;
39            // std::cout << "Nombre d'axes disponibles pour le fichier : "      << saxis.size() << std::endl;
40
41            std::set<const CDomain*>::const_iterator itt;
42            std::set<const CAxis*>::const_iterator it;
43
44            for ( itt = sdom.begin() ; itt != sdom.end(); itt++ )
45            {
46               string domid = ((*itt)->name.hasValue()) ? (string)(*itt)->name : (*itt)->getId();
47               string lonid = (sdom.size() == 1)? string("lon"): string("lon_").append(domid);
48               string latid = (sdom.size() == 1)? string("lat"): string("lat_").append(domid);
49               dims[lonid]  = (*itt)->ni;
50               dims[latid]  = (*itt)->nj;
51            }
52
53            for ( it = saxis.begin() ; it != saxis.end(); it++ )
54            {
55               string axisid = ((*it)->name.hasValue()) ?
56                        (string)(*it)->name : (*it)->getId();
57               dims[axisid]  =  (*it)->size;
58            }
59
60            return (dims);
61         }
62
63         std::string getFileName(void)
64         {
65            return ((getRelFile()->name.hasValue()) ?
66              string(getRelFile()->name): getRelFile()->name.getId());
67         }
68
69         std::string getTimeStamp(void)
70         {
71            const int buffer_size = 100;
72            time_t rawtime;
73            struct tm * timeinfo = NULL;
74            char buffer [buffer_size];
75
76            time ( &rawtime );
77            timeinfo = localtime ( &rawtime );
78            strftime (buffer, buffer_size, "%Y-%b-%d %H:%M:%S %Z", timeinfo);
79
80            return (string(buffer));
81         }
82
83         void writeCoords(const string& id, const Array<double, 1>& cdata)
84         {
85            NcVar *cVar = dataFile->get_var(id.c_str());
86
87            if (!cVar->put(cdata.dataFirst(), cdata.size()))
88               throw XMLIOUndefinedValueException
89                  ("Impossible d'écrire les valeurs de coordonnées "+ id +" !");
90         }
91
92         void setTimeCoords(void)
93         {
94            NcVar *oVar = NULL;
95            Poco::HashMap<string, string> hm;
96            AbstractCalendar * ccal = Context::GetCurrentContext()->getCalendar();
97            const std::set<const FieldOperation*> sfo = this->getRelFile()->getEnabledOperation();
98            std::set<const FieldOperation*>::const_iterator it;
99
100            //std::cout << "Nombre d'opération disponibles pour le fichier : " << sfo.size() << std::endl;
101
102            for ( it = sfo.begin() ; it != sfo.end(); it++ )
103            {
104               const FieldOperation * const   fope = (*it);
105               const Duration& ffope = (*it)->getFreqOp();
106               ostringstream oss; oss << "t_" << *fope << "_" << ffope;
107               const std::string cid(oss.str());
108
109               if (!(oVar = dataFile->add_var(cid.c_str(), ncDouble, dataFile->get_dim(TimeName))))
110                  throw XMLIOUndefinedValueException
111                     ("Impossible d'ajouter la coordonnée temporelle "+cid+" !");
112
113               hm ["standard_name"] = "time";
114               hm ["title"]         = "Time";
115               hm ["long_name"]     = "Time axis";
116               hm ["units"]         = string("seconds since ").append(ccal->getInitDate().toString());
117               hm ["calendar"]      = ccal->getType();
118               hm ["time_origin"]   = ccal->getInitDate().toString();
119               this->addStringAttributesToVar_str(oVar, hm);
120               hm.clear();
121            }
122         }
123
124      protected : /* virtual */
125
126         virtual void initFile(void)
127         {
128            string filename = this->getFileName();
129
130            // Création du fichier ou remplacement si celui-ci existe déjà.
131            dataFile = new NcFile(filename.c_str(), NcFile::Replace);
132            if(!dataFile->is_valid())
133               throw XMLIOUndefinedValueException
134                  ("Impossible d'ouvrir le fichier '"+ filename +"' pour l'écriture des données' !");
135
136            // Ajout de quelques attributs globaux.
137            dataFile->add_att("conventions", "CF-1.1");
138            dataFile->add_att("file_name"  , getFileName().c_str());
139            if (getRelFile()->description.hasValue())
140               dataFile->add_att("file_desc" , string(getRelFile()->description).c_str());
141            dataFile->add_att("production" , "An IPSL model");
142            dataFile->add_att("timeStamp"  , getTimeStamp().c_str()) ;
143         }
144
145         virtual void setDims(void)
146         {
147            bool withTime = true;
148            const Poco::HashMap<string, int>& allDim = this->getAllLimitedDim();
149            Poco::HashMap<string, int>::ConstIterator it;
150
151            // Ajout des dimensions limitées.
152            for (it = allDim.begin() ; it != allDim.end(); it++)
153               if (NULL == dataFile->add_dim((*it).first.c_str(), (*it).second))
154                  throw XMLIOUndefinedValueException("Impossible d'ajouter la dimension "+ (*it).first +" !");
155
156            // Ajout de la dimension temporelle non limitée.
157            if (withTime)
158               if (NULL == dataFile->add_dim(TimeName))
159                  throw XMLIOUndefinedValueException("Impossible d'ajouter la dimension temporelle !");
160         }
161
162         virtual void setCoords(void)
163         {
164            NcVar *latVar = NULL, // variable de latitude.
165                  *lonVar = NULL, // variable de longitude.
166                  *othvar = NULL; // toute autre variable.
167
168            Poco::HashMap<string, string> hm;
169            Poco::HashMap<string, float>  hmf;
170
171            const std::set<const CDomain*> sdom = this->getRelFile()->getEnabledDomains();
172            const std::set<const CAxis*>  saxis = this->getRelFile()->getEnabledAxis();
173
174            std::set<const CDomain*>::const_iterator itt;
175            for ( itt = sdom.begin() ; itt != sdom.end(); itt++ )
176            {
177               string domid = ((*itt)->name.hasValue()) ? (string)(*itt)->name : (*itt)->getId();
178               string lonid = (sdom.size() == 1)? string("lon"): string("lon_").append(domid);
179               string latid = (sdom.size() == 1)? string("lat"): string("lat_").append(domid);
180
181               if (!(latVar = dataFile->add_var(latid.c_str(), ncFloat, dataFile->get_dim(latid.c_str()))))
182                  throw XMLIOUndefinedValueException("Impossible d'ajouter la variable de latitude !");
183
184               if (!(lonVar = dataFile->add_var(lonid.c_str(), ncFloat, dataFile->get_dim(lonid.c_str()))))
185                  throw XMLIOUndefinedValueException("Impossible d'ajouter la variable de longitude !");
186
187               // Attribut de latitude.
188               hm["axis"]          = "Y" ;
189               hm["standard_name"] = "latitude" ;
190               hm["units"]         = "degrees_north";
191               hm["long_name"]     = "Latitude" ;
192               hm["nav_model"]     = domid ;
193               if ((*itt)->latvalue.hasValue())
194               {
195                  hmf["valid_min"] = ((Array<double, 1>)(*itt)->latvalue)(0);
196                  hmf["valid_max"] = ((Array<double, 1>)(*itt)->latvalue)(((Array<double, 1>)(*itt)->latvalue).size()-1);
197                  this->addStringAttributesToVar(latVar, hmf);
198                  hmf.clear();
199               }
200               this->addStringAttributesToVar_str(latVar, hm);
201               hm.clear();
202
203               // Attribut de longitude.
204               hm["axis"]          = "X" ;
205               hm["standard_name"] = "longitude" ;
206               hm["units"]         = "degrees_east";
207               hm["long_name"]     = "Longitude" ;
208               hm["nav_model"]     = domid ;
209               if ((*itt)->lonvalue.hasValue())
210               {
211                  hmf["valid_min"] = ((Array<double, 1>)(*itt)->lonvalue)(0);
212                  hmf["valid_max"] = ((Array<double, 1>)(*itt)->lonvalue)(((Array<double, 1>)(*itt)->lonvalue).size()-1);
213                  this->addStringAttributesToVar(lonVar, hmf);
214                  hmf.clear();
215               }
216               this->addStringAttributesToVar_str(lonVar, hm);
217               hm.clear();
218
219               if ((*itt)->lonvalue.hasValue())
220                  this->writeCoords(lonid, (*itt)->lonvalue);
221               else throw XMLIOUndefinedValueException
222                  ("Les coordonnées de longitude (xvalue) ne sont pas définies pour le domaine \""+domid+"\".");
223
224               if ((*itt)->latvalue.hasValue())
225                  this->writeCoords(latid, (*itt)->latvalue);
226               else throw XMLIOUndefinedValueException
227                  ("Les coordonnées de latitude (yvalue) ne sont pas définies pour le domaine \""+domid+"\".");
228            }
229
230            // Attribut des autres coordonnées.
231            std::set<const CAxis*>::const_iterator it;
232            for ( it = saxis.begin() ; it != saxis.end(); it++ )
233            {
234               string axisid = ((*it)->name.hasValue()) ? (string)(*it)->name : (*it)->getId();
235               if (!(othvar = dataFile->add_var(axisid.c_str(), ncFloat, dataFile->get_dim(axisid.c_str()))))
236                  throw XMLIOUndefinedValueException("Impossible d'ajouter la variable "+ (*it)->getId() +" !");
237
238               hm["axis"]     = "Z" ;
239               hm["title"]    = axisid ;
240               hm["positive"] = "unknown" ;
241
242               if ((*it)->standard_name.hasValue())hm["standard_name"] = (*it)->standard_name ;
243               if ((*it)->long_name.hasValue())    hm["long_name"]     = (*it)->long_name ;
244               if ((*it)->unit.hasValue())         hm["units"]         = (*it)->unit;
245               if ((*it)->zvalue.hasValue())
246               {
247                  hmf["valid_min"] = ((Array<double, 1>)(*it)->zvalue)(0);
248                  hmf["valid_max"] = ((Array<double, 1>)(*it)->zvalue)(((Array<double, 1>)(*it)->zvalue).size()-1);
249                  this->addStringAttributesToVar(othvar, hmf);
250                  hmf.clear();
251               }
252               this->addStringAttributesToVar_str(othvar, hm);
253               hm.clear();
254
255               if ((*it)->zvalue.hasValue())
256                  this->writeCoords(axisid, (*it)->zvalue);
257               else throw XMLIOUndefinedValueException
258                  ("Les coordonnées de l'axe \""+axisid+"\" (value) ne sont pas définies.");
259            }
260
261            this->setTimeCoords();
262         }
263
264         virtual void setVars(void)
265         {
266            NcType tvar = ncFloat;
267            NcVar *var  = NULL;
268
269            Poco::HashMap<string, string> hm;
270            Poco::HashMap<string, float>  hmf;
271
272            bool lonlat = false, wtime = true;
273
274            const std::vector<CField*>& enabledFields = getRelFile()->getEnabledFields();
275            const std::set<const CDomain*> sdom = this->getRelFile()->getEnabledDomains();
276
277            std::vector<CField*>::const_iterator it;
278
279            for ( it = enabledFields.begin() ; it != enabledFields.end(); it++ )
280            {
281                     CField  * const  field = (*it);
282               const CField  * const bfield = (*it)->getBaseObject();
283               const CDomain * const   rdom = field->getGrid()->getRelDomain();
284               const CAxis   * const  raxis = field->getGrid()->getRelAxis();
285               FieldOperation* const   fope = field->operation.getValue();
286               //Duration    * const  ffope = field->freq_op.getValue();
287
288               string fieldid = (field->name.hasValue()) ?
289                                 (string)field->name : bfield->getId();
290               string domid   = (rdom ->name.hasValue()) ?
291                                 (string)rdom ->name : rdom  ->getId();
292               string axisid  = (raxis->name.hasValue()) ?
293                                 (string)raxis->name : raxis ->getId();
294
295               string lonid = (sdom.size() == 1)? string("lon")
296                     : string("lon_").append(domid); // Nom de la coordonnée longitudinale associée.
297               string latid = (sdom.size() == 1)? string("lat")
298                     : string("lat_").append(domid); // Nom de la coordonnée latitudinale  associée.
299
300               lonlat = !field->getGrid()->_hasAxis();
301
302               if (fope != NULL)
303                  if (fope->getId().compare("once") == 0)
304                     wtime = false;
305
306               tvar = ncFloat;
307               if (field->prec.hasValue())
308                  if (field->prec == 8)
309                     tvar = ncDouble;
310
311               if (wtime)
312               {
313                  if (lonlat) // 2D spatio + temps
314                  {
315                     if (!(var = dataFile->add_var(fieldid.c_str(), tvar,
316                                 dataFile->get_dim(TimeName),
317                                 dataFile->get_dim(latid.c_str()),
318                                 dataFile->get_dim(lonid.c_str()))))
319                        throw XMLIOUndefinedValueException("Impossible d'ajouter le champ "+ field->getId() +" !");
320                  }
321                  else // 3D spatio + temps
322                  {
323                     if (!(var = dataFile->add_var(fieldid.c_str(), tvar,
324                                 dataFile->get_dim(TimeName),
325                                 dataFile->get_dim(axisid.c_str()),
326                                 dataFile->get_dim(latid.c_str()),
327                                 dataFile->get_dim(lonid.c_str()))))
328                        throw XMLIOUndefinedValueException("Impossible d'ajouter le champ "+ field->getId() +" !");
329                  }
330               }
331               else
332               {
333                  if (lonlat) // 2D spatio sans temps
334                  {
335                     if (!(var = dataFile->add_var(fieldid.c_str(), tvar,
336                                 dataFile->get_dim(latid.c_str()),
337                                 dataFile->get_dim(lonid.c_str()))))
338                        throw XMLIOUndefinedValueException("Impossible d'ajouter le champ "+ field->getId() +" !");
339                  }
340                  else // 3D spatio sans temps
341                  {
342                     if (!(var = dataFile->add_var(fieldid.c_str(), tvar,
343                                 dataFile->get_dim(axisid.c_str()),
344                                 dataFile->get_dim(latid.c_str()),
345                                 dataFile->get_dim(lonid.c_str()))))
346                        throw XMLIOUndefinedValueException("Impossible d'ajouter le champ "+ field->getId() +" !");
347                  }
348               }
349
350               if (field->standard_name.hasValue())hm["standard_name"]     = field->standard_name ;
351               if (field->long_name.hasValue())    hm["long_name"]         = field->long_name ;
352               if (field->unit.hasValue())         hm["units"]             = field->unit;
353
354               if (field->operation.hasValue())
355               {
356                  hm["online_operation"]   = field->operation.getValue()->getId();
357
358                  if (!field->operation.getValue()->getFreqOp().isNone())
359                     hm["interval_operation"] = field->operation.getValue()->getFreqOp().toString();
360                  else
361                     hm["interval_operation"] = ((Duration)this->getRelFile()->output_freq).toString();
362               }
363               hm["interval_write"] = ((Duration)this->getRelFile()->output_freq).toString();
364
365               this->addStringAttributesToVar_str(var, hm);
366               hm.clear();
367            }
368         }
369
370      private :
371
372         template <class U>
373            void addStringAttributesToVar(NcVar * const var, const Poco::HashMap<string, U>& attr)
374         {
375            typename Poco::HashMap<string, U>::ConstIterator it;
376
377            for ( it = attr.begin() ; it != attr.end(); it++ )
378               if (!var->add_att((*it).first.c_str(), (*it).second))
379                  throw XMLIOUndefinedValueException
380                     ("Impossible d'ajouter l'attribut' "+ (*it).first +" à la variable "+ var->name() +" !");
381         }
382
383         void addStringAttributesToVar_str(NcVar * const var, const Poco::HashMap<string, string>& attr)
384         {
385            Poco::HashMap<string, string>::ConstIterator it;
386
387            for ( it = attr.begin() ; it != attr.end(); it++ )
388               if (!var->add_att((*it).first.c_str(), (*it).second.c_str()))
389                  throw XMLIOUndefinedValueException
390                     ("Impossible d'ajouter l'attribut' "+ (*it).first +" à la variable "+ var->name() +" !");
391         }
392
393         NcFile* dataFile;
394
395   }; //class NetCDF4DataOutput
396
397}// namespace XMLIOSERVER
398
399#endif //__NETCDF4_DATA_OUTPUT__
Note: See TracBrowser for help on using the repository browser.