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

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

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