From 1e7bd1a0983a1cd2b88da8cf848aa58648e6c805 Mon Sep 17 00:00:00 2001
From: rinkk <karsten.rink@ufz.de>
Date: Wed, 17 Jun 2020 11:14:44 +0200
Subject: [PATCH] [meshlib] added flag that allows ignoring nodata values
 during mesh mapping

---
 MeshLib/MeshGenerators/MeshLayerMapper.cpp | 25 ++++++++++++++++------
 MeshLib/MeshGenerators/MeshLayerMapper.h   |  5 ++++-
 2 files changed, 22 insertions(+), 8 deletions(-)

diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index c2ae402cfb3..b333efd6723 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -269,8 +269,10 @@ void MeshLayerMapper::addLayerToMesh(const MeshLib::Mesh &dem_mesh, unsigned lay
 
 bool MeshLayerMapper::layerMapping(MeshLib::Mesh const& mesh,
                                    GeoLib::Raster const& raster,
-                                   double const noDataReplacementValue = 0.0)
+                                   double const nodata_replacement,
+                                   bool const ignore_nodata)
 {
+
     if (mesh.getDimension() != 2)
     {
         ERR("MshLayerMapper::layerMapping() - requires 2D mesh");
@@ -289,18 +291,27 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh const& mesh,
     const std::vector<MeshLib::Node*> &nodes = mesh.getNodes();
     for (unsigned i = 0; i < nNodes; ++i)
     {
-        if (!raster.isPntOnRaster(*nodes[i]))
+        if (!ignore_nodata && !raster.isPntOnRaster(*nodes[i]))
         {
             // use either default value or elevation from layer above
-            nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1], noDataReplacementValue);
+            nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1],
+                                        nodata_replacement);
             continue;
         }
 
-        double elevation (raster.interpolateValueAtPoint(*nodes[i]));
-        if (std::abs(elevation - header.no_data) <
-            std::numeric_limits<double>::epsilon())
+        double elevation (raster.getValueAtPoint(*nodes[i]));
+        constexpr double eps = std::numeric_limits<double>::epsilon();
+        if (std::fabs(elevation - header.no_data) < eps)
+        {
+            if (ignore_nodata)
+            {
+                continue;
+            }
+            elevation = nodata_replacement;
+        }
+        else
         {
-            elevation = noDataReplacementValue;
+            elevation = raster.interpolateValueAtPoint(*nodes[i]);
         }
         nodes[i]->updateCoordinates((*nodes[i])[0], (*nodes[i])[1], elevation);
     }
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.h b/MeshLib/MeshGenerators/MeshLayerMapper.h
index c9f8f2411bd..f7a8429609e 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.h
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.h
@@ -59,7 +59,10 @@ public:
     * locations where no information is given, node elevation is set to
     * noDataReplacementValue.
     */
-    static bool layerMapping(MeshLib::Mesh const &mesh, const GeoLib::Raster &raster, double noDataReplacementValue);
+    static bool layerMapping(MeshLib::Mesh const& mesh,
+                             const GeoLib::Raster& raster,
+                             double nodata_replacement = 0.0,
+                             bool const ignore_nodata = false);
 
     /// Maps the elevation of all mesh nodes to the specified static value.
     static bool mapToStaticValue(MeshLib::Mesh const& mesh, double value);
-- 
GitLab