diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt
index 3cc7122be..6c45b68cb 100644
--- a/libs/MeshKernel/CMakeLists.txt
+++ b/libs/MeshKernel/CMakeLists.txt
@@ -47,6 +47,7 @@ set(
   ${SRC_DIR}/Mesh2D.cpp
   ${SRC_DIR}/Mesh2DGenerateGlobal.cpp
   ${SRC_DIR}/Mesh2DIntersections.cpp
+  ${SRC_DIR}/Mesh2DToCurvilinear.cpp
   ${SRC_DIR}/MeshRefinement.cpp
   ${SRC_DIR}/Network1D.cpp
   ${SRC_DIR}/Operations.cpp
@@ -147,6 +148,7 @@ set(
   ${DOMAIN_INC_DIR}/Mesh2D.hpp
   ${DOMAIN_INC_DIR}/Mesh2DGenerateGlobal.hpp
   ${DOMAIN_INC_DIR}/Mesh2DIntersections.hpp
+  ${DOMAIN_INC_DIR}/Mesh2DToCurvilinear.hpp
   ${DOMAIN_INC_DIR}/MeshConversion.hpp
   ${DOMAIN_INC_DIR}/MeshInterpolation.hpp
   ${DOMAIN_INC_DIR}/MeshRefinement.hpp
diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp
new file mode 100644
index 000000000..776bab077
--- /dev/null
+++ b/libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp
@@ -0,0 +1,89 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2023.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation version 3.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// contact: delft3d.support@deltares.nl
+// Stichting Deltares
+// P.O. Box 177
+// 2600 MH Delft, The Netherlands
+//
+// All indications and logos of, and references to, "Delft3D" and "Deltares"
+// are registered trademarks of Stichting Deltares, and remain the property of
+// Stichting Deltares. All rights reserved.
+//
+//------------------------------------------------------------------------------
+
+#pragma once
+#include <memory>
+
+#include "MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp"
+#include "MeshKernel/Definitions.hpp"
+#include "MeshKernel/Mesh2D.hpp"
+#include "MeshKernel/Point.hpp"
+#include "Utilities/LinearAlgebra.hpp"
+
+using namespace meshkernel::constants;
+
+namespace meshkernel
+{
+    /// @brief Construct a curvilinear grid from an unstructured mesh
+    class Mesh2DToCurvilinear
+    {
+    public:
+        /// @brief Constructor
+        /// @param[in] mesh The input unstructured mesh
+        explicit Mesh2DToCurvilinear(Mesh2D& mesh);
+
+        /// @brief Computes the curvilinear grid starting from a specific point
+        /// @param[in] point The point from which start growing the curvilinear mesh. The point must be inside a quadrangular face
+        /// @returns The curvilinear grid
+        std::unique_ptr<CurvilinearGrid> Compute(const Point& point);
+
+    private:
+        /// @brief Computes the local mapping of the nodes composing the face
+        [[nodiscard]] Eigen::Matrix<UInt, 2, 2> ComputeLocalNodeMapping(UInt face) const;
+
+        /// @brief Computes the node indices of the neighbouring faces
+        [[nodiscard]] UInt ComputeNeighbouringFaceNodes(const UInt face,
+                                                        const Eigen::Matrix<UInt, 2, 2>& localNodeMapping,
+                                                        const UInt d,
+                                                        const std::vector<bool>& visitedFace);
+
+        /// @brief Computes the final curvilinear matrix
+        [[nodiscard]] lin_alg::Matrix<Point> ComputeCurvilinearMatrix();
+
+        Mesh2D& m_mesh;       ///< The mesh to convert
+        std::vector<int> m_i; ///< The i indices of each node on the curvilinear grid
+        std::vector<int> m_j; ///< The j indices of each node on the curvilinear grid
+
+        const std::array<std::array<int, 2>, 4> m_nodeFrom = {{{0, 0},
+                                                               {0, 0},
+                                                               {1, 0},
+                                                               {1, 1}}}; ///< starting edge node indices for each direction in the local mapping
+
+        const std::array<std::array<int, 2>, 4> m_nodeTo = {{{0, 1},
+                                                             {1, 0},
+                                                             {1, 1},
+                                                             {0, 1}}}; ///< ending edge node indices for each direction in the local mapping
+
+        const std::array<std::array<int, 2>, 4> m_directionsDeltas = {{{-1, 0},
+                                                                       {0, -1},
+                                                                       {1, 0},
+                                                                       {0, 1}}}; ///< increments for the new nodes depending on the node direction
+
+        const int n_maxNumRowsColumns = 1000000; ///< The maximum number of allowed rows or columns
+    };
+
+} // namespace meshkernel
diff --git a/libs/MeshKernel/src/Mesh2DToCurvilinear.cpp b/libs/MeshKernel/src/Mesh2DToCurvilinear.cpp
new file mode 100644
index 000000000..b71550568
--- /dev/null
+++ b/libs/MeshKernel/src/Mesh2DToCurvilinear.cpp
@@ -0,0 +1,296 @@
+//---- GPL ---------------------------------------------------------------------
+//
+// Copyright (C)  Stichting Deltares, 2011-2023.
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation version 3.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+// contact: delft3d.support@deltares.nl
+// Stichting Deltares
+// P.O. Box 177
+// 2600 MH Delft, The Netherlands
+//
+// All indications and logos of, and references to, "Delft3D" and "Deltares"
+// are registered trademarks of Stichting Deltares, and remain the property of
+// Stichting Deltares. All rights reserved.
+//
+//------------------------------------------------------------------------------
+
+#include <queue>
+
+#include "MeshKernel/Mesh2DToCurvilinear.hpp"
+
+using namespace meshkernel;
+using namespace constants;
+
+Mesh2DToCurvilinear::Mesh2DToCurvilinear(Mesh2D& mesh) : m_mesh(mesh)
+{
+    if (mesh.GetNumNodes() <= 0)
+    {
+        throw AlgorithmError("Mesh with no nodes");
+    }
+
+    mesh.Administrate();
+
+    if (mesh.GetNumFaces() <= 0)
+    {
+        throw AlgorithmError("Mesh with no faces");
+    }
+}
+
+std::unique_ptr<CurvilinearGrid> Mesh2DToCurvilinear::Compute(const Point& point)
+{
+    // 1. Find the face index
+    m_mesh.BuildTree(Location::Faces);
+    const auto initialFaceIndex = m_mesh.FindLocationIndex(point, Location::Faces);
+    if (m_mesh.GetNumFaceEdges(initialFaceIndex) != geometric::numNodesInQuadrilateral)
+    {
+        throw AlgorithmError("The initial face is not a quadrilateral");
+    }
+
+    // 2. Check if the point is inside the face
+    std::vector<Point> polygonPoints;
+    for (UInt n = 0; n < geometric::numNodesInQuadrilateral; ++n)
+    {
+        const auto node = m_mesh.m_facesNodes[initialFaceIndex][n];
+        polygonPoints.emplace_back(m_mesh.Node(node));
+    }
+    const auto node = m_mesh.m_facesNodes[initialFaceIndex][0];
+    polygonPoints.emplace_back(m_mesh.Node(node));
+
+    Polygon polygon(polygonPoints, m_mesh.m_projection);
+    if (!polygon.Contains(point))
+    {
+        throw AlgorithmError("The initial face does not contain the starting point");
+    }
+
+    // 3. Build the local coordinate system
+    const auto numNodes = m_mesh.GetNumNodes();
+    m_i = std::vector(numNodes, missing::intValue);
+    m_j = std::vector(numNodes, missing::intValue);
+
+    const auto firstEdge = m_mesh.m_facesEdges[initialFaceIndex][0];
+    const auto secondEdge = m_mesh.m_facesEdges[initialFaceIndex][1];
+    const auto thirdEdge = m_mesh.m_facesEdges[initialFaceIndex][2];
+    const auto fourthEdge = m_mesh.m_facesEdges[initialFaceIndex][3];
+
+    const auto firstNodeIndex = m_mesh.FindCommonNode(firstEdge, secondEdge);
+    m_i[firstNodeIndex] = 0;
+    m_j[firstNodeIndex] = 0;
+
+    const auto secondNodeIndex = m_mesh.FindCommonNode(secondEdge, thirdEdge);
+    m_i[secondNodeIndex] = 1;
+    m_j[secondNodeIndex] = 0;
+
+    const auto thirdNodeIndex = m_mesh.FindCommonNode(thirdEdge, fourthEdge);
+    m_i[thirdNodeIndex] = 1;
+    m_j[thirdNodeIndex] = 1;
+
+    const auto fourthNodeIndex = m_mesh.FindCommonNode(fourthEdge, firstEdge);
+    m_i[fourthNodeIndex] = 0;
+    m_j[fourthNodeIndex] = 1;
+
+    // 4. Grow the front using the breath first search algorithm
+    const auto numFaces = m_mesh.GetNumFaces();
+    std::vector visitedFace(numFaces, false);
+    std::queue<UInt> q;
+    q.push(initialFaceIndex);
+
+    while (!q.empty())
+    {
+        const auto face = q.front();
+        q.pop();
+
+        if (visitedFace[face])
+        {
+            continue;
+        }
+        visitedFace[face] = true;
+
+        if (m_mesh.GetNumFaceEdges(face) != geometric::numNodesInQuadrilateral)
+        {
+            continue;
+        }
+
+        const auto localNodeMapping = ComputeLocalNodeMapping(face);
+        for (UInt d = 0; d < geometric::numNodesInQuadrilateral; ++d)
+        {
+            const auto newFaceIndex = ComputeNeighbouringFaceNodes(face, localNodeMapping, d, visitedFace);
+            if (newFaceIndex != missing::uintValue)
+            {
+                q.push(newFaceIndex);
+            }
+        }
+    }
+    return std::make_unique<CurvilinearGrid>(ComputeCurvilinearMatrix(), m_mesh.m_projection);
+}
+
+Eigen::Matrix<UInt, 2, 2> Mesh2DToCurvilinear::ComputeLocalNodeMapping(UInt face) const
+{
+    const auto& faceIndices = m_mesh.m_facesNodes[face];
+
+    const auto node0 = faceIndices[0];
+    const auto node1 = faceIndices[1];
+    const auto node2 = faceIndices[2];
+    const auto node3 = faceIndices[3];
+
+    const std::vector localI{m_i[node0], m_i[node1], m_i[node2], m_i[node3]};
+    const std::vector localJ{m_j[node0], m_j[node1], m_j[node2], m_j[node3]};
+    const std::vector localNodes{node0, node1, node2, node3};
+
+    const auto minI = *std::ranges::min_element(localI);
+    const auto minJ = *std::ranges::min_element(localJ);
+
+    Eigen::Matrix<UInt, 2, 2> matrix;
+    for (UInt i = 0; i < localI.size(); ++i)
+    {
+        if (localI[i] == minI && localJ[i] == minJ)
+        {
+            matrix(0, 0) = localNodes[i];
+        }
+        if (localI[i] == minI + 1 && localJ[i] == minJ)
+        {
+            matrix(1, 0) = localNodes[i];
+        }
+        if (localI[i] == minI && localJ[i] == minJ + 1)
+        {
+            matrix(0, 1) = localNodes[i];
+        }
+        if (localI[i] == minI + 1 && localJ[i] == minJ + 1)
+        {
+            matrix(1, 1) = localNodes[i];
+        }
+    }
+    return matrix;
+}
+
+UInt Mesh2DToCurvilinear::ComputeNeighbouringFaceNodes(const UInt face,
+                                                       const Eigen::Matrix<UInt, 2, 2>& localNodeMapping,
+                                                       const UInt d,
+                                                       const std::vector<bool>& visitedFace)
+{
+
+    const auto firstNode = localNodeMapping(m_nodeFrom[d][0], m_nodeFrom[d][1]);
+    const auto secondNode = localNodeMapping(m_nodeTo[d][0], m_nodeTo[d][1]);
+
+    // find the edge index
+    const auto edgeIndex = m_mesh.FindEdge(firstNode, secondNode);
+    if (edgeIndex == missing::uintValue)
+    {
+        return missing::uintValue;
+    }
+
+    // this edge belongs only to the current face
+    if (m_mesh.m_edgesNumFaces[edgeIndex] < 2)
+    {
+        return missing::uintValue;
+    }
+
+    const auto newFace = face == m_mesh.m_edgesFaces[edgeIndex][0] ? m_mesh.m_edgesFaces[edgeIndex][1] : m_mesh.m_edgesFaces[edgeIndex][0];
+
+    if (visitedFace[newFace])
+    {
+        return missing::uintValue;
+    }
+
+    if (m_mesh.GetNumFaceEdges(newFace) != geometric::numNodesInQuadrilateral)
+    {
+        return missing::uintValue;
+    }
+
+    int edgeIndexInNewFace = 0;
+    for (UInt e = 0u; e < geometric::numNodesInQuadrilateral; ++e)
+    {
+        if (m_mesh.m_facesEdges[newFace][e] == edgeIndex)
+        {
+            edgeIndexInNewFace = static_cast<int>(e);
+            break;
+        }
+    }
+    auto nextEdgeIndexInNewFace = edgeIndexInNewFace + 1;
+    nextEdgeIndexInNewFace = nextEdgeIndexInNewFace == 4 ? 0 : nextEdgeIndexInNewFace;
+    const auto nextEdgeInNewFace = m_mesh.m_facesEdges[newFace][nextEdgeIndexInNewFace];
+    const auto firstCommonNode = m_mesh.FindCommonNode(edgeIndex, nextEdgeInNewFace);
+    const auto firstOtherNode = OtherNodeOfEdge(m_mesh.GetEdge(nextEdgeInNewFace), firstCommonNode);
+    const auto iFirstOtherNode = m_i[firstCommonNode] + m_directionsDeltas[d][0];
+    const auto jFirstOtherNode = m_j[firstCommonNode] + m_directionsDeltas[d][1];
+
+    auto previousEdgeIndexInNewFace = edgeIndexInNewFace - 1;
+    previousEdgeIndexInNewFace = previousEdgeIndexInNewFace == -1 ? 3 : previousEdgeIndexInNewFace;
+    const auto previousEdgeInNewFace = m_mesh.m_facesEdges[newFace][previousEdgeIndexInNewFace];
+    const auto secondCommonNode = m_mesh.FindCommonNode(edgeIndex, previousEdgeInNewFace);
+    const auto secondOtherNode = OtherNodeOfEdge(m_mesh.GetEdge(previousEdgeInNewFace), secondCommonNode);
+    const auto iSecondCommonNode = m_i[secondCommonNode] + m_directionsDeltas[d][0];
+    const auto jSecondCommonNode = m_j[secondCommonNode] + m_directionsDeltas[d][1];
+
+    const auto invalid = (m_i[firstOtherNode] != missing::intValue && m_i[firstOtherNode] != iFirstOtherNode) ||
+                         (m_j[firstOtherNode] != missing::intValue && m_j[firstOtherNode] != jFirstOtherNode) ||
+                         (m_i[secondOtherNode] != missing::intValue && m_i[secondOtherNode] != iSecondCommonNode) ||
+                         (m_j[secondOtherNode] != missing::intValue && m_j[secondOtherNode] != jSecondCommonNode);
+
+    if (!invalid)
+    {
+        m_i[firstOtherNode] = iFirstOtherNode;
+        m_j[firstOtherNode] = jFirstOtherNode;
+        m_i[secondOtherNode] = iSecondCommonNode;
+        m_j[secondOtherNode] = jSecondCommonNode;
+        return newFace;
+    }
+    return missing::uintValue;
+}
+
+lin_alg::Matrix<Point> Mesh2DToCurvilinear::ComputeCurvilinearMatrix()
+{
+
+    int minI = n_maxNumRowsColumns;
+    int minJ = n_maxNumRowsColumns;
+    const auto numNodes = m_mesh.GetNumNodes();
+    for (UInt n = 0; n < numNodes; ++n)
+    {
+        if (m_i[n] != missing::intValue)
+        {
+            minI = std::min(minI, m_i[n]);
+        }
+        if (m_j[n] != missing::intValue)
+        {
+            minJ = std::min(minJ, m_j[n]);
+        }
+    }
+
+    int maxI = -n_maxNumRowsColumns;
+    int maxJ = -n_maxNumRowsColumns;
+    for (UInt n = 0; n < numNodes; ++n)
+    {
+        if (m_i[n] != missing::intValue)
+        {
+            m_i[n] = m_i[n] - minI;
+            maxI = std::max(maxI, m_i[n]);
+        }
+        if (m_j[n] != missing::intValue)
+        {
+            m_j[n] = m_j[n] - minJ;
+            maxJ = std::max(maxJ, m_j[n]);
+        }
+    }
+
+    lin_alg::Matrix<Point> result(maxI + 1, maxJ + 1);
+    for (UInt n = 0; n < numNodes; ++n)
+    {
+        const auto i = m_i[n];
+        const auto j = m_j[n];
+        if (i != missing::intValue && j != missing::intValue)
+        {
+            result(i, j) = m_mesh.Node(n);
+        }
+    }
+    return result;
+}
diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp
index 6907512cc..f86bc0ed6 100644
--- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp
+++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp
@@ -1,20 +1,20 @@
-#include "MeshKernel/Mesh2DIntersections.hpp"
-
 #include <chrono>
 #include <gmock/gmock-matchers.h>
 #include <gtest/gtest.h>
 #include <random>
 
-#include <MeshKernel/Constants.hpp>
-#include <MeshKernel/Entities.hpp>
-#include <MeshKernel/Mesh.hpp>
-#include <MeshKernel/Mesh2D.hpp>
-#include <MeshKernel/Polygons.hpp>
-#include <MeshKernel/RemoveDisconnectedRegions.hpp>
-#include <TestUtils/Definitions.hpp>
-#include <TestUtils/MakeMeshes.hpp>
+#include "MeshKernel/Constants.hpp"
+#include "MeshKernel/Entities.hpp"
+#include "MeshKernel/Mesh.hpp"
+#include "MeshKernel/Mesh2D.hpp"
+#include "MeshKernel/Mesh2DIntersections.hpp"
+#include "MeshKernel/Mesh2DToCurvilinear.hpp"
+#include "MeshKernel/Operations.hpp"
+#include "MeshKernel/Polygons.hpp"
+#include "MeshKernel/RemoveDisconnectedRegions.hpp"
 
-#include <MeshKernel/Operations.hpp>
+#include "TestUtils/Definitions.hpp"
+#include "TestUtils/MakeMeshes.hpp"
 
 TEST(Mesh2D, OneQuadTestConstructor)
 {
@@ -1230,3 +1230,111 @@ TEST(Mesh2D, DeleteMesh_WithPolygonAndIncludedCircumcenters_ShouldDeleteInnerFac
         EXPECT_EQ(originalEdges[i].second, mesh->GetEdge(i).second);
     }
 }
+
+TEST(Mesh2D, Mesh2DToCurvilinear_WithRectangularMesh_ShouldCreateFullCurvilinearMesh)
+{
+    // Prepare
+    const auto mesh = MakeRectangularMeshForTesting(3,
+                                                    3,
+                                                    10.0,
+                                                    10.0,
+                                                    meshkernel::Projection::cartesian);
+    meshkernel::Mesh2DToCurvilinear mesh2DToCurvilinear(*mesh);
+
+    // Execute
+    const meshkernel::Point point(5.0, 5.0);
+    const auto curvilinearGrid = mesh2DToCurvilinear.Compute(point);
+
+    // Assert
+    ASSERT_EQ(3, curvilinearGrid->NumM());
+    ASSERT_EQ(3, curvilinearGrid->NumN());
+    ASSERT_EQ(9, curvilinearGrid->GetNumNodes());
+
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 2).x);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 2).y);
+
+    ASSERT_EQ(5.0, curvilinearGrid->GetNode(0, 1).x);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 1).y);
+
+    ASSERT_EQ(10.0, curvilinearGrid->GetNode(0, 0).x);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 0).y);
+}
+
+TEST(Mesh2D, Mesh2DToCurvilinear_WithStartingPointOutsideMesh_ShouldThrowAnException)
+{
+    // Prepare
+    const auto mesh = MakeRectangularMeshForTesting(3,
+                                                    3,
+                                                    10.0,
+                                                    10.0,
+                                                    meshkernel::Projection::cartesian);
+    meshkernel::Mesh2DToCurvilinear mesh2DToCurvilinear(*mesh);
+
+    // Execute
+    const meshkernel::Point point(-20.0, -20.0);
+
+    // Assert
+    EXPECT_THROW(mesh2DToCurvilinear.Compute(point), meshkernel::AlgorithmError);
+}
+
+TEST(Mesh2D, Mesh2DToCurvilinear_WithMixedMesh_ShouldCreatePartialCurvilinearMesh)
+{
+    // Prepare a mixed mesh with two triangles at the boundary
+    std::vector<meshkernel::Point> nodes;
+    nodes.push_back({-10.0, 0.0});
+    nodes.push_back({0.0, 0.0});
+    nodes.push_back({10.0, 0.0});
+    nodes.push_back({20.0, 0.0});
+    nodes.push_back({0.0, 10.0});
+    nodes.push_back({10.0, 10.0});
+    nodes.push_back({20.0, 10.0});
+    nodes.push_back({0.0, 20.0});
+    nodes.push_back({10.0, 20.0});
+
+    std::vector<meshkernel::Edge> edges;
+    edges.push_back({0, 4});
+    edges.push_back({0, 1});
+    edges.push_back({1, 4});
+    edges.push_back({1, 2});
+    edges.push_back({2, 5});
+    edges.push_back({2, 3});
+    edges.push_back({3, 6});
+    edges.push_back({4, 7});
+    edges.push_back({4, 5});
+    edges.push_back({5, 8});
+    edges.push_back({5, 6});
+    edges.push_back({7, 8});
+    edges.push_back({6, 8});
+
+    // 2 Execute
+    auto mesh = meshkernel::Mesh2D(edges, nodes, meshkernel::Projection::cartesian);
+    meshkernel::Mesh2DToCurvilinear mesh2DToCurvilinear(mesh);
+    const meshkernel::Point point(5.0, 5.0);
+    const auto curvilinearGrid = mesh2DToCurvilinear.Compute(point);
+
+    // Assert
+    ASSERT_EQ(3, curvilinearGrid->NumM());
+    ASSERT_EQ(3, curvilinearGrid->NumN());
+    ASSERT_EQ(9, curvilinearGrid->GetNumNodes());
+
+    ASSERT_EQ(20.0, curvilinearGrid->GetNode(0, 0).x);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 0).y);
+    ASSERT_EQ(20.0, curvilinearGrid->GetNode(1, 0).x);
+    ASSERT_EQ(10.0, curvilinearGrid->GetNode(1, 0).y);
+    ASSERT_EQ(-999.0, curvilinearGrid->GetNode(2, 0).x);
+    ASSERT_EQ(-999.0, curvilinearGrid->GetNode(2, 0).y);
+
+    ASSERT_EQ(10.0, curvilinearGrid->GetNode(0, 1).x);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 1).y);
+    ASSERT_EQ(10.0, curvilinearGrid->GetNode(1, 1).x);
+    ASSERT_EQ(10.0, curvilinearGrid->GetNode(1, 1).y);
+    ASSERT_EQ(10.0, curvilinearGrid->GetNode(2, 1).x);
+    ASSERT_EQ(20.0, curvilinearGrid->GetNode(2, 1).y);
+
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 2).x);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(0, 2).y);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(1, 2).x);
+    ASSERT_EQ(10.0, curvilinearGrid->GetNode(1, 2).y);
+    ASSERT_EQ(0.0, curvilinearGrid->GetNode(2, 2).x);
+    ASSERT_EQ(20.0, curvilinearGrid->GetNode(2, 2).y);
+}
diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp
index ddf518974..17c871892 100644
--- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp
+++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp
@@ -855,6 +855,14 @@ namespace meshkernelapi
         /// @returns Error code
         MKERNEL_API int mkernel_mesh2d_convert_projection(int meshKernelId, int projectionType, const char* const zoneString);
 
+        /// @brief Converts a mesh to a curvilinear mesh
+        ///
+        /// @param[in] meshKernelId The id of the mesh state
+        /// @param[in] xPointCoordinate The x coordinate of the point where to start the conversion point coordinate
+        /// @param[in] yPointCoordinate The y point coordinate
+        /// @returns Error code
+        MKERNEL_API int mkernel_mesh2d_convert_to_curvilinear(int meshKernelId, double xPointCoordinate, double yPointCoordinate);
+
         /// @brief Count the number of hanging edges in a mesh2d.
         /// An hanging edge is an edge where one of the two nodes is not connected.
         /// @param[in] meshKernelId The id of the mesh state
diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp
index 61d942c88..a6a2c1571 100644
--- a/libs/MeshKernelApi/src/MeshKernel.cpp
+++ b/libs/MeshKernelApi/src/MeshKernel.cpp
@@ -27,6 +27,7 @@
 
 #include "MeshKernel/CurvilinearGrid/CurvilinearGridDeleteExterior.hpp"
 #include "MeshKernel/Mesh2DIntersections.hpp"
+#include "MeshKernel/Mesh2DToCurvilinear.hpp"
 #include "MeshKernel/SamplesHessianCalculator.hpp"
 
 #include "MeshKernel/CurvilinearGrid/CurvilinearGridDeleteInterior.hpp"
@@ -770,6 +771,28 @@ namespace meshkernelapi
         return lastExitCode;
     }
 
+    MKERNEL_API int mkernel_mesh2d_convert_to_curvilinear(int meshKernelId, double xPointCoordinate, double yPointCoordinate)
+    {
+        lastExitCode = meshkernel::ExitCode::Success;
+        try
+        {
+            if (!meshKernelState.contains(meshKernelId))
+            {
+                throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist.");
+            }
+
+            meshkernel::Mesh2DToCurvilinear mesh2DToCurvilinear(*meshKernelState[meshKernelId].m_mesh2d);
+
+            meshKernelState[meshKernelId].m_curvilinearGrid = mesh2DToCurvilinear.Compute({xPointCoordinate, yPointCoordinate});
+            meshKernelState[meshKernelId].m_mesh2d = std::make_unique<meshkernel::Mesh2D>(meshKernelState[meshKernelId].m_projection);
+        }
+        catch (...)
+        {
+            lastExitCode = HandleException();
+        }
+        return lastExitCode;
+    }
+
     MKERNEL_API int mkernel_mesh2d_count_hanging_edges(int meshKernelId, int& numHangingEdges)
     {
         lastExitCode = meshkernel::ExitCode::Success;
diff --git a/libs/MeshKernelApi/tests/src/ApiTest.cpp b/libs/MeshKernelApi/tests/src/ApiTest.cpp
index edefd6241..74273d570 100644
--- a/libs/MeshKernelApi/tests/src/ApiTest.cpp
+++ b/libs/MeshKernelApi/tests/src/ApiTest.cpp
@@ -3248,3 +3248,44 @@ TEST_P(MeshLocationIndexTests, GetLocationIndex_OnAMesh_ShouldGetTheLocationInde
     ASSERT_EQ(locationIndex, expectedIndex);
 }
 INSTANTIATE_TEST_SUITE_P(LocationIndexParametrizedTests, MeshLocationIndexTests, ::testing::ValuesIn(MeshLocationIndexTests::GetData()));
+
+TEST(Mesh2D, ConvertToCurvilinear_ShouldConvertMeshToCurvilinear)
+{
+    // create first mesh
+    auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(10, 10, 1.0, meshkernel::Point(0.0, 0.0));
+
+    meshkernelapi::Mesh2D mesh2d{};
+    mesh2d.num_nodes = static_cast<int>(num_nodes);
+    mesh2d.num_edges = static_cast<int>(num_edges);
+    mesh2d.node_x = node_x.data();
+    mesh2d.node_y = node_y.data();
+    mesh2d.edge_nodes = edge_nodes.data();
+
+    // allocate state
+    int meshKernelId = 0;
+    int errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    // first initialise using the first mesh, mesh2d
+    errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    // Execute
+    errorCode = meshkernelapi::mkernel_mesh2d_convert_to_curvilinear(meshKernelId, 5.0, 5.0);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    meshkernelapi::Mesh2D mesh2dOut{};
+    errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2dOut);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    meshkernelapi::CurvilinearGrid curvilinearOut{};
+    errorCode = mkernel_curvilinear_get_dimensions(meshKernelId, curvilinearOut);
+    ASSERT_EQ(meshkernel::ExitCode::Success, errorCode);
+
+    // Assert
+    ASSERT_EQ(0, mesh2dOut.num_nodes);
+    ASSERT_EQ(0, mesh2dOut.num_edges);
+
+    ASSERT_EQ(11, curvilinearOut.num_m);
+    ASSERT_EQ(11, curvilinearOut.num_n);
+}