Skip to content

Commit

Permalink
update MDAL to 0.0.7 (3di support, projections support, bugfixes)
Browse files Browse the repository at this point in the history
  • Loading branch information
@PeterPetrik @nyalldawson
PeterPetrik authored and nyalldawson committed Sep 10, 2018
1 parent d17912c commit 97c9580
Show file tree
Hide file tree
Showing 19 changed files with 1,147 additions and 20 deletions.
5 changes: 5 additions & 0 deletions external/mdal/api/mdal.h
View file
Expand Up @@ -83,6 +83,11 @@ MDAL_EXPORT MeshH MDAL_LoadMesh( const char *meshFile );
//! Closes mesh, frees the memory
MDAL_EXPORT void MDAL_CloseMesh( MeshH mesh );

/**
* Returns mesh projection
* not thread-safe and valid only till next call
*/
MDAL_EXPORT const char *MDAL_M_projection( MeshH mesh );
//! Returns vertex count for the mesh
MDAL_EXPORT int MDAL_M_vertexCount( MeshH mesh );
//! Returns vertex X coord for the mesh
Expand Down
31 changes: 31 additions & 0 deletions external/mdal/cmake/FindNetCDF.cmake
View file
@@ -0,0 +1,31 @@
# - Find NetCDF
# Find the native NetCDF includes and library
# Copyright (c) 2018, Peter Petrik <zilolv at gmail dot com>
#
# This module returns these variables for the rest of the project to use.
#
# NETCDF_FOUND - True if NetCDF found including required interfaces (see below)
# NETCDF_LIBRARY - All netcdf related libraries
# NETCDF_INCLUDE_DIR - All directories to include

IF (NETCDF_INCLUDE_DIR AND NETCDF_LIBRARY)
# Already in cache, be silent
SET (NETCDF_FIND_QUIETLY TRUE)
ENDIF ()

FIND_PACKAGE(PkgConfig)
PKG_CHECK_MODULES(PC_NETCDF QUIET netcdf)
SET(NETCDF_DEFINITIONS ${PC_NETCDF_CFLAGS_OTHER})

FIND_PATH (NETCDF_INCLUDE_DIR netcdf.h
HINTS ${PC_NETCDF_INCLUDEDIR} ${PC_NETCDF_INCLUDE_DIRS} ${NETCDF_PREFIX}/include
PATH_SUFFIXES libnetcdf )

FIND_LIBRARY (NETCDF_LIBRARY
NAMES netcdf libnetcdf
HINTS HINTS ${PC_NETCDF_LIBDIR} ${PC_NETCDF_LIBRARY_DIRS} ${NETCDF_PREFIX}/lib)

INCLUDE (FindPackageHandleStandardArgs)
FIND_PACKAGE_HANDLE_STANDARD_ARGS (NetCDF
DEFAULT_MSG NETCDF_LIBRARY NETCDF_INCLUDE_DIR)

2 changes: 2 additions & 0 deletions external/mdal/cmake_templates/mdal_config.hpp.in
View file
Expand Up @@ -7,6 +7,8 @@

#cmakedefine HAVE_GDAL

#cmakedefine HAVE_NETCDF

#endif // MDAL_CONFIG_HPP


233 changes: 233 additions & 0 deletions external/mdal/frmts/mdal_3di.cpp
View file
@@ -0,0 +1,233 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2018 Peter Petrik (zilolv at gmail dot com)
*/

#include "mdal_3di.hpp"

MDAL::Loader3Di::Loader3Di( const std::string &fileName )
: LoaderCF( fileName )
{
}

MDAL::CFDimensions MDAL::Loader3Di::populateDimensions()
{
CFDimensions dims;
size_t count;
int ncid;

// 2D Mesh
mNcFile.getDimension( "nMesh2D_nodes", &count, &ncid );
dims.setDimension( CFDimensions::Face2D, count, ncid );

mNcFile.getDimension( "nCorner_Nodes", &count, &ncid );
dims.setDimension( CFDimensions::MaxVerticesInFace, count, ncid );

// Vertices count is populated later in populateFacesAndVertices
// it is not known from the array dimensions

// Time
mNcFile.getDimension( "time", &count, &ncid );
dims.setDimension( CFDimensions::Time, count, ncid );

return dims;
}

void MDAL::Loader3Di::populateFacesAndVertices( MDAL::Mesh *mesh )
{
assert( mesh );
size_t faceCount = mDimensions.size( CFDimensions::Face2D );
mesh->faces.resize( faceCount );
size_t verticesInFace = mDimensions.size( CFDimensions::MaxVerticesInFace );
size_t arrsize = faceCount * verticesInFace;
std::map<std::string, size_t> xyToVertex2DId;

// X coordinate
int ncidX = mNcFile.getVarId( "Mesh2DContour_x" );
double fillX = mNcFile.getFillValue( ncidX );
std::vector<double> faceVerticesX( arrsize );
if ( nc_get_var_double( mNcFile.handle(), ncidX, faceVerticesX.data() ) ) throw MDAL_Status::Err_UnknownFormat;

// Y coordinate
int ncidY = mNcFile.getVarId( "Mesh2DContour_y" );
double fillY = mNcFile.getFillValue( ncidY );
std::vector<double> faceVerticesY( arrsize );
if ( nc_get_var_double( mNcFile.handle(), ncidY, faceVerticesY.data() ) ) throw MDAL_Status::Err_UnknownFormat;

// now populate create faces and backtrack which vertices
// are used in multiple faces
for ( size_t faceId = 0; faceId < faceCount; ++faceId )
{
Face face;

for ( size_t faceVertexId = 0; faceVertexId < verticesInFace; ++faceVertexId )
{
size_t arrId = faceId * verticesInFace + faceVertexId;
Vertex vertex;
vertex.x = faceVerticesX[arrId];
vertex.y = faceVerticesY[arrId];
vertex.z = 0;

if ( MDAL::equals( vertex.x, fillX ) || MDAL::equals( vertex.y, fillY ) )
break;


size_t vertexId;

std::string key = std::to_string( vertex.x ) + "," + std::to_string( vertex.y );
const auto it = xyToVertex2DId.find( key );
if ( it == xyToVertex2DId.end() )
{
// new vertex
vertexId = mesh->vertices.size();
xyToVertex2DId[key] = vertexId;
mesh->vertices.push_back( vertex );
}
else
{
// existing vertex
vertexId = it->second;
}

face.push_back( vertexId );

}

mesh->faces[faceId] = face;
}

// Only now we have number of vertices, since we identified vertices that
// are used in multiple faces
mDimensions.setDimension( CFDimensions::Vertex2D, mesh->vertices.size() );
}

void MDAL::Loader3Di::addBedElevation( MDAL::Mesh *mesh )
{
assert( mesh );
if ( mesh->faces.empty() )
return;

size_t faceCount = mesh->faces.size();

// read Z coordinate of 3di computation nodes centers
int ncidZ = mNcFile.getVarId( "Mesh2DFace_zcc" );
double fillZ = mNcFile.getFillValue( ncidZ );
std::vector<double> coordZ( faceCount );
if ( nc_get_var_double( mNcFile.handle(), ncidZ, coordZ.data() ) )
return; //error reading the array


std::shared_ptr<DatasetGroup> group = std::make_shared< DatasetGroup >();
group->isOnVertices = false;
group->isScalar = true;
group->setName( "Bed Elevation" );
group->uri = mesh->uri;
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< Dataset >();
dataset->time = 0.0;
dataset->values.resize( faceCount );
dataset->active.resize( faceCount );
dataset->parent = group.get();
for ( size_t i = 0; i < faceCount; ++i )
{
dataset->values[i].x = MDAL::safeValue( coordZ[i], fillZ );
}
group->datasets.push_back( dataset );
mesh->datasetGroups.push_back( group );
}

std::string MDAL::Loader3Di::getCoordinateSystemVariableName()
{
return "projected_coordinate_system";
}

std::set<std::string> MDAL::Loader3Di::ignoreNetCDFVariables()
{
std::set<std::string> ignore_variables;

ignore_variables.insert( "projected_coordinate_system" );
ignore_variables.insert( "time" );

std::vector<std::string> meshes;
meshes.push_back( "Mesh1D" );
meshes.push_back( "Mesh2D" );

for ( const std::string &mesh : meshes )
{
ignore_variables.insert( mesh );
ignore_variables.insert( mesh + "Node_id" );
ignore_variables.insert( mesh + "Node_type" );

ignore_variables.insert( mesh + "Face_xcc" );
ignore_variables.insert( mesh + "Face_ycc" );
ignore_variables.insert( mesh + "Face_zcc" );
ignore_variables.insert( mesh + "Contour_x" );
ignore_variables.insert( mesh + "Contour_y" );
ignore_variables.insert( mesh + "Face_sumax" );

ignore_variables.insert( mesh + "Line_id" );
ignore_variables.insert( mesh + "Line_xcc" );
ignore_variables.insert( mesh + "Line_ycc" );
ignore_variables.insert( mesh + "Line_zcc" );
ignore_variables.insert( mesh + "Line_type" );
}

return ignore_variables;
}

std::string MDAL::Loader3Di::nameSuffix( MDAL::CFDimensions::Type type )
{
MDAL_UNUSED( type );
return "";
}

void MDAL::Loader3Di::parseNetCDFVariableMetadata( int varid, const std::string &variableName, std::string &name, bool *is_vector, bool *is_x )
{
*is_vector = false;
*is_x = true;

std::string long_name = mNcFile.getAttrStr( "long_name", varid );
if ( long_name.empty() )
{
std::string standard_name = mNcFile.getAttrStr( "standard_name", varid );
if ( standard_name.empty() )
{
name = variableName;
}
else
{
if ( MDAL::contains( standard_name, "_x_" ) )
{
*is_vector = true;
name = MDAL::replace( standard_name, "_x_", "" );
}
else if ( MDAL::contains( standard_name, "_y_" ) )
{
*is_vector = true;
*is_x = false;
name = MDAL::replace( standard_name, "_y_", "" );
}
else
{
name = standard_name;
}
}
}
else
{
if ( MDAL::contains( long_name, " in x direction" ) )
{
*is_vector = true;
name = MDAL::replace( long_name, " in x direction", "" );
}
else if ( MDAL::contains( long_name, " in y direction" ) )
{
*is_vector = true;
*is_x = false;
name = MDAL::replace( long_name, " in y direction", "" );
}
else
{
name = long_name;
}
}
}
62 changes: 62 additions & 0 deletions external/mdal/frmts/mdal_3di.hpp
View file
@@ -0,0 +1,62 @@
/*
MDAL - Mesh Data Abstraction Library (MIT License)
Copyright (C) 2018 Peter Petrik (zilolv at gmail dot com)
*/

#ifndef MDAL_3DI_HPP
#define MDAL_3DI_HPP

#include <map>
#include <string>
#include <stddef.h>

#include "mdal_cf.hpp"

namespace MDAL
{

/**
* Loader of 3Di file format.
*
* The result 3Di NetCDF file is based on CF-conventions with some additions.
* It is unstructured grid with data stored in NetCDF/HDF5 file format.
* A division is made between a 1D and 2D which can be distinguished through
* the prefixes “MESH2D” and “MESH1D”. For both meshes the information is present
* in coordinate, id, and type variables.
*
* A version of the data scheme is not present yet.
*
* The 2D Mesh consists of calculation Nodes that represents centers of Faces.
* There is no concept of Vertices in the file. The vertices that forms a face
* are specified by X,Y coordinates in the "Face Contours" arrays. The "lines"
* represent the face's edges and are again specified by X,Y coordinate of the
* line center. Data is specified on calculation nodes (i.e. dataset defined on faces)
* and on lines (i.e. dataset defined on edges - not implemented yet)
*
* The 1D Mesh is present too, but not parsed yet.
*/
class Loader3Di: public LoaderCF
{
public:
Loader3Di( const std::string &fileName );
~Loader3Di() override = default;

private:
CFDimensions populateDimensions() override;
void populateFacesAndVertices( MDAL::Mesh *mesh ) override;
void addBedElevation( MDAL::Mesh *mesh ) override;
std::string getCoordinateSystemVariableName() override;
std::set<std::string> ignoreNetCDFVariables() override;
std::string nameSuffix( CFDimensions::Type type ) override;
void parseNetCDFVariableMetadata( int varid, const std::string &variableName,
std::string &name, bool *is_vector, bool *is_x ) override;

//! Returns number of vertices
size_t parse2DMesh();

void addBedElevationDatasetOnFaces();
};

} // namespace MDAL

#endif // MDAL_3DI_HPP
8 changes: 4 additions & 4 deletions external/mdal/frmts/mdal_ascii_dat.cpp
View file
Expand Up @@ -69,7 +69,7 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
oldFormat = true;
isVector = ( line == "VECTOR" );

group.reset( new DatasetGroup() );
group = std::make_shared< DatasetGroup >();
group->uri = mDatFile;
group->setName( name );
group->isScalar = !isVector;
Expand Down Expand Up @@ -118,7 +118,7 @@ void MDAL::LoaderAsciiDat::load( MDAL::Mesh *mesh, MDAL_Status *status )
}
isVector = cardType == "BEGVEC";

group.reset( new DatasetGroup() );
group = std::make_shared< DatasetGroup >();
group->uri = mDatFile;
group->setName( name );
group->isScalar = !isVector;
Expand Down Expand Up @@ -197,7 +197,7 @@ void MDAL::LoaderAsciiDat::readVertexTimestep(
size_t vertexCount = mesh->vertices.size();
size_t faceCount = mesh->faces.size();

std::shared_ptr<MDAL::Dataset> dataset( new MDAL::Dataset );
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< MDAL::Dataset >();
dataset->time = t / 3600.; // TODO read TIMEUNITS
dataset->values.resize( vertexCount );
dataset->active.resize( faceCount );
Expand Down Expand Up @@ -267,7 +267,7 @@ void MDAL::LoaderAsciiDat::readFaceTimestep(

size_t faceCount = mesh->faces.size();

std::shared_ptr<MDAL::Dataset> dataset( new MDAL::Dataset );
std::shared_ptr<MDAL::Dataset> dataset = std::make_shared< MDAL::Dataset >();
dataset->time = t / 3600.;
dataset->values.resize( faceCount );
dataset->parent = group.get();
Expand Down

0 comments on commit 97c9580

Please sign in to comment.