1 /*
2  MDAL - Mesh Data Abstraction Library (MIT License)
3  Copyright (C) 2018 Lutra Consulting Ltd.
4 */
5 
6 #include <stddef.h>
7 #include <iosfwd>
8 #include <iostream>
9 #include <fstream>
10 #include <sstream>
11 #include <string>
12 #include <vector>
13 #include <map>
14 #include <cassert>
15 #include <memory>
16 
17 #include "mdal_binary_dat.hpp"
18 #include "mdal.h"
19 #include "mdal_utils.hpp"
20 #include "mdal_logger.hpp"
21 
22 #include <math.h>
23 
24 static const int CT_VERSION   = 3000;
25 static const int CT_OBJTYPE   = 100;
26 static const int CT_SFLT      = 110;
27 static const int CT_SFLG      = 120;
28 static const int CT_BEGSCL    = 130;
29 static const int CT_BEGVEC    = 140;
30 static const int CT_VECTYPE   = 150;
31 static const int CT_OBJID     = 160;
32 static const int CT_NUMDATA   = 170;
33 static const int CT_NUMCELLS  = 180;
34 static const int CT_NAME      = 190;
35 static const int CT_TS        = 200;
36 static const int CT_ENDDS     = 210;
37 static const int CT_RT_JULIAN = 240;
38 static const int CT_TIMEUNITS = 250;
39 static const int CT_2D_MESHES = 3;
40 static const int CT_FLOAT_SIZE = 4;
41 static const int CF_FLAG_SIZE = 1;
42 static const int CF_FLAG_INT_SIZE = 4;
43 
44 
exit_with_error(MDAL_Status error,const std::string msg)45 static void exit_with_error( MDAL_Status error, const std::string msg )
46 {
47   MDAL::Log::error( error, "BINARY_DAT", msg );
48 }
49 
read(std::ifstream & in,char * s,int n)50 static bool read( std::ifstream &in, char *s, int n )
51 {
52   in.read( s, n );
53   if ( !in )
54     return true; //error
55   else
56     return false; //OK
57 }
58 
readIStat(std::ifstream & in,int sflg,char * flag)59 static bool readIStat( std::ifstream &in, int sflg, char *flag )
60 {
61   if ( sflg == CF_FLAG_SIZE )
62   {
63     in.read( flag, sflg );
64     if ( !in )
65       return true; // error
66   }
67   else
68   {
69     int istat;
70     in.read( reinterpret_cast< char * >( &istat ), sflg );
71     if ( !in )
72       return true; // error
73     else
74       *flag = ( istat == 1 );
75   }
76   return false;
77 }
78 
DriverBinaryDat()79 MDAL::DriverBinaryDat::DriverBinaryDat():
80   Driver( "BINARY_DAT",
81           "Binary DAT",
82           "*.dat",
83           Capability::ReadDatasets | Capability::WriteDatasetsOnVertices
84         )
85 {
86 }
87 
create()88 MDAL::DriverBinaryDat *MDAL::DriverBinaryDat::create()
89 {
90   return new DriverBinaryDat();
91 }
92 
93 MDAL::DriverBinaryDat::~DriverBinaryDat() = default;
94 
canReadDatasets(const std::string & uri)95 bool MDAL::DriverBinaryDat::canReadDatasets( const std::string &uri )
96 {
97   std::ifstream in = MDAL::openInputFile( uri, std::ifstream::in | std::ifstream::binary );
98   int version;
99 
100   if ( read( in, reinterpret_cast< char * >( &version ), 4 ) )
101     return false;
102 
103   if ( version != CT_VERSION ) // Version should be 3000
104     return false;
105 
106   return true;
107 }
108 
109 /**
110  * The DAT format contains "datasets" and each dataset has N-outputs. One output
111  * represents data for all vertices/faces for one timestep
112  *
113  * in TUFLOW results there could be also a special timestep (99999) with maximums
114  * we will put it into a separate dataset with name suffixed with "/Maximums"
115  *
116  * In MDAL we convert one output to one MDAL dataset;
117  *
118  */
load(const std::string & datFile,MDAL::Mesh * mesh)119 void MDAL::DriverBinaryDat::load( const std::string &datFile, MDAL::Mesh *mesh )
120 {
121   mDatFile = datFile;
122   MDAL::Log::resetLastStatus();
123 
124   if ( !MDAL::fileExists( mDatFile ) )
125   {
126     MDAL::Log::error( MDAL_Status::Err_FileNotFound, name(), "File could not be found " + mDatFile );
127     return;
128   }
129 
130   std::ifstream in = MDAL::openInputFile( mDatFile, std::ifstream::in | std::ifstream::binary );
131 
132   // implementation based on information from:
133   // http://www.xmswiki.com/wiki/SMS:Binary_Dataset_Files_*.dat
134   if ( !in ) return exit_with_error( MDAL_Status::Err_FileNotFound, "Couldn't open the file" );
135 
136   size_t vertexCount = mesh->verticesCount();
137   size_t elemCount = mesh->facesCount();
138 
139   int card = 0;
140   int version;
141   int objecttype;
142   int sflt;
143   int sflg = 0;
144   int vectype;
145   int objid;
146   int numdata;
147   int numcells;
148   char groupName[40];
149   double referenceTime;
150   int timeUnit = 0;
151   std::string timeUnitStr;
152   char istat;
153   float time;
154 
155   if ( read( in, reinterpret_cast< char * >( &version ), 4 ) ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "Unable to read version" );
156 
157   if ( version != CT_VERSION ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "Invalid version " );
158 
159   std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
160                                           name(),
161                                           mesh,
162                                           mDatFile
163                                         ); // DAT datasets
164   group->setDataLocation( MDAL_DataLocation::DataOnVertices );
165 
166   // in TUFLOW results there could be also a special timestep (99999) with maximums
167   // we will put it into a separate dataset
168   std::shared_ptr<DatasetGroup> groupMax = std::make_shared< DatasetGroup >(
169         name(),
170         mesh,
171         mDatFile
172       );
173   groupMax->setDataLocation( MDAL_DataLocation::DataOnVertices );
174 
175   while ( card != CT_ENDDS )
176   {
177     if ( read( in, reinterpret_cast< char * >( &card ), 4 ) )
178     {
179       // We've reached the end of the file and there was no ends card
180       break;
181     }
182 
183     switch ( card )
184     {
185       case CT_OBJTYPE:
186         // Object type
187         if ( read( in, reinterpret_cast< char * >( &objecttype ), 4 ) || objecttype != CT_2D_MESHES ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "Invalid object type" );
188         break;
189 
190       case CT_SFLT:
191         // Float size
192         if ( read( in, reinterpret_cast< char * >( &sflt ), 4 ) || sflt != CT_FLOAT_SIZE ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read float size" );
193         break;
194 
195       case CT_SFLG:
196         // Flag size
197         if ( read( in, reinterpret_cast< char * >( &sflg ), 4 ) )
198           if ( sflg != CF_FLAG_SIZE && sflg != CF_FLAG_INT_SIZE )
199             return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read flag size" );
200         break;
201 
202       case CT_BEGSCL:
203         group->setIsScalar( true );
204         groupMax->setIsScalar( true );
205         break;
206 
207       case CT_BEGVEC:
208         group->setIsScalar( false );
209         groupMax->setIsScalar( false );
210         break;
211 
212       case CT_VECTYPE:
213         // Vector type
214         if ( read( in, reinterpret_cast< char * >( &vectype ), 4 ) || vectype != 0 ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read vector type" );
215         break;
216 
217       case CT_OBJID:
218         // Object id
219         if ( read( in, reinterpret_cast< char * >( &objid ), 4 ) ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read object id" );
220         break;
221 
222       case CT_NUMDATA:
223         // Num data
224         if ( read( in, reinterpret_cast< char * >( &numdata ), 4 ) ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read num data" );
225         if ( numdata != static_cast< int >( vertexCount ) ) return exit_with_error( MDAL_Status::Err_IncompatibleMesh, "invalid num data" );
226         break;
227 
228       case CT_NUMCELLS:
229         // Num data
230         if ( read( in, reinterpret_cast< char * >( &numcells ), 4 ) ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read num cells" );
231         if ( numcells != static_cast< int >( elemCount ) ) return exit_with_error( MDAL_Status::Err_IncompatibleMesh, "invalid num cells" );
232         break;
233 
234       case CT_NAME:
235         // Name
236         if ( read( in, reinterpret_cast< char * >( &groupName ), 40 ) ) return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read name" );
237         if ( groupName[39] != 0 )
238           groupName[39] = 0;
239         group->setName( trim( std::string( groupName ) ) );
240         groupMax->setName( group->name() + "/Maximums" );
241         break;
242 
243       case CT_RT_JULIAN:
244         // Reference time
245         if ( readIStat( in, sflg, &istat ) )
246           return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read reference time" );
247 
248         if ( read( in, reinterpret_cast< char * >( &time ), 8 ) )
249           return exit_with_error( MDAL_Status::Err_UnknownFormat, "unable to read reference time" );
250 
251         referenceTime = static_cast<double>( time );
252         group->setReferenceTime( DateTime( referenceTime, DateTime::JulianDay ) );
253         break;
254 
255       case CT_TIMEUNITS:
256         // Time unit
257         if ( read( in, reinterpret_cast< char * >( &timeUnit ), 4 ) )
258           return exit_with_error( MDAL_Status::Err_UnknownFormat, "Unable to read time units" );
259 
260         switch ( timeUnit )
261         {
262           case 0:
263             timeUnitStr = "hours";
264             break;
265           case 1:
266             timeUnitStr = "minutes";
267             break;
268           case 2:
269             timeUnitStr = "seconds";
270             break;
271           case 4:
272             timeUnitStr = "days";
273             break;
274           default:
275             timeUnitStr = "unknown";
276             break;
277         }
278         group->setMetadata( "TIMEUNITS", timeUnitStr );
279         break;
280 
281       case CT_TS:
282         // Time step!
283         if ( readIStat( in, sflg, &istat ) )
284           return exit_with_error( MDAL_Status::Err_UnknownFormat, "Invalid time step" );
285 
286         if ( read( in, reinterpret_cast< char * >( &time ), 4 ) )
287           return exit_with_error( MDAL_Status::Err_UnknownFormat, "Invalid time step" );
288 
289         double rawTime = static_cast<double>( time );
290         MDAL::RelativeTimestamp t( rawTime, MDAL::parseDurationTimeUnit( timeUnitStr ) );
291 
292         if ( readVertexTimestep( mesh, group, groupMax, t, istat, sflg, in ) )
293           return exit_with_error( MDAL_Status::Err_UnknownFormat, "Unable to read vertex timestep" );
294 
295         break;
296     }
297   }
298 
299   if ( !group || group->datasets.size() == 0 )
300     return exit_with_error( MDAL_Status::Err_UnknownFormat, "No datasets" );
301 
302   group->setStatistics( MDAL::calculateStatistics( group ) );
303   mesh->datasetGroups.push_back( group );
304 
305   if ( groupMax && groupMax->datasets.size() > 0 )
306   {
307     groupMax->setStatistics( MDAL::calculateStatistics( groupMax ) );
308     mesh->datasetGroups.push_back( groupMax );
309   }
310 }
311 
readVertexTimestep(const MDAL::Mesh * mesh,std::shared_ptr<DatasetGroup> group,std::shared_ptr<DatasetGroup> groupMax,MDAL::RelativeTimestamp time,bool hasStatus,int sflg,std::ifstream & in)312 bool MDAL::DriverBinaryDat::readVertexTimestep(
313   const MDAL::Mesh *mesh,
314   std::shared_ptr<DatasetGroup> group,
315   std::shared_ptr<DatasetGroup> groupMax,
316   MDAL::RelativeTimestamp time,
317   bool hasStatus,
318   int sflg,
319   std::ifstream &in )
320 {
321   assert( group && groupMax && ( group->isScalar() == groupMax->isScalar() ) );
322   bool isScalar = group->isScalar();
323 
324   size_t vertexCount = mesh->verticesCount();
325   size_t faceCount = mesh->facesCount();
326 
327   std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MDAL::MemoryDataset2D >( group.get(), hasStatus );
328 
329   bool active = true;
330   for ( size_t i = 0; i < faceCount; ++i )
331   {
332     if ( hasStatus )
333     {
334       if ( readIStat( in, sflg, reinterpret_cast< char * >( &active ) ) )
335         return true; //error
336       dataset->setActive( i, active );
337     }
338   }
339 
340   for ( size_t i = 0; i < vertexCount; ++i )
341   {
342     if ( !isScalar )
343     {
344       float x, y;
345 
346       if ( read( in, reinterpret_cast< char * >( &x ), 4 ) )
347         return true; //error
348       if ( read( in, reinterpret_cast< char * >( &y ), 4 ) )
349         return true; //error
350 
351       dataset->setVectorValue( i, static_cast< double >( x ), static_cast< double >( y ) );
352     }
353     else
354     {
355       float scalar;
356 
357       if ( read( in, reinterpret_cast< char * >( &scalar ), 4 ) )
358         return true; //error
359 
360       dataset->setScalarValue( i, static_cast< double >( scalar ) );
361     }
362   }
363 
364   if ( MDAL::equals( time.value( MDAL::RelativeTimestamp::hours ), 99999.0 ) ) // Special TUFLOW dataset with maximus
365   {
366     dataset->setTime( time );
367     dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
368     groupMax->datasets.push_back( dataset );
369   }
370   else
371   {
372     dataset->setTime( time );
373     dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
374     group->datasets.push_back( dataset );
375   }
376   return false; //OK
377 }
378 
379 // ////////////////////////////////////////////
380 // WRITE
381 // ////////////////////////////////////////////
382 
writeRawData(std::ofstream & out,const char * s,int n)383 static bool writeRawData( std::ofstream &out, const char *s, int n )
384 {
385   out.write( s, n );
386   if ( !out )
387     return true; //error
388   else
389     return false; //OK
390 }
391 
persist(MDAL::DatasetGroup * group)392 bool MDAL::DriverBinaryDat::persist( MDAL::DatasetGroup *group )
393 {
394   assert( group->dataLocation() == MDAL_DataLocation::DataOnVertices );
395 
396   std::ofstream out = MDAL::openOutputFile( group->uri(), std::ofstream::out | std::ofstream::binary );
397 
398   // implementation based on information from:
399   // http://www.xmswiki.com/wiki/SMS:Binary_Dataset_Files_*.dat
400   if ( !out )
401     return true; // Couldn't open the file
402 
403   const Mesh *mesh = group->mesh();
404   size_t nodeCount = mesh->verticesCount();
405   size_t elemCount = mesh->facesCount();
406 
407   // version card
408   writeRawData( out, reinterpret_cast< const char * >( &CT_VERSION ), 4 );
409 
410   // objecttype
411   writeRawData( out, reinterpret_cast< const char * >( &CT_OBJTYPE ), 4 );
412   writeRawData( out, reinterpret_cast< const char * >( &CT_2D_MESHES ), 4 );
413 
414   // float size
415   writeRawData( out, reinterpret_cast< const char * >( &CT_SFLT ), 4 );
416   writeRawData( out, reinterpret_cast< const char * >( &CT_FLOAT_SIZE ), 4 );
417 
418   // Flag size
419   writeRawData( out, reinterpret_cast< const char * >( &CT_SFLG ), 4 );
420   writeRawData( out, reinterpret_cast< const char * >( &CF_FLAG_SIZE ), 4 );
421 
422   // Dataset Group Type
423   if ( group->isScalar() )
424   {
425     writeRawData( out, reinterpret_cast< const char * >( &CT_BEGSCL ), 4 );
426   }
427   else
428   {
429     writeRawData( out, reinterpret_cast< const char * >( &CT_BEGVEC ), 4 );
430   }
431 
432   // Object id (ignored)
433   int ignored_val = 1;
434   writeRawData( out, reinterpret_cast< const char * >( &CT_OBJID ), 4 );
435   writeRawData( out, reinterpret_cast< const char * >( &ignored_val ), 4 );
436 
437   // Num nodes
438   writeRawData( out, reinterpret_cast< const char * >( &CT_NUMDATA ), 4 );
439   writeRawData( out, reinterpret_cast< const char * >( &nodeCount ), 4 );
440 
441   // Num cells
442   writeRawData( out, reinterpret_cast< const char * >( &CT_NUMCELLS ), 4 );
443   writeRawData( out, reinterpret_cast< const char * >( &elemCount ), 4 );
444 
445   // Name
446   writeRawData( out, reinterpret_cast< const char * >( &CT_NAME ), 4 );
447   writeRawData( out, MDAL::leftJustified( group->name(), 39 ).c_str(), 40 );
448 
449   // Time steps
450   int istat = 1; // include if elements are active
451 
452   for ( size_t time_index = 0; time_index < group->datasets.size(); ++ time_index )
453   {
454     const std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::dynamic_pointer_cast<MDAL::MemoryDataset2D>( group->datasets[time_index] );
455 
456     writeRawData( out, reinterpret_cast< const char * >( &CT_TS ), 4 );
457     writeRawData( out, reinterpret_cast< const char * >( &istat ), 1 );
458     float ftime = static_cast<float>( dataset->time( RelativeTimestamp::hours ) );
459     writeRawData( out, reinterpret_cast< const char * >( &ftime ), 4 );
460 
461     if ( istat )
462     {
463       // Write status flags
464       for ( size_t i = 0; i < elemCount; i++ )
465       {
466         bool active = static_cast<bool>( dataset->active( i ) );
467         writeRawData( out, reinterpret_cast< const char * >( &active ), 1 );
468       }
469     }
470 
471     for ( size_t i = 0; i < nodeCount; i++ )
472     {
473       // Read values flags
474       if ( !group->isScalar() )
475       {
476         float x = static_cast<float>( dataset->valueX( i ) );
477         float y = static_cast<float>( dataset->valueY( i ) );
478         writeRawData( out, reinterpret_cast< const char * >( &x ), 4 );
479         writeRawData( out, reinterpret_cast< const char * >( &y ), 4 );
480       }
481       else
482       {
483         float val = static_cast<float>( dataset->scalarValue( i ) );
484         writeRawData( out, reinterpret_cast< const char * >( &val ), 4 );
485       }
486     }
487   }
488 
489   if ( writeRawData( out, reinterpret_cast< const char * >( &CT_ENDDS ), 4 ) ) return true;
490 
491   return false;
492 }
493 
writeDatasetOnFileSuffix() const494 std::string MDAL::DriverBinaryDat::writeDatasetOnFileSuffix() const
495 {
496   return "dat";
497 }
498