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