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