1 /*
2 MDAL - Mesh Data Abstraction Library (MIT License)
3 Copyright (C) 2018 Peter Petrik (zilolv at gmail dot com)
4 */
5
6 #include "mdal_3di.hpp"
7 #include "mdal_logger.hpp"
8 #include <cmath>
9 #include <netcdf.h>
10 #include <assert.h>
11 #include "mdal_sqlite3.hpp"
12
13
Driver3Di()14 MDAL::Driver3Di::Driver3Di()
15 : DriverCF(
16 "3Di",
17 "3Di Results",
18 "results_3di.nc",
19 Capability::ReadMesh )
20 {
21 }
22
create()23 MDAL::Driver3Di *MDAL::Driver3Di::create()
24 {
25 return new Driver3Di();
26 }
27
buildUri(const std::string & meshFile)28 std::string MDAL::Driver3Di::buildUri( const std::string &meshFile )
29 {
30 mNcFile.reset( new NetCDFFile );
31
32 try
33 {
34 mNcFile->openFile( meshFile );
35 }
36 catch ( MDAL::Error &err )
37 {
38 err.setDriver( name() );
39 MDAL::Log::error( err );
40 return std::string();
41 }
42
43 std::vector<std::string> meshNames;
44
45 CFDimensions dims;
46 bool sqliteFileExist = check1DConnection( meshFile );
47 if ( sqliteFileExist )
48 {
49 populate1DMeshDimensions( dims );
50 if ( dims.size( CFDimensions::Vertex ) > 0 && dims.size( CFDimensions::Edge ) > 0 )
51 {
52 meshNames.push_back( "Mesh1D" );
53 }
54 }
55
56 populate2DMeshDimensions( dims );
57 if ( dims.size( CFDimensions::Face ) > 0 )
58 meshNames.push_back( "Mesh2D" );
59
60 if ( !meshNames.size() )
61 {
62 MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), "No meshes found in file" + meshFile );
63 return std::string( "" );
64 }
65
66 return MDAL::buildAndMergeMeshUris( meshFile, meshNames, name() );
67 }
68
populate2DMeshDimensions(MDAL::CFDimensions & dims)69 void MDAL::Driver3Di::populate2DMeshDimensions( MDAL::CFDimensions &dims )
70 {
71 size_t count;
72 int ncid;
73
74 // 2D Mesh
75 mNcFile->getDimension( "nMesh2D_nodes", &count, &ncid );
76 dims.setDimension( CFDimensions::Face, count, ncid );
77
78 mNcFile->getDimension( "nCorner_Nodes", &count, &ncid );
79 dims.setDimension( CFDimensions::MaxVerticesInFace, count, ncid );
80
81 // Vertices count is populated later in populateElements
82 // it is not known from the array dimensions
83 }
84
populateDimensions()85 MDAL::CFDimensions MDAL::Driver3Di::populateDimensions( )
86 {
87 CFDimensions dims;
88 size_t count;
89 int ncid;
90
91 if ( mRequestedMeshName == "Mesh1D" )
92 {
93 populate1DMeshDimensions( dims );
94 }
95 else
96 {
97 // Default loaded mesh is the 2D mesh
98 populate2DMeshDimensions( dims );
99 }
100
101 // Time
102 mNcFile->getDimension( "time", &count, &ncid );
103 dims.setDimension( CFDimensions::Time, count, ncid );
104
105 return dims;
106 }
107
populateElements(Vertices & vertices,Edges & edges,Faces & faces)108 void MDAL::Driver3Di::populateElements( Vertices &vertices, Edges &edges, Faces &faces )
109 {
110 if ( mRequestedMeshName == "Mesh1D" )
111 populateMesh1DElements( vertices, edges );
112 else
113 populateMesh2DElements( vertices, faces );
114 }
115
populateMesh2DElements(MDAL::Vertices & vertices,MDAL::Faces & faces)116 void MDAL::Driver3Di::populateMesh2DElements( MDAL::Vertices &vertices, MDAL::Faces &faces )
117 {
118 assert( vertices.empty() );
119 size_t faceCount = mDimensions.size( CFDimensions::Face );
120 faces.resize( faceCount );
121 size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
122 size_t arrsize = faceCount * verticesInFace;
123 std::map<std::string, size_t> xyToVertex2DId;
124
125 // X coordinate
126 int ncidX = mNcFile->getVarId( "Mesh2DContour_x" );
127 double fillX = mNcFile->getFillValue( ncidX );
128 std::vector<double> faceVerticesX( arrsize );
129 if ( nc_get_var_double( mNcFile->handle(), ncidX, faceVerticesX.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
130
131 // Y coordinate
132 int ncidY = mNcFile->getVarId( "Mesh2DContour_y" );
133 double fillY = mNcFile->getFillValue( ncidY );
134 std::vector<double> faceVerticesY( arrsize );
135 if ( nc_get_var_double( mNcFile->handle(), ncidY, faceVerticesY.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
136
137 // now populate create faces and backtrack which vertices
138 // are used in multiple faces
139 for ( size_t faceId = 0; faceId < faceCount; ++faceId )
140 {
141 Face face;
142
143 for ( size_t faceVertexId = 0; faceVertexId < verticesInFace; ++faceVertexId )
144 {
145 size_t arrId = faceId * verticesInFace + faceVertexId;
146 Vertex vertex;
147 vertex.x = faceVerticesX[arrId];
148 vertex.y = faceVerticesY[arrId];
149 vertex.z = 0;
150
151 if ( MDAL::equals( vertex.x, fillX ) || MDAL::equals( vertex.y, fillY ) )
152 break;
153
154
155 size_t vertexId;
156
157 std::string key = std::to_string( vertex.x ) + "," + std::to_string( vertex.y );
158 const auto it = xyToVertex2DId.find( key );
159 if ( it == xyToVertex2DId.end() )
160 {
161 // new vertex
162 vertexId = vertices.size();
163 xyToVertex2DId[key] = vertexId;
164 vertices.push_back( vertex );
165 }
166 else
167 {
168 // existing vertex
169 vertexId = it->second;
170 }
171
172 face.push_back( vertexId );
173
174 }
175
176 faces[faceId] = face;
177 }
178
179 // Only now we have number of vertices, since we identified vertices that
180 // are used in multiple faces
181 mDimensions.setDimension( CFDimensions::Vertex, vertices.size() );
182 }
183
addBedElevation(MemoryMesh * mesh)184 void MDAL::Driver3Di::addBedElevation( MemoryMesh *mesh )
185 {
186 assert( mesh );
187 if ( 0 == mesh->facesCount() )
188 return;
189
190 size_t faceCount = mesh->facesCount();
191
192 // read Z coordinate of 3di computation nodes centers
193 int ncidZ = mNcFile->getVarId( "Mesh2DFace_zcc" );
194 double fillZ = mNcFile->getFillValue( ncidZ );
195 std::vector<double> coordZ( faceCount );
196 if ( nc_get_var_double( mNcFile->handle(), ncidZ, coordZ.data() ) )
197 return; //error reading the array
198
199
200 std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >(
201 name(),
202 mesh,
203 mesh->uri(),
204 "Bed Elevation"
205 );
206
207 group->setDataLocation( MDAL_DataLocation::DataOnFaces );
208 group->setIsScalar( true );
209
210 std::shared_ptr<MDAL::MemoryDataset2D> dataset = std::make_shared< MemoryDataset2D >( group.get() );
211 dataset->setTime( MDAL::RelativeTimestamp() );
212 for ( size_t i = 0; i < faceCount; ++i )
213 {
214 dataset->setScalarValue( i, MDAL::safeValue( coordZ[i], fillZ ) );
215 }
216 dataset->setStatistics( MDAL::calculateStatistics( dataset ) );
217 group->setStatistics( MDAL::calculateStatistics( group ) );
218 group->datasets.push_back( dataset );
219 mesh->datasetGroups.push_back( group );
220 }
221
getCoordinateSystemVariableName()222 std::string MDAL::Driver3Di::getCoordinateSystemVariableName()
223 {
224 return "projected_coordinate_system";
225 }
226
getTimeVariableName() const227 std::string MDAL::Driver3Di::getTimeVariableName() const
228 {
229 return "time";
230 }
231
ignoreNetCDFVariables()232 std::set<std::string> MDAL::Driver3Di::ignoreNetCDFVariables()
233 {
234 std::set<std::string> ignore_variables;
235
236 ignore_variables.insert( "projected_coordinate_system" );
237 ignore_variables.insert( "time" );
238
239 std::vector<std::string> meshes;
240 meshes.push_back( "Mesh1D" );
241 meshes.push_back( "Mesh2D" );
242
243 for ( const std::string &mesh : meshes )
244 {
245 ignore_variables.insert( mesh );
246 ignore_variables.insert( mesh + "Node_id" );
247 ignore_variables.insert( mesh + "Node_type" );
248
249 ignore_variables.insert( mesh + "Face_xcc" );
250 ignore_variables.insert( mesh + "Face_ycc" );
251 ignore_variables.insert( mesh + "Face_zcc" );
252 ignore_variables.insert( mesh + "Contour_x" );
253 ignore_variables.insert( mesh + "Contour_y" );
254 ignore_variables.insert( mesh + "Face_sumax" );
255
256 ignore_variables.insert( mesh + "Line_id" );
257 ignore_variables.insert( mesh + "Line_xcc" );
258 ignore_variables.insert( mesh + "Line_ycc" );
259 ignore_variables.insert( mesh + "Line_zcc" );
260 ignore_variables.insert( mesh + "Line_type" );
261 }
262
263 return ignore_variables;
264 }
265
parseNetCDFVariableMetadata(int varid,std::string & variableName,std::string & name,bool * is_vector,bool * isPolar,bool * invertedDirection,bool * is_x)266 void MDAL::Driver3Di::parseNetCDFVariableMetadata( int varid,
267 std::string &variableName,
268 std::string &name,
269 bool *is_vector,
270 bool *isPolar,
271 bool *invertedDirection,
272 bool *is_x )
273 {
274 *is_vector = false;
275 *is_x = true;
276 *isPolar = false;
277 MDAL_UNUSED( invertedDirection )
278
279 std::string long_name = mNcFile->getAttrStr( "long_name", varid );
280 if ( long_name.empty() )
281 {
282 std::string standard_name = mNcFile->getAttrStr( "standard_name", varid );
283 if ( standard_name.empty() )
284 {
285 name = variableName;
286 }
287 else
288 {
289 variableName = standard_name;
290 if ( MDAL::contains( standard_name, "_x_" ) )
291 {
292 *is_vector = true;
293 name = MDAL::replace( standard_name, "_x_", "" );
294 }
295 else if ( MDAL::contains( standard_name, "_y_" ) )
296 {
297 *is_vector = true;
298 *is_x = false;
299 name = MDAL::replace( standard_name, "_y_", "" );
300 }
301 else
302 {
303 name = standard_name;
304 }
305 }
306 }
307 else
308 {
309 variableName = long_name;
310 if ( MDAL::contains( long_name, " in x direction" ) )
311 {
312 *is_vector = true;
313 name = MDAL::replace( long_name, " in x direction", "" );
314 }
315 else if ( MDAL::contains( long_name, " in y direction" ) )
316 {
317 *is_vector = true;
318 *is_x = false;
319 name = MDAL::replace( long_name, " in y direction", "" );
320 }
321 else
322 {
323 name = long_name;
324 }
325 }
326 }
327
parseClassification(int varid) const328 std::vector<std::pair<double, double>> MDAL::Driver3Di::parseClassification( int varid ) const
329 {
330 MDAL_UNUSED( varid );
331 return std::vector<std::pair<double, double>>();
332 }
333
populate1DMeshDimensions(MDAL::CFDimensions & dims)334 void MDAL::Driver3Di::populate1DMeshDimensions( MDAL::CFDimensions &dims )
335 {
336 size_t count;
337 int ncid;
338
339 mNcFile->getDimension( "nMesh1D_nodes", &count, &ncid );
340 dims.setDimension( CFDimensions::Vertex, count, ncid );
341
342 mNcFile->getDimension( "nMesh1D_lines", &count, &ncid );
343 dims.setDimension( CFDimensions::Edge, count, ncid );
344 }
345
populateMesh1DElements(MDAL::Vertices & vertices,MDAL::Edges & edges)346 void MDAL::Driver3Di::populateMesh1DElements( MDAL::Vertices &vertices, MDAL::Edges &edges )
347 {
348 assert( vertices.empty() && edges.empty() );
349 size_t vertexCount = mDimensions.size( CFDimensions::Vertex );
350 size_t edgesCount = mDimensions.size( CFDimensions::Edge );
351 vertices.resize( vertexCount );
352 edges.resize( edgesCount );
353
354 // X coordinate
355 int ncidX = mNcFile->getVarId( "Mesh1DNode_xcc" );
356 double fillX = mNcFile->getFillValue( ncidX );
357 std::vector<double> verticesX( vertexCount );
358 if ( nc_get_var_double( mNcFile->handle(), ncidX, verticesX.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
359
360 // Y coordinate
361 int ncidY = mNcFile->getVarId( "Mesh1DNode_ycc" );
362 double fillY = mNcFile->getFillValue( ncidY );
363 std::vector<double> verticesY( vertexCount );
364 if ( nc_get_var_double( mNcFile->handle(), ncidY, verticesY.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
365
366 // Z coordinate
367 int ncidZ = mNcFile->getVarId( "Mesh1DNode_zcc" );
368 double fillZ = mNcFile->getFillValue( ncidZ );
369 std::vector<double> verticesZ( vertexCount );
370 if ( nc_get_var_double( mNcFile->handle(), ncidZ, verticesZ.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
371
372 // Node Id
373 int ncidNodeId = mNcFile->getVarId( "Mesh1DNode_id" );
374 std::vector<int> verticesId( vertexCount );
375 if ( nc_get_var_int( mNcFile->handle(), ncidNodeId, verticesId.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
376
377 // Edge Id
378 int ncidEgeId = mNcFile->getVarId( "Mesh1DLine_id" );
379 std::vector<int> edgesId( edgesCount );
380 if ( nc_get_var_int( mNcFile->handle(), ncidEgeId, edgesId.data() ) ) throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unknown format" );
381
382 for ( size_t i = 0; i < vertexCount; ++i )
383 {
384 Vertex vertex;
385 vertex.x = verticesX[i];
386 if ( vertex.x == fillX )
387 vertex.x = std::numeric_limits<double>::quiet_NaN();
388 vertex.y = verticesY[i];
389 if ( vertex.y == fillY )
390 vertex.y = std::numeric_limits<double>::quiet_NaN();
391 vertex.z = verticesZ[i];
392 if ( vertex.z == fillZ )
393 vertex.z = std::numeric_limits<double>::quiet_NaN();
394 vertices[i] = vertex;
395 }
396
397 parse1DConnection( verticesId, edgesId, edges );
398 }
399
400
check1DConnection(std::string fileName)401 bool MDAL::Driver3Di::check1DConnection( std::string fileName )
402 {
403 std::string sqliteFile = dirName( fileName ) + "/gridadmin.sqlite";
404
405 if ( !fileExists( sqliteFile ) )
406 return false;
407
408 Sqlite3Db sqliteDb;
409 return sqliteDb.open( sqliteFile );
410
411 }
412
parse1DConnection(const std::vector<int> & nodesId,const std::vector<int> & edgesId,Edges & edges)413 void MDAL::Driver3Di::parse1DConnection( const std::vector<int> &nodesId,
414 const std::vector<int> &edgesId,
415 Edges &edges )
416 {
417 std::string sqliteFileName = dirName( mNcFile->getFileName() ) + "/gridadmin.sqlite";
418
419 //construct id map
420 std::map<int, size_t> edgeMap;
421 std::map<int, size_t> nodeMap;
422 for ( size_t edgeIndex = 0; edgeIndex < edges.size(); ++edgeIndex )
423 edgeMap[edgesId.at( edgeIndex )] = edgeIndex;
424
425 for ( size_t nodeIndex = 0; nodeIndex < nodesId.size(); ++nodeIndex )
426 nodeMap[nodesId.at( nodeIndex )] = nodeIndex;
427
428
429 Sqlite3Db db;
430 if ( !db.open( sqliteFileName ) || !( db.get() ) )
431 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unable to open sqlite database" );
432
433 Sqlite3Statement stmt;
434
435 if ( ! stmt.prepare( &db, "SELECT id, start_node_idx, end_node_idx FROM flowlines" ) )
436 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unable to read edges connectivity from sqlite database" );
437
438 if ( stmt.columnCount() < 0 || stmt.columnCount() != 3 )
439 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Invalid edges connectivity schema in sqlite database" );
440
441 while ( stmt.next() )
442 {
443 int idEdge, idStartNode, idEndNode;
444 idEdge = stmt.getInt( 0 );
445 idStartNode = stmt.getInt( 1 );
446 idEndNode = stmt.getInt( 2 );
447
448 auto itEdge = edgeMap.find( idEdge );
449 auto itStart = nodeMap.find( idStartNode );
450 auto itEnd = nodeMap.find( idEndNode );
451
452 if ( itEdge != edgeMap.end() && itStart != nodeMap.end() && itEnd != nodeMap.end() )
453 {
454 edges[itEdge->second].startVertex = itStart->second;
455 edges[itEdge->second].endVertex = itEnd->second;
456 }
457 }
458 }
459