GRIDEDIT-1558 Fixed bug when snapping polygon with hole to land boundary

libs/MeshKernel/tests/src/PolygonsTests.cpp

@@ -1114,3 +1114,56 @@ TEST(Polygons, LinearRefinePolygonOutOfBoundsCheck)
     // polygon end node index incorrect
     EXPECT_THROW(auto refinedPolygon = polygons.LinearRefinePolygon(0, 1, 9), meshkernel::ConstraintError);
 }
+
+TEST(Polygons, SnapPolygonWithHoleToLandBoundary)
+{
+    // Test the algorithm for snapping multi-polygons to land boundaries.
+
+    //constexpr double tolerance = 1.0e-8;
+
+    // The land boundary to which the polygon is to be snapped.
+    // This land boundary is made up of two sections
+    std::vector<meshkernel::Point> landBoundaryPoints{{-1.0, -0.2},
+                                                      {11.0, -0.2}};
+
+    meshkernel::LandBoundary landBoundary(landBoundaryPoints);
+
+    // The original polygon points.
+    // The points make up two distinct polygons
+    std::vector<meshkernel::Point> polygonPoints{{0.0, 0.0},
+                                                 {10.0, 0.0},
+                                                 {10.0, 10.0},
+                                                 {0.0, 10.0},
+                                                 {0.0, 0.0},
+                                                 {meshkernel::constants::missing::innerOuterSeparator, meshkernel::constants::missing::innerOuterSeparator},
+                                                 {2.0, 2.0},
+                                                 {7.0, 2.0},
+                                                 {7.0, 7.0},
+                                                 {2.0, 7.0},
+                                                 {2.0, 2.0}};
+
+    meshkernel::Polygons polygon(polygonPoints, meshkernel::Projection::cartesian);
+
+    // Snap the polygon to the land boundary
+    // Only snaps the first polygon
+    polygon.SnapToLandBoundary(landBoundary, 0, 0);
+
+    const meshkernel::Polygon& firstPolygon = polygon.Enclosure(0).Outer();
+    const meshkernel::Polygon& innerPolygon = polygon.Enclosure(0).Inner(0);
+
+    std::cout << "first poly ";
+
+    for (size_t i = 0; i < firstPolygon.Size (); ++i) {
+        std::cout << "{" << firstPolygon.Node (i).x << "  " << firstPolygon.Node (i).y << "} ";
+    }
+
+    std::cout << std::endl;
+
+    std::cout << "inner poly ";
+
+    for (size_t i = 0; i < innerPolygon.Size (); ++i) {
+        std::cout << "{" << innerPolygon.Node (i).x << "  " << innerPolygon.Node (i).y << "} ";
+    }
+
+    std::cout << std::endl;
+}

libs/MeshKernelApi/src/MeshKernel.cpp

@@ -3373,19 +3373,36 @@ namespace meshkernelapi
             meshkernel::LandBoundary landBoundary(landBoundaryPoints);
             meshkernel::Polygons polygons(polygonPoints, meshKernelState[meshKernelId].m_mesh2d->m_projection);
 
-            polygons.SnapToLandBoundary(landBoundary, startIndex, endIndex);
             const auto [enclosureIndex, polygonStartNode, polygonEndNode] = polygons.PolygonIndex(startIndex, endIndex);
+            polygons.SnapToLandBoundary(landBoundary, startIndex, endIndex);
 
-            //--------------------------------
-            // Now copy back the polygon values
+            meshkernel::UInt enclosureStartPoint = 0;
+            meshkernel::UInt enclosurePreviousStartPoint = 0;
+            meshkernel::UInt enclosureCount = 0;
+
+            // Find the index of the first point in the polygnal enclosure with index enclosureIndex
+            for (meshkernel::UInt i = 0; i < polygonPoints.size (); ++i)
+            {
+                if ((polygonPoints [i].x == meshkernel::constants::missing::doubleValue && polygonPoints [i].y == meshkernel::constants::missing::doubleValue) || i == polygonPoints.size () - 1)
+                {
+                    if (enclosureCount == enclosureIndex) {
+                        enclosureStartPoint = enclosurePreviousStartPoint;
+                        break;
+                    }
+
+                    enclosurePreviousStartPoint = i + 1;
+                    ++enclosureCount;
+                }
+            }
 
             const std::vector<meshkernel::Point>& snappedPolygonPoints = polygons.Enclosure(enclosureIndex).Outer().Nodes();
 
-            for (int i = 0; i < polygon.num_coordinates; ++i)
+            for (size_t i = 0; i < snappedPolygonPoints.size (); ++i)
             {
-                polygon.coordinates_x[i] = snappedPolygonPoints[i].x;
-                polygon.coordinates_y[i] = snappedPolygonPoints[i].y;
+                polygon.coordinates_x[i + enclosureStartPoint] = snappedPolygonPoints[i].x;
+                polygon.coordinates_y[i + enclosureStartPoint] = snappedPolygonPoints[i].y;
             }
+
         }
         catch (...)
         {

libs/MeshKernelApi/tests/src/PolygonTests.cpp

@@ -38,3 +38,57 @@ TEST(PolygonTests, PolygonRefinementTests)
     ASSERT_NEAR(10.0, yPolygon[3], tolerance);
     ASSERT_NEAR(10.0, yPolygon[4], tolerance);
 }
+
+TEST(PolygonTests, PolygonWithHoleSnappingTest)
+{
+    constexpr double tolerance = 1.0e-8;
+
+    int meshKernelId;
+    auto errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    meshkernelapi::GeometryList land;
+    land.geometry_separator = meshkernel::constants::missing::doubleValue;
+    std::vector xLand{-1.0, 31.0};
+    std::vector yLand{-0.25, -0.25};
+    land.coordinates_x = xLand.data();
+    land.coordinates_y = yLand.data();
+    land.num_coordinates = static_cast<int>(xLand.size());
+
+    meshkernelapi::GeometryList polygon;
+    polygon.geometry_separator = meshkernel::constants::missing::doubleValue;
+    polygon.inner_outer_separator = meshkernel::constants::missing::innerOuterSeparator;
+
+    std::vector xPolygon{0.0, 10.0, 10.0, 0.0, 0.0, meshkernel::constants::missing::innerOuterSeparator, 2.0, 7.0, 7.0, 2.0, 2.0, meshkernel::constants::missing::doubleValue, 20.0, 30.0, 30.0, 20.0, 20.0};
+    std::vector yPolygon{0.0, 0.0, 10.0, 10.0, 0.0, meshkernel::constants::missing::innerOuterSeparator, 2.0, 2.0, 7.0, 7.0, 2.0, meshkernel::constants::missing::doubleValue, 0.0, 0.0, 10.0, 10.0, 0.0};
+
+    polygon.coordinates_x = xPolygon.data();
+    polygon.coordinates_y = yPolygon.data();
+    polygon.num_coordinates = static_cast<int>(xPolygon.size());
+
+    errorCode = mkernel_polygon_snap_to_landboundary(meshKernelId, land, polygon, 4, 1);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    const std::vector<double> expectedXAfterFirstSnapping{0.0, 10.0, 10.0, 0.0, 0.0, -998.0, 2.0, 7.0, 7.0, 2.0, 2.0, -999.0, 20.0, 30.0, 30.0, 20.0, 20.0};
+    const std::vector<double> expectedYAfterFirstSnapping{-0.25, -0.25, 10.0, 10.0, -0.25, -998.0, 2.0, 2.0, 7.0, 7.0, 2.0, -999.0, 0.0, 0.0, 10.0, 10.0, 0.0};
+
+    ASSERT_EQ(polygon.num_coordinates, static_cast<int>(expectedXAfterFirstSnapping.size()));
+
+    for (size_t i = 0; i < expectedXAfterFirstSnapping.size(); ++i)
+    {
+        EXPECT_NEAR(polygon.coordinates_x[i], expectedXAfterFirstSnapping[i], tolerance);
+        EXPECT_NEAR(polygon.coordinates_y[i], expectedYAfterFirstSnapping[i], tolerance);
+    }
+
+    errorCode = mkernel_polygon_snap_to_landboundary(meshKernelId, land, polygon, 16, 13);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    const std::vector<double> expectedXAfterSecondSnapping{0.0, 10.0, 10.0, 0.0, 0.0, -998.0, 2.0, 7.0, 7.0, 2.0, 2.0, -999.0, 20.0, 30.0, 30.0, 20.0, 20.0};
+    const std::vector<double> expectedYAfterSecondSnapping{-0.25, -0.25, 10.0, 10.0, -0.25, -998.0, 2.0, 2.0, 7.0, 7.0, 2.0, -999.0, -0.25, -0.25, 10.0, 10.0, -0.25};
+
+    for (size_t i = 0; i < expectedXAfterSecondSnapping.size(); ++i)
+    {
+        EXPECT_NEAR(polygon.coordinates_x[i], expectedXAfterSecondSnapping[i], tolerance);
+        EXPECT_NEAR(polygon.coordinates_y[i], expectedYAfterSecondSnapping[i], tolerance);
+    }
+}