1 /*
2 MDAL - Mesh Data Abstraction Library (MIT License)
3 Copyright (C) 2019 Peter Petrik (zilolv at gmail dot com)
4 */
5
6 #include "mdal_ugrid.hpp"
7 #include "mdal_utils.hpp"
8 #include "mdal_logger.hpp"
9
10 #include <netcdf.h>
11 #include <assert.h>
12 #include <algorithm>
13 #include <cmath>
14
15 #define FILL_COORDINATES_VALUE -999.0
16 #define FILL_FACE2D_VALUE -999
17
DriverUgrid()18 MDAL::DriverUgrid::DriverUgrid()
19 : DriverCF(
20 "Ugrid",
21 "UGRID",
22 "*.nc",
23 Capability::ReadMesh | Capability::SaveMesh )
24 {}
25
create()26 MDAL::DriverUgrid *MDAL::DriverUgrid::create()
27 {
28 return new DriverUgrid();
29 }
30
findMeshesNames() const31 std::vector<std::string> MDAL::DriverUgrid::findMeshesNames() const
32 {
33 std::vector<std::string> meshesInFile;
34
35 const std::vector<std::string> variables = mNcFile->readArrNames();
36 for ( const std::string &var : variables )
37 {
38 bool isMeshTopology = mNcFile->getAttrStr( var, "cf_role" ) == "mesh_topology";
39 if ( isMeshTopology )
40 {
41 // file can include more meshes
42 meshesInFile.push_back( var );
43 }
44 }
45
46 return meshesInFile;
47 }
48
buildUri(const std::string & meshFile)49 std::string MDAL::DriverUgrid::buildUri( const std::string &meshFile )
50 {
51 mNcFile.reset( new NetCDFFile );
52
53 try
54 {
55 mNcFile->openFile( meshFile );
56 }
57 catch ( MDAL::Error &err )
58 {
59 err.setDriver( name() );
60 MDAL::Log::error( err );
61 return std::string();
62 }
63
64 std::vector<std::string> meshNames = findMeshesNames();
65 if ( !meshNames.size() )
66 {
67 MDAL::Log::error( MDAL_Status::Err_UnknownFormat, name(), "No meshes found in file" + meshFile );
68 return std::string( "" );
69 }
70
71 // ignore network variable
72 std::vector<std::string>::iterator position = std::find( meshNames.begin(), meshNames.end(), "network" );
73 if ( position != meshNames.end() )
74 meshNames.erase( position );
75
76 return MDAL::buildAndMergeMeshUris( meshFile, meshNames, name() );
77 }
78
nodeZVariableName() const79 std::string MDAL::DriverUgrid::nodeZVariableName() const
80 {
81 const std::vector<std::string> variables = mNcFile->readArrNames();
82 for ( const std::string &varName : variables )
83 {
84 const std::string stdName = mNcFile->getAttrStr( varName, "standard_name" );
85 const std::string meshName = mNcFile->getAttrStr( varName, "mesh" );
86 const std::string location = mNcFile->getAttrStr( varName, "location" );
87
88 if ( stdName == "altitude" && meshName == mMeshName && location == "node" )
89 {
90 return varName;
91 }
92 }
93
94 // not found, the file in non UGRID standard conforming,
95 // but lets try the common name
96 return mMeshName + "_node_z";
97 }
98
populateDimensions()99 MDAL::CFDimensions MDAL::DriverUgrid::populateDimensions( )
100 {
101 CFDimensions dims;
102 size_t count;
103 int ncid;
104
105 mAllMeshNames = findMeshesNames();
106
107 if ( mAllMeshNames.empty() )
108 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, name(), "File " + mFileName + " does not contain any valid mesh definition" );
109
110 if ( !mRequestedMeshName.empty() )
111 {
112 if ( std::find( std::begin( mAllMeshNames ), std::end( mAllMeshNames ), mRequestedMeshName ) != std::end( mAllMeshNames ) )
113 mMeshName = mRequestedMeshName;
114 else
115 throw MDAL::Error( MDAL_Status::Err_InvalidData, "No such mesh with name: " + mRequestedMeshName, name() );
116 }
117 else
118 {
119 if ( mAllMeshNames.size() == 1 )
120 mMeshName = mAllMeshNames.at( 0 );
121 else // there are more meshes in file
122 {
123 if ( MDAL::contains( mAllMeshNames.at( 0 ), "network" ) ) // ignore the network variable for a moment
124 mMeshName = mAllMeshNames.at( 1 );
125 else
126 mMeshName = mAllMeshNames.at( 0 );
127
128 MDAL::Log::warning( MDAL_Status::Warn_MultipleMeshesInFile, name(), "Found multiple meshes in file, working with: " + mMeshName );
129 }
130 }
131
132 if ( mMeshName.empty() ) throw MDAL::Error( MDAL_Status::Err_InvalidData, "Unable to parse mesh name from file" );
133
134 mMeshDimension = mNcFile->getAttrInt( mMeshName, "topology_dimension" );
135
136 if ( ( mMeshDimension < 1 ) || ( mMeshDimension > 2 ) )
137 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, name(), "Unable to parse topology dimension from mesh or mesh is 3D" );
138
139 MDAL::Log::info( "Parsing " + std::to_string( mMeshDimension ) + "D mesh with name: " + mMeshName );
140
141 std::string nodeXVariable, nodeYVariable;
142 if ( mMeshDimension == 1 )
143 parseCoordinatesFrom1DMesh( mMeshName, "node_coordinates", nodeXVariable, nodeYVariable );
144 else
145 parse2VariablesFromAttribute( mMeshName, "node_coordinates", nodeXVariable, nodeYVariable, false );
146
147 std::vector<size_t> nodeDimension;
148 std::vector<int> nodeDimensionId;
149 mNcFile->getDimensions( nodeXVariable, nodeDimension, nodeDimensionId );
150
151 if ( nodeDimension.size() != 1 )
152 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Error while parsing dimensions" );
153
154 dims.setDimension( CFDimensions::Vertex, nodeDimension.at( 0 ), nodeDimensionId.at( 0 ) );
155
156 /* continue parsing dimension dependent variables */
157 if ( mMeshDimension == 1 )
158 populate1DMeshDimensions( dims );
159 else
160 populate2DMeshDimensions( dims, ncid );
161
162 /* Time variable - not required for UGRID format */
163 if ( mNcFile->hasDimension( "time" ) )
164 {
165 mNcFile->getDimension( "time", &count, &ncid );
166 dims.setDimension( CFDimensions::Time, count, ncid );
167 }
168 else
169 {
170 dims.setDimension( CFDimensions::Time, 0, -1 );
171 }
172
173 return dims;
174 }
175
populate1DMeshDimensions(MDAL::CFDimensions & dims)176 void MDAL::DriverUgrid::populate1DMeshDimensions( MDAL::CFDimensions &dims )
177 {
178 /* Parse number of edges ( dimension ) from mesh */
179 std::string edgeConnectivityVariableName = mNcFile->getAttrStr( mMeshName, "edge_node_connectivity" );
180 if ( edgeConnectivityVariableName.empty() )
181 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Did not find edge node connectivity attribute" );
182
183 std::vector<size_t> edgeDimension;
184 std::vector<int> edgeDimensionId;
185 mNcFile->getDimensions( edgeConnectivityVariableName, edgeDimension, edgeDimensionId );
186 if ( edgeDimension.size() != 2 )
187 throw MDAL::Error( MDAL_Status::Err_InvalidData, name(), "Unable to parse dimensions for edge_nodes_connectivity variable" );
188
189 size_t edgesCount = edgeDimension.at( 0 ); // Only interested in first value, edge will always have only 2 nodes
190 int edgesCountId = edgeDimensionId.at( 0 );
191
192 dims.setDimension( CFDimensions::Edge, edgesCount, edgesCountId );
193 }
194
populate2DMeshDimensions(MDAL::CFDimensions & dims,int & ncid)195 void MDAL::DriverUgrid::populate2DMeshDimensions( MDAL::CFDimensions &dims, int &ncid )
196 {
197 // face dimension location is retrieved from the face_node_connectivity variable
198 // if face_dimension is defined as attribute, the dimension at this location help to desambiguate vertex per faces and number of faces
199 std::string faceConnectivityVariablesName = mNcFile->getAttrStr( mMeshName, "face_node_connectivity" );
200 std::string faceDimensionLocation = mNcFile->getAttrStr( mMeshName, "face_dimension" );
201 if ( faceConnectivityVariablesName == "" )
202 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Did not find face connectivity attribute" );
203
204 size_t facesCount;
205 size_t maxVerticesPerFace;
206 size_t count;
207
208 std::vector<size_t> faceDimension;
209 std::vector<int> faceDimensionId;
210 int facesIndexDimensionId;
211 int maxVerticesPerFaceDimensionId;
212 mNcFile->getDimensions( faceConnectivityVariablesName, faceDimension, faceDimensionId );
213 if ( faceDimension.size() != 2 )
214 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Face dimension is 2D" );
215
216 // if face_dimension is not present in file, get it from dimension element
217 if ( faceDimensionLocation != "" )
218 {
219 mNcFile->getDimension( faceDimensionLocation, &facesCount, &ncid );
220 if ( facesCount == faceDimension.at( 0 ) )
221 {
222 facesIndexDimensionId = faceDimensionId.at( 0 );
223 maxVerticesPerFaceDimensionId = faceDimensionId.at( 1 );
224 maxVerticesPerFace = faceDimension.at( 1 );
225 }
226 else
227 {
228 facesIndexDimensionId = faceDimensionId.at( 1 );
229 maxVerticesPerFaceDimensionId = faceDimensionId.at( 0 );
230 maxVerticesPerFace = faceDimension.at( 0 );
231 }
232 }
233 else
234 {
235 facesIndexDimensionId = faceDimensionId.at( 0 );
236 facesCount = faceDimension.at( 0 );
237 maxVerticesPerFaceDimensionId = faceDimensionId.at( 1 );
238 maxVerticesPerFace = faceDimension.at( 1 );
239 }
240
241 dims.setDimension( CFDimensions::Face, facesCount, facesIndexDimensionId );
242 dims.setDimension( CFDimensions::MaxVerticesInFace, maxVerticesPerFace, maxVerticesPerFaceDimensionId );
243
244 // number of edges in the mesh, not required for UGRID format
245 const std::string mesh2dEdge = mNcFile->getAttrStr( mMeshName, "edge_dimension" );
246 if ( mNcFile->hasDimension( mesh2dEdge ) )
247 {
248 mNcFile->getDimension( mesh2dEdge, &count, &ncid );
249 dims.setDimension( CFDimensions::Face2DEdge, count, ncid );
250 }
251 else
252 {
253 dims.setDimension( CFDimensions::Face2DEdge, 0, -1 );
254 }
255 }
256
populateElements(Vertices & vertices,Edges & edges,Faces & faces)257 void MDAL::DriverUgrid::populateElements( Vertices &vertices, Edges &edges, Faces &faces )
258 {
259 populateVertices( vertices );
260
261 if ( mMeshDimension == 1 )
262 populateEdges( edges ); // 1D mesh
263 else
264 populateFaces( faces ); // 2D mesh
265 }
266
populateVertices(MDAL::Vertices & vertices)267 void MDAL::DriverUgrid::populateVertices( MDAL::Vertices &vertices )
268 {
269 assert( vertices.empty() );
270 size_t vertexCount = mDimensions.size( CFDimensions::Vertex );
271 vertices.resize( vertexCount );
272 Vertex *vertexPtr = vertices.data();
273
274 // node_coordinates should be something like Mesh2D_node_x Mesh2D_node_y
275 std::string verticesXName, verticesYName;
276 if ( mMeshDimension == 1 )
277 parseCoordinatesFrom1DMesh( mMeshName, "node_coordinates", verticesXName, verticesYName );
278 else
279 parse2VariablesFromAttribute( mMeshName, "node_coordinates", verticesXName, verticesYName, false );
280
281 const std::vector<double> verticesX = mNcFile->readDoubleArr( verticesXName, vertexCount );
282 const std::vector<double> verticesY = mNcFile->readDoubleArr( verticesYName, vertexCount );
283
284 std::vector<double> verticesZ;
285 if ( mNcFile->hasArr( nodeZVariableName() ) )
286 {
287 verticesZ = mNcFile->readDoubleArr( nodeZVariableName(), vertexCount );
288 }
289
290 if ( verticesX.size() == 1 &&
291 verticesY.size() == 1 &&
292 verticesZ.size() == 1 &&
293 verticesX.at( 0 ) == FILL_COORDINATES_VALUE &&
294 verticesY.at( 0 ) == FILL_COORDINATES_VALUE &&
295 verticesZ.at( 0 ) == FILL_COORDINATES_VALUE )
296 {
297 vertexCount = 0;
298 vertices.clear();
299 }
300
301
302 for ( size_t i = 0; i < vertexCount; ++i, ++vertexPtr )
303 {
304 vertexPtr->x = verticesX[i];
305 vertexPtr->y = verticesY[i];
306 if ( !verticesZ.empty() )
307 vertexPtr->z = verticesZ[i];
308 }
309 }
310
populateEdges(MDAL::Edges & edges)311 void MDAL::DriverUgrid::populateEdges( MDAL::Edges &edges )
312 {
313 assert( edges.empty() );
314
315 // number of edges
316 size_t edgesCount = mDimensions.size( CFDimensions::Edge );
317 edges.resize( edgesCount );
318
319 const std::string edgeNodeConnectivityVar = mNcFile->getAttrStr( mMeshName, "edge_node_connectivity" );
320 if ( edgeNodeConnectivityVar.empty() )
321 MDAL::Log::error( MDAL_Status::Err_MissingDriver, "Unable to find edge_node_connectivity attribute of " + mMeshName );
322
323 // load edges
324 std::vector<int> edgeNodesIdxs = mNcFile->readIntArr( edgeNodeConnectivityVar, edgesCount * 2 ); // two nodes per edge
325 int startIndex = mNcFile->getAttrInt( edgeNodeConnectivityVar, "start_index" );
326
327 // iterate over all edge_nodes coordinates - those are indexes for nodes
328 for ( size_t i = 0; i < edgesCount; ++i )
329 {
330
331 int startEdgeIx = MDAL::toInt( i ) * 2;
332 int endEdgeIx = MDAL::toInt( i ) * 2 + 1;
333
334 edges[i].startVertex = edgeNodesIdxs[startEdgeIx] - startIndex;
335 edges[i].endVertex = edgeNodesIdxs[endEdgeIx] - startIndex;
336 }
337 }
338
populateFaces(MDAL::Faces & faces)339 void MDAL::DriverUgrid::populateFaces( MDAL::Faces &faces )
340 {
341 assert( faces.empty() );
342 size_t faceCount = mDimensions.size( CFDimensions::Face );
343 faces.resize( faceCount );
344
345 // Parse 2D Mesh
346 // face_node_connectivity is usually something like Mesh2D_face_nodes
347 const std::string mesh2dFaceNodeConnectivity = mNcFile->getAttrStr( mMeshName, "face_node_connectivity" );
348
349 size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
350 int fillVal = -1;
351 if ( mNcFile->hasAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" ) )
352 fillVal = mNcFile->getAttrInt( mesh2dFaceNodeConnectivity, "_FillValue" );
353 int startIndex = mNcFile->getAttrInt( mesh2dFaceNodeConnectivity, "start_index" );
354 std::vector<int> faceNodesConn = mNcFile->readIntArr( mesh2dFaceNodeConnectivity, faceCount * verticesInFace );
355
356 for ( size_t i = 0; i < faceCount; ++i )
357 {
358 std::vector<size_t> idxs;
359
360 for ( size_t j = 0; j < verticesInFace; ++j )
361 {
362 size_t idx = verticesInFace * i + j;
363 int val = faceNodesConn[idx];
364
365 if ( fillVal == val )
366 {
367 // found fill val
368 break;
369 }
370 else
371 {
372 idxs.push_back( static_cast<size_t>( val - startIndex ) );
373 }
374 }
375 faces[i] = idxs;
376 }
377
378 if ( faces.size() == 1 && faces.at( 0 ).size() == 0 )
379 faces.clear();
380 }
381
addBedElevation(MDAL::MemoryMesh * mesh)382 void MDAL::DriverUgrid::addBedElevation( MDAL::MemoryMesh *mesh )
383 {
384 if ( mNcFile->hasArr( nodeZVariableName() ) ) MDAL::addBedElevationDatasetGroup( mesh, mesh->vertices() );
385 }
386
getCoordinateSystemVariableName()387 std::string MDAL::DriverUgrid::getCoordinateSystemVariableName()
388 {
389 std::string coordinateSystemVariable;
390
391 // first try to get the coordinate system variable from grid definition
392 std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile->getAttrStr( mMeshName, "node_coordinates" ), ' ' );
393 if ( nodeVariablesName.size() > 1 )
394 {
395 if ( mNcFile->hasArr( nodeVariablesName[0] ) )
396 {
397 coordinateSystemVariable = mNcFile->getAttrStr( nodeVariablesName[0], "grid_mapping" );
398 }
399 }
400
401
402 // if automatic discovery fails, try to check some hardcoded common variables that store projection
403 if ( coordinateSystemVariable.empty() )
404 {
405 if ( mNcFile->hasArr( "projected_coordinate_system" ) )
406 coordinateSystemVariable = "projected_coordinate_system";
407 else if ( mNcFile->hasArr( "wgs84" ) )
408 coordinateSystemVariable = "wgs84";
409 }
410
411 // return, may be empty
412 return coordinateSystemVariable;
413 }
414
ignoreNetCDFVariables()415 std::set<std::string> MDAL::DriverUgrid::ignoreNetCDFVariables()
416 {
417 std::set<std::string> ignoreVariables;
418
419 ignoreVariables.insert( "projected_coordinate_system" );
420 ignoreVariables.insert( "time" );
421 ignoreVariables.insert( "timestep" );
422
423 for ( const std::string &mesh : mAllMeshNames )
424 {
425 ignoreVariables.insert( mesh );
426
427 int dim = mNcFile->getAttrInt( mesh, "topology_dimension" );
428 if ( dim == 1 )
429 ignore1DMeshVariables( mesh, ignoreVariables );
430 else
431 ignore2DMeshVariables( mesh, ignoreVariables );
432 }
433
434 return ignoreVariables;
435 }
436
ignore1DMeshVariables(const std::string & mesh,std::set<std::string> & ignoreVariables)437 void MDAL::DriverUgrid::ignore1DMeshVariables( const std::string &mesh, std::set<std::string> &ignoreVariables )
438 {
439 // ignore all variables with network in name
440 // network topology does not contain any data
441 if ( MDAL::contains( mesh, "network" ) )
442 {
443 std::vector<std::string> variables = mNcFile->readArrNames();
444 for ( const std::string &var : variables )
445 {
446 if ( MDAL::contains( var, "network" ) )
447 ignoreVariables.insert( var );
448 }
449 return;
450 }
451
452 ignoreVariables.insert( mNcFile->getAttrStr( mesh, "edge_node_connectivity" ) );
453 ignoreVariables.insert( mNcFile->getAttrStr( mesh, "node_id" ) );
454 ignoreVariables.insert( mNcFile->getAttrStr( mesh, "node_long_name" ) );
455
456 std::vector<std::string> coordinateVarsToIgnore {"node_coordinates", "edge_coordinates"};
457
458 for ( const std::string &coordinateIt : coordinateVarsToIgnore )
459 {
460 std::string coordinatesVar = mNcFile->getAttrStr( mesh, coordinateIt );
461 std::vector<std::string> allCoords = MDAL::split( coordinatesVar, " " );
462
463 for ( const std::string &var : allCoords )
464 {
465 ignoreVariables.insert( var );
466 ignoreVariables.insert( mNcFile->getAttrStr( var, "bounds" ) );
467 }
468 }
469 }
470
ignore2DMeshVariables(const std::string & mesh,std::set<std::string> & ignoreVariables)471 void MDAL::DriverUgrid::ignore2DMeshVariables( const std::string &mesh, std::set<std::string> &ignoreVariables )
472 {
473 std::string xName, yName;
474 parse2VariablesFromAttribute( mesh, "node_coordinates", xName, yName, true );
475 ignoreVariables.insert( xName );
476 ignoreVariables.insert( yName );
477 ignoreVariables.insert( nodeZVariableName() );
478 ignoreVariables.insert( mNcFile->getAttrStr( mesh, "edge_node_connectivity" ) );
479 parse2VariablesFromAttribute( mesh, "edge_coordinates", xName, yName, true );
480
481 if ( !xName.empty() )
482 {
483 ignoreVariables.insert( xName );
484 ignoreVariables.insert( mNcFile->getAttrStr( xName, "bounds" ) );
485 }
486 if ( !yName.empty() )
487 {
488 ignoreVariables.insert( yName );
489 ignoreVariables.insert( mNcFile->getAttrStr( yName, "bounds" ) );
490 }
491
492 ignoreVariables.insert( mNcFile->getAttrStr( mesh, "face_node_connectivity" ) );
493 parse2VariablesFromAttribute( mesh, "face_coordinates", xName, yName, true );
494 if ( !xName.empty() )
495 {
496 ignoreVariables.insert( xName );
497 ignoreVariables.insert( mNcFile->getAttrStr( xName, "bounds" ) );
498 }
499 if ( !yName.empty() )
500 {
501 ignoreVariables.insert( yName );
502 ignoreVariables.insert( mNcFile->getAttrStr( yName, "bounds" ) );
503 }
504 ignoreVariables.insert( mNcFile->getAttrStr( mesh, "edge_face_connectivity" ) );
505 }
506
parseNetCDFVariableMetadata(int varid,std::string & variableName,std::string & name,bool * isVector,bool * isPolar,bool * invertedDirection,bool * isX)507 void MDAL::DriverUgrid::parseNetCDFVariableMetadata( int varid,
508 std::string &variableName,
509 std::string &name,
510 bool *isVector,
511 bool *isPolar,
512 bool *invertedDirection,
513 bool *isX )
514 {
515 *isVector = false;
516 *isX = true;
517 *isPolar = false;
518 *invertedDirection = false;
519
520 std::string longName = mNcFile->getAttrStr( "long_name", varid );
521 if ( longName.empty() )
522 {
523 std::string standardName = mNcFile->getAttrStr( "standard_name", varid );
524 if ( standardName.empty() )
525 {
526 name = variableName;
527 }
528 else
529 {
530 variableName = standardName;
531 if ( MDAL::contains( standardName, "_x_" ) )
532 {
533 *isVector = true;
534 name = MDAL::replace( standardName, "_x_", "" );
535 }
536 else if ( MDAL::contains( standardName, "_y_" ) )
537 {
538 *isVector = true;
539 *isX = false;
540 name = MDAL::replace( standardName, "_y_", "" );
541 }
542 else if ( MDAL::contains( standardName, "_from_direction" ) )
543 {
544 *isVector = true;
545 *isPolar = true;
546 *isX = false;
547 *invertedDirection = true;
548 name = MDAL::replace( standardName, "_speed", "_velocity" );
549 name = MDAL::replace( name, "_from_direction", "" );
550 }
551 else if ( MDAL::contains( standardName, "_to_direction" ) )
552 {
553 *isVector = true;
554 *isPolar = true;
555 *isX = false;
556 name = MDAL::replace( standardName, "_speed", "_velocity" );
557 name = MDAL::replace( name, "_to_direction", "" );
558 }
559 else
560 {
561 name = standardName;
562 }
563 }
564 }
565 else
566 {
567 variableName = longName;
568 if ( MDAL::contains( longName, ", x-component" ) || MDAL::contains( longName, "u component of " ) )
569 {
570 *isVector = true;
571 name = MDAL::replace( longName, ", x-component", "" );
572 name = MDAL::replace( name, "u component of ", "" );
573 }
574 else if ( MDAL::contains( longName, ", y-component" ) || MDAL::contains( longName, "v component of " ) )
575 {
576 *isVector = true;
577 *isX = false;
578 name = MDAL::replace( longName, ", y-component", "" );
579 name = MDAL::replace( name, "v component of ", "" );
580 }
581 else if ( MDAL::contains( longName, " magnitude" ) )
582 {
583 *isVector = true;
584 *isPolar = true;
585 *isX = true;
586 name = MDAL::replace( longName, "speed", "velocity" );
587 name = MDAL::removeFrom( name, " magnitude" );
588 }
589 else if ( MDAL::contains( longName, "direction" ) )
590 {
591 *isVector = true;
592 *isPolar = true;
593 *isX = false;
594
595 // check from_/to_direction in standard_name
596 std::string standardName = mNcFile->getAttrStr( "standard_name", varid );
597 *invertedDirection = MDAL::contains( longName, "from direction" );
598
599 name = MDAL::replace( longName, "speed", "velocity" );
600 name = MDAL::removeFrom( name, " from direction" );
601 name = MDAL::removeFrom( name, " to direction" );
602 name = MDAL::removeFrom( name, " direction" );
603 }
604 else
605 {
606 name = longName;
607 }
608 }
609 }
610
getTimeVariableName() const611 std::string MDAL::DriverUgrid::getTimeVariableName() const
612 {
613 return "time";
614 }
615
parseCoordinatesFrom1DMesh(const std::string & meshName,const std::string & attr_name,std::string & var1,std::string & var2)616 void MDAL::DriverUgrid::parseCoordinatesFrom1DMesh( const std::string &meshName, const std::string &attr_name, std::string &var1, std::string &var2 )
617 {
618 std::vector<std::string> nodeVariablesName = MDAL::split( mNcFile->getAttrStr( meshName, attr_name ), ' ' );
619
620 if ( nodeVariablesName.size() < 2 )
621 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Error while parsing node coordinates" );
622 else if ( nodeVariablesName.size() > 3 )
623 {
624 // format does not follow UGRID convention and have extra variables in coordinate attribute
625 // tring to parse coordinate variables, it usually ends with _x and _y, e.g. mesh1d_node_x, mesh1d_node_y
626
627 MDAL::Log::warning( MDAL_Status::Warn_InvalidElements, name(),
628 "Node coordinates consists of more than 3 variables, taking variable with _x in name by default" );
629
630 for ( const auto &nodeVar : nodeVariablesName )
631 {
632 if ( MDAL::contains( nodeVar, "_x" ) )
633 {
634 var1 = nodeVar;
635 }
636 else if ( MDAL::contains( nodeVar, "_y" ) )
637 {
638 var2 = nodeVar;
639 }
640 }
641 if ( var1.empty() || var2.empty() )
642 throw MDAL::Error( MDAL_Status::Err_InvalidData, name(), "Could not parse node coordinates from mesh" );
643 }
644 else // 2 variables as node coordinates
645 {
646 var1 = nodeVariablesName.at( 0 );
647 var2 = nodeVariablesName.at( 1 );
648 }
649 }
650
parse2VariablesFromAttribute(const std::string & name,const std::string & attr_name,std::string & var1,std::string & var2,bool optional) const651 void MDAL::DriverUgrid::parse2VariablesFromAttribute( const std::string &name, const std::string &attr_name,
652 std::string &var1, std::string &var2, bool optional ) const
653 {
654 const std::string mesh2dNodeCoordinates = mNcFile->getAttrStr( name, attr_name );
655 const std::vector<std::string> chunks = MDAL::split( mesh2dNodeCoordinates, ' ' );
656
657 if ( chunks.size() != 2 )
658 {
659 if ( optional )
660 {
661 var1 = "";
662 var2 = "";
663 }
664 else
665 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Unable to parse variables from attribute" );
666 }
667 else
668 {
669 var1 = chunks[0];
670 var2 = chunks[1];
671 }
672 }
673
save(const std::string & fileName,const std::string &,MDAL::Mesh * mesh)674 void MDAL::DriverUgrid::save( const std::string &fileName, const std::string &, MDAL::Mesh *mesh )
675 {
676 mFileName = fileName;
677
678 try
679 {
680 // Create file
681 mNcFile.reset( new NetCDFFile );
682 mNcFile->createFile( mFileName );
683
684 // Write globals
685 writeGlobals( );
686
687 // Write variables
688 writeVariables( mesh );
689 }
690 catch ( MDAL_Status error )
691 {
692 MDAL::Log::error( error, name(), "could not save file " + fileName );
693 }
694 catch ( MDAL::Error err )
695 {
696 MDAL::Log::error( err, name() );
697 }
698 }
699
saveMeshOnFileSuffix() const700 std::string MDAL::DriverUgrid::saveMeshOnFileSuffix() const
701 {
702 return "nc";
703 }
704
writeVariables(MDAL::Mesh * mesh)705 void MDAL::DriverUgrid::writeVariables( MDAL::Mesh *mesh )
706 {
707 // Global dimensions
708 ;
709 int dimNodeCountId = mNcFile->defineDimension( "nmesh2d_node", mesh->verticesCount() == 0 ? 1 : mesh->verticesCount() ); //if no vertices, set 1 since 0==NC_UNLIMITED
710 int dimFaceCountId = mNcFile->defineDimension( "nmesh2d_face", mesh->facesCount() == 0 ? 1 : mesh->facesCount() ); //if no vertices, set 1 since 0==NC_UNLIMITED
711 mNcFile->defineDimension( "nmesh2d_edge", 1 ); // no data on edges, cannot be 0, since 0==NC_UNLIMITED
712 int dimTimeId = mNcFile->defineDimension( "time", NC_UNLIMITED );
713 int dimMaxNodesPerFaceId = mNcFile->defineDimension( "max_nmesh2d_face_nodes",
714 mesh->faceVerticesMaximumCount() == 0 ? 1 : mesh->faceVerticesMaximumCount() ); //if 0, set 1 since 0==NC_UNLIMITED
715
716 // Mesh 2D Definition
717 int mesh2dId = mNcFile->defineVar( "mesh2d", NC_INT, 0, nullptr );
718 mNcFile->putAttrStr( mesh2dId, "cf_role", "mesh_topology" );
719 mNcFile->putAttrStr( mesh2dId, "long_name", "Topology data of 2D network" );
720 mNcFile->putAttrInt( mesh2dId, "topology_dimension", 2 );
721 mNcFile->putAttrStr( mesh2dId, "node_coordinates", "mesh2d_node_x mesh2d_node_y" );
722 mNcFile->putAttrStr( mesh2dId, "node_dimension", "nmesh2d_node" );
723 mNcFile->putAttrStr( mesh2dId, "edge_dimension", "nmesh2d_edge" );
724 mNcFile->putAttrStr( mesh2dId, "max_face_nodes_dimension", "max_nmesh2d_face_nodes" );
725 mNcFile->putAttrStr( mesh2dId, "face_node_connectivity", "mesh2d_face_nodes" );
726 mNcFile->putAttrStr( mesh2dId, "face_dimension", "nmesh2d_face" );
727
728 // Nodes X coordinate
729 int mesh2dNodeXId = mNcFile->defineVar( "mesh2d_node_x", NC_DOUBLE, 1, &dimNodeCountId );
730 mNcFile->putAttrStr( mesh2dNodeXId, "standard_name", "projection_x_coordinate" );
731 mNcFile->putAttrStr( mesh2dNodeXId, "long_name", "x-coordinate of mesh nodes" );
732 mNcFile->putAttrStr( mesh2dNodeXId, "mesh", "mesh2d" );
733 mNcFile->putAttrStr( mesh2dNodeXId, "location", "node" );
734
735 // Nodes Y coordinate
736 int mesh2dNodeYId = mNcFile->defineVar( "mesh2d_node_y", NC_DOUBLE, 1, &dimNodeCountId );
737 mNcFile->putAttrStr( mesh2dNodeYId, "standard_name", "projection_y_coordinate" );
738 mNcFile->putAttrStr( mesh2dNodeYId, "long_name", "y-coordinate of mesh nodes" );
739 mNcFile->putAttrStr( mesh2dNodeYId, "mesh", "mesh2d" );
740 mNcFile->putAttrStr( mesh2dNodeYId, "location", "node" );
741
742 // Nodes Z coordinate
743 int mesh2dNodeZId = mNcFile->defineVar( "mesh2d_node_z", NC_DOUBLE, 1, &dimNodeCountId );
744 mNcFile->putAttrStr( mesh2dNodeZId, "mesh", "mesh2d" );
745 mNcFile->putAttrStr( mesh2dNodeZId, "location", "node" );
746 mNcFile->putAttrStr( mesh2dNodeZId, "coordinates", "mesh2d_node_x mesh2d_node_y" );
747 mNcFile->putAttrStr( mesh2dNodeZId, "standard_name", "altitude" );
748 mNcFile->putAttrStr( mesh2dNodeZId, "long_name", "z-coordinate of mesh nodes" );
749 mNcFile->putAttrStr( mesh2dNodeZId, "grid_mapping", "projected_coordinate_system" );
750 mNcFile->putAttrDouble( mesh2dNodeZId, "_FillValue", FILL_COORDINATES_VALUE );
751
752 // Faces 2D Variable
753 int mesh2FaceNodesId_dimIds [] { dimFaceCountId, dimMaxNodesPerFaceId };
754 int mesh2FaceNodesId = mNcFile->defineVar( "mesh2d_face_nodes", NC_INT, 2, mesh2FaceNodesId_dimIds );
755 mNcFile->putAttrStr( mesh2FaceNodesId, "cf_role", "face_node_connectivity" );
756 mNcFile->putAttrStr( mesh2FaceNodesId, "mesh", "mesh2d" );
757 mNcFile->putAttrStr( mesh2FaceNodesId, "location", "face" );
758 mNcFile->putAttrStr( mesh2FaceNodesId, "long_name", "Mapping from every face to its corner nodes (counterclockwise)" );
759 mNcFile->putAttrInt( mesh2FaceNodesId, "start_index", 0 );
760 mNcFile->putAttrInt( mesh2FaceNodesId, "_FillValue", FILL_FACE2D_VALUE );
761
762 // Projected Coordinate System
763 int pcsId = mNcFile->defineVar( "projected_coordinate_system", NC_INT, 0, nullptr );
764
765 if ( mesh->crs() == "" )
766 {
767 mNcFile->putAttrInt( pcsId, "epsg", 0 );
768 mNcFile->putAttrStr( pcsId, "EPSG_code", "epgs:0" );
769 }
770 else
771 {
772 std::vector<std::string> words = MDAL::split( mesh->crs(), ":" );
773
774 if ( words[0] == "EPSG" && words.size() > 1 )
775 {
776 mNcFile->putAttrInt( pcsId, "epsg", std::stoi( words[1] ) );
777 mNcFile->putAttrStr( pcsId, "EPSG_code", mesh->crs() );
778 }
779 else
780 {
781 mNcFile->putAttrStr( pcsId, "wkt", mesh->crs() );
782 }
783 }
784
785 // Time array
786 int timeId = mNcFile->defineVar( "time", NC_DOUBLE, 1, &dimTimeId );
787 mNcFile->putAttrStr( timeId, "units", "hours since 2000-01-01 00:00:00" );
788
789 // Turning off define mode - allows data write
790 nc_enddef( mNcFile->handle() );
791
792 // Write vertices
793
794 const size_t maxBufferSize = 1000;
795 const size_t bufferSize = std::min( mesh->verticesCount(), maxBufferSize );
796 const size_t verticesCoordCount = bufferSize * 3;
797
798 std::vector<double> verticesCoordinates( verticesCoordCount );
799 std::unique_ptr<MDAL::MeshVertexIterator> vertexIterator = mesh->readVertices();
800
801 if ( mesh->verticesCount() == 0 )
802 {
803 // if there is no vertices fill the first fake vertex, see global dimension
804 mNcFile->putDataDouble( mesh2dNodeXId, 0, FILL_COORDINATES_VALUE );
805 mNcFile->putDataDouble( mesh2dNodeYId, 0, FILL_COORDINATES_VALUE );
806 mNcFile->putDataDouble( mesh2dNodeZId, 0, FILL_COORDINATES_VALUE );
807 }
808 else
809 {
810 size_t vertexIndex = 0;
811 size_t vertexFileIndex = 0;
812 while ( vertexIndex < mesh->verticesCount() )
813 {
814 size_t verticesRead = vertexIterator->next( bufferSize, verticesCoordinates.data() );
815 if ( verticesRead == 0 )
816 break;
817
818 for ( size_t i = 0; i < verticesRead; i++ )
819 {
820 mNcFile->putDataDouble( mesh2dNodeXId, vertexFileIndex, verticesCoordinates[3 * i] );
821 mNcFile->putDataDouble( mesh2dNodeYId, vertexFileIndex, verticesCoordinates[3 * i + 1] );
822 if ( std::isnan( verticesCoordinates[3 * i + 2] ) )
823 mNcFile->putDataDouble( mesh2dNodeZId, vertexFileIndex, FILL_COORDINATES_VALUE );
824 else
825 mNcFile->putDataDouble( mesh2dNodeZId, vertexFileIndex, verticesCoordinates[3 * i + 2] );
826 vertexFileIndex++;
827 }
828 vertexIndex += verticesRead;
829 }
830 }
831
832 // Write faces
833 std::unique_ptr<MDAL::MeshFaceIterator> faceIterator = mesh->readFaces();
834 const size_t faceVerticesMax = mesh->faceVerticesMaximumCount();
835 const size_t facesCount = mesh->facesCount();
836 const size_t faceOffsetsBufferLen = std::min( facesCount, maxBufferSize );
837 const size_t vertexIndicesBufferLen = faceOffsetsBufferLen * faceVerticesMax;
838
839 std::vector<int> faceOffsetsBuffer( faceOffsetsBufferLen );
840 std::vector<int> vertexIndicesBuffer( vertexIndicesBufferLen );
841
842 size_t faceIndex = 0;
843
844 if ( facesCount == 0 )
845 {
846 // if there is no vertices fill the first fake vertex, see global dimension
847 int fillValue = FILL_FACE2D_VALUE;
848 mNcFile->putDataArrayInt( 0, 0, 1, &fillValue );
849 }
850 else
851 {
852 while ( faceIndex < facesCount )
853 {
854 size_t facesRead = faceIterator->next(
855 faceOffsetsBufferLen,
856 faceOffsetsBuffer.data(),
857 vertexIndicesBufferLen,
858 vertexIndicesBuffer.data() );
859 if ( facesRead == 0 )
860 break;
861
862 for ( size_t i = 0; i < facesRead; i++ )
863 {
864 std::vector<int> verticesFaceData( faceVerticesMax, FILL_FACE2D_VALUE );
865 int startIndex = 0;
866 if ( i > 0 )
867 startIndex = faceOffsetsBuffer[ i - 1 ];
868 int endIndex = faceOffsetsBuffer[ i ];
869
870 size_t k = 0;
871 for ( int j = startIndex; j < endIndex; ++j )
872 {
873 int vertexIndex = vertexIndicesBuffer[ static_cast<size_t>( j ) ];
874 verticesFaceData[k++] = vertexIndex;
875 }
876 mNcFile->putDataArrayInt( mesh2FaceNodesId, faceIndex + i, faceVerticesMax, verticesFaceData.data() );
877 }
878 faceIndex += facesRead;
879 }
880 }
881
882 // Time values (not implemented)
883 mNcFile->putDataDouble( timeId, 0, 0.0 );
884
885 // Turning on define mode
886 nc_redef( mNcFile->handle() );
887 }
888
writeGlobals()889 void MDAL::DriverUgrid::writeGlobals()
890 {
891 mNcFile->putAttrStr( NC_GLOBAL, "source", "MDAL " + std::string( MDAL_Version() ) );
892 mNcFile->putAttrStr( NC_GLOBAL, "date_created", MDAL::getCurrentTimeStamp() );
893 mNcFile->putAttrStr( NC_GLOBAL, "Conventions", "CF-1.6 UGRID-1.0" );
894 }
895
parseClassification(int varid) const896 std::vector<std::pair<double, double>> MDAL::DriverUgrid::parseClassification( int varid ) const
897 {
898 std::vector<std::pair<double, double>> classes;
899 std::string flagBoundVarName = mNcFile->getAttrStr( "flag_bounds", varid );
900 if ( !flagBoundVarName.empty() )
901 {
902 try
903 {
904 int boundsVarId = mNcFile->getVarId( flagBoundVarName );
905 std::vector<size_t> classDims;
906 std::vector<int> classDimIds;
907 mNcFile->getDimensions( flagBoundVarName, classDims, classDimIds );
908 std::vector<double> boundValues = mNcFile->readDoubleArr( boundsVarId, 0, 0, classDims[0], classDims[1] );
909
910 if ( classDims[1] != 2 || classDims[0] <= 0 )
911 throw MDAL::Error( MDAL_Status::Err_UnknownFormat, "Invalid classification dimension" );
912
913 std::pair<std::string, std::string> classificationMeta;
914 classificationMeta.first = "classification";
915 std::string classification;
916 for ( size_t i = 0; i < classDims[0]; ++i )
917 {
918 std::pair<double, double> classBound;
919 classBound.first = boundValues[i * 2];
920 classBound.second = boundValues[i * 2 + 1];
921 classes.push_back( classBound );
922 }
923 }
924 catch ( MDAL::Error &err )
925 {
926 MDAL::Log::warning( err.status, err.driver, "Error when parsing class bounds: " + err.mssg + ", classification ignored" );
927 }
928 }
929
930 return classes;
931 }
932