diff --git a/libs/MeshKernel/include/MeshKernel/Definitions.hpp b/libs/MeshKernel/include/MeshKernel/Definitions.hpp
index d6cade34fe..a521e4ee65 100644
--- a/libs/MeshKernel/include/MeshKernel/Definitions.hpp
+++ b/libs/MeshKernel/include/MeshKernel/Definitions.hpp
@@ -113,70 +113,4 @@ namespace meshkernel
     /// @brief Get the string representation of the CurvilinearDirection enumeration values.
     const std::string& CurvilinearDirectionToString(CurvilinearDirection direction);
 
-    /// @brief The concept specifies that the array type must have an access operator returning the array element type or can be converted to one
-    template <typename ArrayType, typename ResultType>
-    concept ArrayConstAccessConcept = requires(const ArrayType& array, const size_t i) {
-                                          {
-                                              array[i]
-                                              } -> std::convertible_to<ResultType>;
-                                      };
-
-    /// @brief The concept specifies that the array type must have an access operator returning a reference to the array element type
-    template <typename ArrayType, typename ResultType>
-    concept ArrayNonConstAccessConcept = requires(ArrayType& array, const size_t i) {
-                                             {
-                                                 array[i]
-                                                 } -> std::same_as<ResultType&>;
-                                         };
-
-    /// @brief A concept that specifies that the array must have a size function return the number of elements in the array
-    template <typename ArrayType>
-    concept ArraySizeConcept = requires(const ArrayType& array) {
-                                   {
-                                       array.size()
-                                       } -> std::same_as<size_t>;
-                               };
-
-    /// @brief A concept that specifies that the array must have a begin and end function.
-    ///
-    /// Would like to also specify the return type here, but span needs some c++23 functionality here.
-    /// Then change all iterator usage to cbegin and cend returning a const_iterator
-    /// std::same_as<typename ArrayType::const_iterator>
-    template <typename ArrayType>
-    concept ArrayConstIteratorsConcept = requires(const ArrayType& array) {
-                                             {
-                                                 array.begin()
-                                             };
-                                             {
-                                                 array.end()
-                                             };
-                                         };
-
-    /// @brief A concept that specifies that the array must have a begin and end function.
-    ///
-    /// Would like to also specify the return type here, but span needs some c++23 functionality here.
-    /// Then change all iterator usage to cbegin and cend returning a const_iterator
-    /// std::same_as<typename ArrayType::const_iterator>
-    template <typename ArrayType>
-    concept ArrayNonConstIteratorsConcept = requires(ArrayType& array) {
-                                                {
-                                                    array.begin()
-                                                    } -> std::same_as<typename ArrayType::iterator>;
-                                                {
-                                                    array.end()
-                                                    } -> std::same_as<typename ArrayType::iterator>;
-                                            };
-
-    /// @brief A concept that specifies all the functionality required to be usable as a constant array of doubles.
-    template <typename ArrayType>
-    concept ValidConstDoubleArray = ArrayConstAccessConcept<ArrayType, double> &&
-                                    ArrayConstIteratorsConcept<ArrayType> &&
-                                    ArraySizeConcept<ArrayType>;
-
-    /// @brief A concept that specifies all the functionality required to be usable as a constant array of doubles.
-    template <typename ArrayType>
-    concept ValidNonConstDoubleArray = ArrayNonConstAccessConcept<ArrayType, double> &&
-                                       ArrayNonConstIteratorsConcept<ArrayType> &&
-                                       ArraySizeConcept<ArrayType>;
-
 } // namespace meshkernel
diff --git a/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp b/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp
index a39260197a..b83cbddf56 100644
--- a/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp
+++ b/libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp
@@ -30,6 +30,7 @@
 #include <array>
 #include <functional>
 #include <memory>
+#include <span>
 #include <string>
 #include <vector>
 
@@ -49,9 +50,6 @@ namespace meshkernel
     class BoundedArray
     {
     public:
-        /// @brief Constructor
-        BoundedArray() : m_size(0) {}
-
         /// @brief Number of elements in the array
         UInt size() const
         {
@@ -66,7 +64,6 @@ namespace meshkernel
                 throw ConstraintError("array already at size limit: {}", Dimension);
             }
 
-            // m_indices.push_back(index);
             m_indices[m_size] = index;
             ++m_size;
         }
@@ -117,14 +114,12 @@ namespace meshkernel
     {
     public:
         /// @brief Constructor with array of points
-        template <meshkernel::ValidConstPointArray PointVector>
-        MeshTriangulation(const PointVector& nodes,
+        MeshTriangulation(const std::span<const Point> nodes,
                           const Projection projection);
 
         /// @brief Constructor with separate arrays of x- and y-coordinates
-        template <meshkernel::ValidConstDoubleArray VectorType>
-        MeshTriangulation(const VectorType& xNodes,
-                          const VectorType& yNodes,
+        MeshTriangulation(const std::span<const double> xNodes,
+                          const std::span<const double> yNodes,
                           const Projection projection);
 
         /// @brief Get the projection type used in the triangulation
@@ -199,62 +194,6 @@ inline meshkernel::Projection meshkernel::MeshTriangulation::GetProjection() con
     return m_projection;
 }
 
-template <meshkernel::ValidConstPointArray PointVector>
-meshkernel::MeshTriangulation::MeshTriangulation(const PointVector& nodes,
-                                                 const Projection projection)
-    : m_nodes(nodes.begin(), nodes.end()),
-      m_projection(projection)
-{
-    if (nodes.size() < constants::geometric::numNodesInTriangle)
-    {
-        throw ConstraintError("Not enough nodes for a single triangle: {}", nodes.size());
-    }
-
-    std::vector<double> xNodes(nodes.size());
-    std::vector<double> yNodes(nodes.size());
-
-    std::transform(nodes.begin(), nodes.end(), xNodes.begin(),
-                   [](const Point& p)
-                   { return p.x; });
-    std::transform(nodes.begin(), nodes.end(), yNodes.begin(),
-                   [](const Point& p)
-                   { return p.y; });
-
-    std::span<double> xNodesSpan(xNodes.data(), xNodes.size());
-    std::span<double> yNodesSpan(yNodes.data(), yNodes.size());
-
-    Compute(xNodesSpan, yNodesSpan);
-}
-
-template <meshkernel::ValidConstDoubleArray VectorType>
-meshkernel::MeshTriangulation::MeshTriangulation(const VectorType& xNodes,
-                                                 const VectorType& yNodes,
-                                                 const Projection projection)
-    : m_nodes(xNodes.size()),
-      m_projection(projection)
-{
-    if (xNodes.size() < constants::geometric::numNodesInTriangle)
-    {
-        throw ConstraintError("Not enough nodes for a single triangle: {}", xNodes.size());
-    }
-
-    if (xNodes.size() != yNodes.size())
-    {
-        throw ConstraintError("Size mismatch between x- and y-node values: {} /= {}",
-                              xNodes.size(), yNodes.size());
-    }
-
-    for (size_t i = 0; i < xNodes.size(); ++i)
-    {
-        m_nodes[i] = Point{xNodes[i], yNodes[i]};
-    }
-
-    std::span<const double> xNodesSpan(xNodes.data(), xNodes.size());
-    std::span<const double> yNodesSpan(yNodes.data(), yNodes.size());
-
-    Compute(xNodesSpan, yNodesSpan);
-}
-
 inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfNodes() const
 {
     return static_cast<UInt>(m_nodes.size());
diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp
index ccd3ea4509..a34bf5296b 100644
--- a/libs/MeshKernel/include/MeshKernel/Operations.hpp
+++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp
@@ -332,9 +332,8 @@ namespace meshkernel
     /// @param[in] triangleNodes A series of node forming an open triangle
     /// @param[in] projection The coordinate system projection.
     /// @returns If point is inside the designated triangle
-    template <class PointVector>
     [[nodiscard]] bool IsPointInTriangle(const Point& point,
-                                         const PointVector& triangleNodes,
+                                         const std::span<const Point> triangleNodes,
                                          const Projection& projection);
 
     /// @brief Computes three base components
@@ -595,8 +594,7 @@ namespace meshkernel
     /// @param[in] points The point series.
     /// @param[in] projection The projection to use.
     /// @return The average coordinate.
-    template <ValidConstPointArray PointVector>
-    [[nodiscard]] Point ComputeAverageCoordinate(const PointVector& points, const Projection& projection);
+    [[nodiscard]] Point ComputeAverageCoordinate(const std::span<const Point> points, const Projection& projection);
 
     /// @brief Cartesian projection of a point on a segment defined by other two points
     /// @param firstNode The first node of the segment
@@ -675,120 +673,3 @@ inline meshkernel::Cartesian3DPoint meshkernel::operator*(const Cartesian3DPoint
 {
     return value * p;
 }
-
-namespace meshkernel
-{
-    template <ValidConstPointArray PointVector>
-    [[nodiscard]] Point ComputeAverageCoordinate(const PointVector& points, const Projection& projection)
-    {
-        size_t validCount = std::ranges::count_if(points, [](const Point& p)
-                                                  { return p.IsValid(); });
-
-        if (projection == Projection::sphericalAccurate)
-        {
-
-            UInt firstValidPoint = 0;
-
-            if (validCount != points.size())
-            {
-                auto iterator = std::ranges::find_if(points, [](const Point& p)
-                                                     { return p.IsValid(); });
-                firstValidPoint = static_cast<UInt>(iterator - points.begin());
-            }
-
-            Cartesian3DPoint averagePoint3D{0.0, 0.0, 0.0};
-            for (const auto& point : points)
-            {
-                if (!point.IsValid())
-                {
-                    continue;
-                }
-
-                const Cartesian3DPoint point3D{SphericalToCartesian3D(point)};
-                averagePoint3D.x += point3D.x;
-                averagePoint3D.y += point3D.y;
-                averagePoint3D.z += point3D.z;
-            }
-            averagePoint3D.x = averagePoint3D.x / static_cast<double>(validCount);
-            averagePoint3D.y = averagePoint3D.y / static_cast<double>(validCount);
-            averagePoint3D.z = averagePoint3D.z / static_cast<double>(validCount);
-
-            return Cartesian3DToSpherical(averagePoint3D, points[firstValidPoint].x);
-        }
-
-        auto result = std::accumulate(points.begin(), points.end(), Point{0.0, 0.0}, [](const Point& sum, const Point& current)
-                                      { return current.IsValid() ? sum + current : sum; });
-        result.x = result.x / static_cast<double>(validCount);
-        result.y = result.y / static_cast<double>(validCount);
-        return result;
-    }
-} // namespace meshkernel
-
-template <class PointVector>
-[[nodiscard]] bool meshkernel::IsPointInTriangle(const Point& point,
-                                                 const PointVector& triangleNodes,
-                                                 const Projection& projection)
-{
-    if (triangleNodes.empty())
-    {
-        return true;
-    }
-
-    bool isInTriangle = false;
-
-    if (projection == Projection::cartesian || projection == Projection::spherical)
-    {
-        int windingNumber = 0;
-
-        for (UInt n = 0; n < constants::geometric::numNodesInTriangle; n++)
-        {
-            UInt endIndex = n == 2 ? 0 : n + 1;
-
-            const auto crossProductValue = crossProduct(triangleNodes[n], triangleNodes[endIndex], triangleNodes[n], point, Projection::cartesian);
-
-            if (IsEqual(crossProductValue, 0.0))
-            {
-                // point on the line
-                return true;
-            }
-
-            if (triangleNodes[n].y <= point.y) // an upward crossing
-            {
-                if (triangleNodes[endIndex].y > point.y && crossProductValue > 0.0)
-                {
-                    ++windingNumber; // have  a valid up intersect
-                }
-            }
-            else
-            {
-                if (triangleNodes[endIndex].y <= point.y && crossProductValue < 0.0) // a downward crossing
-                {
-                    --windingNumber; // have  a valid down intersect
-                }
-            }
-        }
-
-        isInTriangle = windingNumber == 0 ? false : true;
-    }
-
-    if (projection == Projection::sphericalAccurate)
-    {
-        std::vector<Point> triangleCopy;
-        triangleCopy.reserve(constants::geometric::numNodesInTriangle + 1);
-        Point centre(0.0, 0.0);
-
-        for (UInt i = 0; i < constants::geometric::numNodesInTriangle; ++i)
-        {
-            centre += triangleNodes[i];
-            triangleCopy.emplace_back(triangleNodes[i]);
-        }
-
-        triangleCopy.emplace_back(triangleNodes[0]);
-
-        centre *= 1.0 / 3.0;
-
-        isInTriangle = IsPointInPolygonNodes(point, triangleCopy, projection, centre);
-    }
-
-    return isInTriangle;
-}
diff --git a/libs/MeshKernel/include/MeshKernel/Point.hpp b/libs/MeshKernel/include/MeshKernel/Point.hpp
index 89ccc864ae..2fd4859081 100644
--- a/libs/MeshKernel/include/MeshKernel/Point.hpp
+++ b/libs/MeshKernel/include/MeshKernel/Point.hpp
@@ -265,12 +265,6 @@ namespace meshkernel
     /// @brief Get the delta-x and -y in Cartesian coordinate system
     Vector GetDeltaCartesian(const Point& p1, const Point& p2);
 
-    /// @brief A concept that specifies all the functionality required to be usable as an array of points.
-    template <typename ArrayType>
-    concept ValidConstPointArray = ArrayConstAccessConcept<ArrayType, Point> &&
-                                   ArrayConstIteratorsConcept<ArrayType> &&
-                                   ArraySizeConcept<ArrayType>;
-
 } // namespace meshkernel
 
 inline void meshkernel::Point::SetInvalid()
diff --git a/libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp b/libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp
index 458edb67a6..c97b14e162 100644
--- a/libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp
+++ b/libs/MeshKernel/include/MeshKernel/SampleAveragingInterpolator.hpp
@@ -70,45 +70,28 @@ namespace meshkernel
     {
     public:
         /// @brief Constructor.
-        ///
-        /// The VectorType can be any array type of double precision values, e.g. std::vector, std::span.
-        template <meshkernel::ValidConstDoubleArray VectorType>
-        SampleAveragingInterpolator(const VectorType& xNodes,
-                                    const VectorType& yNodes,
+        SampleAveragingInterpolator(const std::span<const double> xNodes,
+                                    const std::span<const double> yNodes,
                                     const Projection projection,
-                                    const InterpolationParameters& interpolationParameters)
-            : m_samplePoints(CombineCoordinates(xNodes, yNodes)),
-              m_projection(projection),
-              m_interpolationParameters(interpolationParameters),
-              m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method,
-                                                                                   interpolationParameters.m_minimumNumberOfSamples,
-                                                                                   projection)),
-              m_nodeRTree(RTreeFactory::Create(projection))
-        {
-            m_nodeRTree->BuildTree(m_samplePoints);
-        }
+                                    const InterpolationParameters& interpolationParameters);
 
         /// @brief Constructor.
-        ///
-        /// The VectorType can be any array type of double precision values, e.g. std::vector, std::span.
-        template <meshkernel::ValidConstPointArray PointVector>
-        SampleAveragingInterpolator(const PointVector& nodes,
+        SampleAveragingInterpolator(const std::span<const Point> nodes,
                                     const Projection projection,
-                                    const InterpolationParameters& interpolationParameters)
-            : m_samplePoints(nodes.begin(), nodes.end()),
-              m_projection(projection),
-              m_interpolationParameters(interpolationParameters),
-              m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method,
-                                                                                   interpolationParameters.m_minimumNumberOfSamples,
-                                                                                   projection)),
-              m_nodeRTree(RTreeFactory::Create(projection))
-        {
-            m_nodeRTree->BuildTree(m_samplePoints);
-        }
+                                    const InterpolationParameters& interpolationParameters);
 
         /// @brief Get the number of nodes of size of the sample data.
         UInt Size() const override;
 
+        /// @brief Set sample data from std::span object
+        void SetData(const int propertyId, const std::span<const double> sampleData) override;
+
+        /// @brief Interpolate the sample data set at the interpolation nodes.
+        void Interpolate(const int propertyId, const std::span<const Point> iterpolationNodes, std::span<double> result) const override;
+
+        /// @brief Interpolate the sample data set at the locationd defined.
+        void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double> result) const override;
+
         /// @brief Interpolate the sample data set at a single interpolation point.
         ///
         /// If interpolation at multiple points is required then better performance
@@ -122,28 +105,7 @@ namespace meshkernel
         static constexpr UInt MaximumNumberOfEdgesPerNode = 16; ///< Maximum number of edges per node
 
         /// @brief Combine x- and y-coordinate arrays to a point array
-        template <meshkernel::ValidConstDoubleArray VectorType>
-        static std::vector<Point> CombineCoordinates(const VectorType& xNodes, const VectorType& yNodes)
-        {
-            std::vector<Point> result(xNodes.size());
-
-            for (size_t i = 0; i < xNodes.size(); ++i)
-            {
-                result[i].x = xNodes[i];
-                result[i].y = yNodes[i];
-            }
-
-            return result;
-        }
-
-        /// @brief Set sample data from std::span object
-        void SetDataSpan(const int propertyId, const std::span<const double>& sampleData) override;
-
-        /// @brief Interpolate the sample data set at the interpolation nodes.
-        void InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const override;
-
-        /// @brief Interpolate the sample data set at the locationd defined.
-        void InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const override;
+        static std::vector<Point> CombineCoordinates(const std::span<const double> xNodes, const std::span<const double> yNodes);
 
         /// @brief Interpolate the sample data on the element at the interpolation point.
         double InterpolateOnElement(const UInt elementId, const Point& interpolationPoint, const std::vector<double>& sampleValues) const;
diff --git a/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp b/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp
index 4ba1348dbe..7c6660a6a1 100644
--- a/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp
+++ b/libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp
@@ -32,7 +32,7 @@
 #include <string>
 #include <vector>
 
-#include "MeshKernel/AveragingInterpolation.hpp" //  Only for the enum Method. should move to Definitions.hpp
+#include "MeshKernel/AveragingInterpolation.hpp"
 #include "MeshKernel/AveragingStrategies/AveragingStrategy.hpp"
 #include "MeshKernel/AveragingStrategies/AveragingStrategyFactory.hpp"
 #include "MeshKernel/Constants.hpp"
@@ -58,29 +58,14 @@ namespace meshkernel
         virtual UInt Size() const = 0;
 
         /// @brief Set sample data
-        template <meshkernel::ValidConstDoubleArray VectorType>
-        void SetData(const int propertyId, const VectorType& sampleData)
-        {
-            const std::span<const double> spanSampleData(sampleData.data(), sampleData.size());
-            SetDataSpan(propertyId, sampleData);
-        }
+        virtual void SetData(const int propertyId, const std::span<const double> sampleData) = 0;
 
         /// @brief Interpolate the sample data set at the interpolation nodes.
-        template <meshkernel::ValidConstPointArray PointVectorType, meshkernel::ValidConstDoubleArray ScalarVectorType>
-        void Interpolate(const int propertyId, const PointVectorType& iterpolationNodes, ScalarVectorType& result) const
-        {
-            const std::span<const Point> spanInterpolationNodes(iterpolationNodes.data(), iterpolationNodes.size());
-            std::span<double> spanResult(result.data(), result.size());
-            InterpolateSpan(propertyId, spanInterpolationNodes, spanResult);
-        }
+        virtual void Interpolate(const int propertyId, const std::span<const Point> iterpolationNodes, std::span<double> result) const = 0;
 
-        /// @brief Interpolate the sample data set at the interpolation point from a mesh.
-        template <meshkernel::ValidConstDoubleArray ScalarVectorType>
-        void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, ScalarVectorType& result) const
-        {
-            std::span<double> spanResult(result.data(), result.size());
-            InterpolateSpan(propertyId, mesh, location, spanResult);
-        }
+        /// @brief Interpolate the sample data.
+        // TODO may need a polygon too: const Polygons& polygon
+        virtual void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double> result) const = 0;
 
         /// @brief Interpolate the sample data set at a single interpolation point.
         ///
@@ -92,16 +77,6 @@ namespace meshkernel
         virtual bool Contains(const int propertyId) const = 0;
 
     private:
-        /// @brief Set sample data contained in an std::span object
-        virtual void SetDataSpan(const int propertyId, const std::span<const double>& sampleData) = 0;
-
-        /// @brief Interpolate the sample data set at the interpolation nodes.
-        virtual void InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const = 0;
-
-        /// @brief Interpolate the sample data.
-        // TODO may need a polygon too: const Polygons& polygon
-        virtual void InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const = 0;
-
         // TODO the map -> sample-data can be here
         // with protected members to access the data.
     };
diff --git a/libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp b/libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp
index ec4150471c..7243c32d47 100644
--- a/libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp
+++ b/libs/MeshKernel/include/MeshKernel/SampleTriangulationInterpolator.hpp
@@ -51,28 +51,27 @@ namespace meshkernel
     {
     public:
         /// @brief Constructor.
-        ///
-        /// The VectorType can be any array type of double precision values, e.g. std::vector, std::span.
-        template <meshkernel::ValidConstDoubleArray VectorType>
-        SampleTriangulationInterpolator(const VectorType& xNodes,
-                                        const VectorType& yNodes,
+        SampleTriangulationInterpolator(const std::span<const double> xNodes,
+                                        const std::span<const double> yNodes,
                                         const Projection projection)
             : m_triangulation(xNodes, yNodes, projection) {}
 
         /// @brief Constructor.
-        ///
-        /// The VectorType can be any array type of double precision values, e.g. std::vector, std::span.
-        template <meshkernel::ValidConstPointArray PointVector>
-        SampleTriangulationInterpolator(const PointVector& nodes,
+        SampleTriangulationInterpolator(const std::span<const Point> nodes,
                                         const Projection projection)
             : m_triangulation(nodes, projection) {}
 
         /// @brief Get the number of nodes of size of the sample data.
         UInt Size() const override;
 
-        /// @brief Set sample data
-        // template <meshkernel::ValidConstDoubleArray VectorType>
-        // void SetData(const int propertyId, const VectorType& sampleData) override;
+        /// @brief Set sample data from std::span object
+        void SetData(const int propertyId, const std::span<const double> sampleData) override;
+
+        /// @brief Interpolate the sample data at the points for the location (nodes, edges, faces)
+        void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double> result) const override;
+
+        /// @brief Interpolate the sample data set at the interpolation nodes.
+        void Interpolate(const int propertyId, const std::span<const Point> iterpolationNodes, std::span<double> result) const override;
 
         /// @brief Interpolate the sample data set at a single interpolation point.
         ///
@@ -84,15 +83,6 @@ namespace meshkernel
         bool Contains(const int propertyId) const override;
 
     private:
-        /// @brief Set sample data from std::span object
-        void SetDataSpan(const int propertyId, const std::span<const double>& sampleData) override;
-
-        /// @brief Interpolate the sample data at the points for the location (nodes, edges, faces)
-        void InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const override;
-
-        /// @brief Interpolate the sample data set at the interpolation nodes.
-        void InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const override;
-
         /// @brief Interpolate the sample data on the element at the interpolation point.
         double InterpolateOnElement(const UInt elementId, const Point& interpolationPoint, const std::vector<double>& sampleValues) const;
 
diff --git a/libs/MeshKernel/src/MeshTriangulation.cpp b/libs/MeshKernel/src/MeshTriangulation.cpp
index 362ae54bc0..c1ae8faac3 100644
--- a/libs/MeshKernel/src/MeshTriangulation.cpp
+++ b/libs/MeshKernel/src/MeshTriangulation.cpp
@@ -52,6 +52,54 @@ extern "C"
                        double trisize);
 }
 
+meshkernel::MeshTriangulation::MeshTriangulation(const std::span<const Point> nodes,
+                                                 const Projection projection)
+    : m_nodes(nodes.begin(), nodes.end()),
+      m_projection(projection)
+{
+    if (nodes.size() < constants::geometric::numNodesInTriangle)
+    {
+        throw ConstraintError("Not enough nodes for a single triangle: {}", nodes.size());
+    }
+
+    std::vector<double> xNodes(nodes.size());
+    std::vector<double> yNodes(nodes.size());
+
+    std::transform(nodes.begin(), nodes.end(), xNodes.begin(),
+                   [](const Point& p)
+                   { return p.x; });
+    std::transform(nodes.begin(), nodes.end(), yNodes.begin(),
+                   [](const Point& p)
+                   { return p.y; });
+
+    Compute(xNodes, yNodes);
+}
+
+meshkernel::MeshTriangulation::MeshTriangulation(const std::span<const double> xNodes,
+                                                 const std::span<const double> yNodes,
+                                                 const Projection projection)
+    : m_nodes(xNodes.size()),
+      m_projection(projection)
+{
+    if (xNodes.size() < constants::geometric::numNodesInTriangle)
+    {
+        throw ConstraintError("Not enough nodes for a single triangle: {}", xNodes.size());
+    }
+
+    if (xNodes.size() != yNodes.size())
+    {
+        throw ConstraintError("Size mismatch between x- and y-node values: {} /= {}",
+                              xNodes.size(), yNodes.size());
+    }
+
+    for (size_t i = 0; i < xNodes.size(); ++i)
+    {
+        m_nodes[i] = Point{xNodes[i], yNodes[i]};
+    }
+
+    Compute(xNodes, yNodes);
+}
+
 void meshkernel::MeshTriangulation::Compute(const std::span<const double>& xNodes,
                                             const std::span<const double>& yNodes)
 {
diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp
index cb820f6cad..98d538333c 100644
--- a/libs/MeshKernel/src/Operations.cpp
+++ b/libs/MeshKernel/src/Operations.cpp
@@ -286,6 +286,74 @@ namespace meshkernel
         return dx1 * dy2 - dy1 * dx2;
     }
 
+    bool IsPointInTriangle(const Point& point,
+                           const std::span<const Point> triangleNodes,
+                           const Projection& projection)
+    {
+        if (triangleNodes.empty())
+        {
+            return true;
+        }
+
+        bool isInTriangle = false;
+
+        if (projection == Projection::cartesian || projection == Projection::spherical)
+        {
+            int windingNumber = 0;
+
+            for (UInt n = 0; n < constants::geometric::numNodesInTriangle; n++)
+            {
+                UInt endIndex = n == 2 ? 0 : n + 1;
+
+                const auto crossProductValue = crossProduct(triangleNodes[n], triangleNodes[endIndex], triangleNodes[n], point, Projection::cartesian);
+
+                if (IsEqual(crossProductValue, 0.0))
+                {
+                    // point on the line
+                    return true;
+                }
+
+                if (triangleNodes[n].y <= point.y) // an upward crossing
+                {
+                    if (triangleNodes[endIndex].y > point.y && crossProductValue > 0.0)
+                    {
+                        ++windingNumber; // have  a valid up intersect
+                    }
+                }
+                else
+                {
+                    if (triangleNodes[endIndex].y <= point.y && crossProductValue < 0.0) // a downward crossing
+                    {
+                        --windingNumber; // have  a valid down intersect
+                    }
+                }
+            }
+
+            isInTriangle = windingNumber == 0 ? false : true;
+        }
+
+        if (projection == Projection::sphericalAccurate)
+        {
+            std::vector<Point> triangleCopy;
+            triangleCopy.reserve(constants::geometric::numNodesInTriangle + 1);
+            Point centre(0.0, 0.0);
+
+            for (UInt i = 0; i < constants::geometric::numNodesInTriangle; ++i)
+            {
+                centre += triangleNodes[i];
+                triangleCopy.emplace_back(triangleNodes[i]);
+            }
+
+            triangleCopy.emplace_back(triangleNodes[0]);
+
+            centre *= 1.0 / 3.0;
+
+            isInTriangle = IsPointInPolygonNodes(point, triangleCopy, projection, centre);
+        }
+
+        return isInTriangle;
+    }
+
     bool IsPointInPolygonNodes(const Point& point,
                                const std::vector<Point>& polygonNodes,
                                const Projection& projection,
@@ -1579,6 +1647,50 @@ namespace meshkernel
         return result;
     }
 
+    meshkernel::Point ComputeAverageCoordinate(const std::span<const Point> points, const Projection& projection)
+    {
+        size_t validCount = std::ranges::count_if(points, [](const Point& p)
+                                                  { return p.IsValid(); });
+
+        if (projection == Projection::sphericalAccurate)
+        {
+
+            UInt firstValidPoint = 0;
+
+            if (validCount != points.size())
+            {
+                auto iterator = std::ranges::find_if(points, [](const Point& p)
+                                                     { return p.IsValid(); });
+                firstValidPoint = static_cast<UInt>(iterator - points.begin());
+            }
+
+            Cartesian3DPoint averagePoint3D{0.0, 0.0, 0.0};
+            for (const auto& point : points)
+            {
+                if (!point.IsValid())
+                {
+                    continue;
+                }
+
+                const Cartesian3DPoint point3D{SphericalToCartesian3D(point)};
+                averagePoint3D.x += point3D.x;
+                averagePoint3D.y += point3D.y;
+                averagePoint3D.z += point3D.z;
+            }
+            averagePoint3D.x = averagePoint3D.x / static_cast<double>(validCount);
+            averagePoint3D.y = averagePoint3D.y / static_cast<double>(validCount);
+            averagePoint3D.z = averagePoint3D.z / static_cast<double>(validCount);
+
+            return Cartesian3DToSpherical(averagePoint3D, points[firstValidPoint].x);
+        }
+
+        auto result = std::accumulate(points.begin(), points.end(), Point{0.0, 0.0}, [](const Point& sum, const Point& current)
+                                      { return current.IsValid() ? sum + current : sum; });
+        result.x = result.x / static_cast<double>(validCount);
+        result.y = result.y / static_cast<double>(validCount);
+        return result;
+    }
+
     std::tuple<Point, double, bool> OrthogonalProjectionOnSegment(Point const& firstNode,
                                                                   Point const& secondNode,
                                                                   Point const& pointToProject)
diff --git a/libs/MeshKernel/src/SampleAveragingInterpolator.cpp b/libs/MeshKernel/src/SampleAveragingInterpolator.cpp
index 22a9a26d84..5885cc8798 100644
--- a/libs/MeshKernel/src/SampleAveragingInterpolator.cpp
+++ b/libs/MeshKernel/src/SampleAveragingInterpolator.cpp
@@ -28,7 +28,50 @@
 #include "MeshKernel/SampleAveragingInterpolator.hpp"
 #include "MeshKernel/Operations.hpp"
 
-void meshkernel::SampleAveragingInterpolator::SetDataSpan(const int propertyId, const std::span<const double>& sampleData)
+std::vector<meshkernel::Point> meshkernel::SampleAveragingInterpolator::CombineCoordinates(const std::span<const double> xNodes,
+                                                                                           const std::span<const double> yNodes)
+{
+    std::vector<Point> result(xNodes.size());
+
+    for (size_t i = 0; i < xNodes.size(); ++i)
+    {
+        result[i].x = xNodes[i];
+        result[i].y = yNodes[i];
+    }
+
+    return result;
+}
+
+meshkernel::SampleAveragingInterpolator::SampleAveragingInterpolator(const std::span<const double> xNodes,
+                                                                     const std::span<const double> yNodes,
+                                                                     const Projection projection,
+                                                                     const InterpolationParameters& interpolationParameters)
+    : m_samplePoints(CombineCoordinates(xNodes, yNodes)),
+      m_projection(projection),
+      m_interpolationParameters(interpolationParameters),
+      m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method,
+                                                                           interpolationParameters.m_minimumNumberOfSamples,
+                                                                           projection)),
+      m_nodeRTree(RTreeFactory::Create(projection))
+{
+    m_nodeRTree->BuildTree(m_samplePoints);
+}
+
+meshkernel::SampleAveragingInterpolator::SampleAveragingInterpolator(const std::span<const Point> nodes,
+                                                                     const Projection projection,
+                                                                     const InterpolationParameters& interpolationParameters)
+    : m_samplePoints(nodes.begin(), nodes.end()),
+      m_projection(projection),
+      m_interpolationParameters(interpolationParameters),
+      m_strategy(averaging::AveragingStrategyFactory::GetAveragingStrategy(interpolationParameters.m_method,
+                                                                           interpolationParameters.m_minimumNumberOfSamples,
+                                                                           projection)),
+      m_nodeRTree(RTreeFactory::Create(projection))
+{
+    m_nodeRTree->BuildTree(m_samplePoints);
+}
+
+void meshkernel::SampleAveragingInterpolator::SetData(const int propertyId, const std::span<const double> sampleData)
 {
     if (m_samplePoints.size() != sampleData.size())
     {
@@ -39,7 +82,7 @@ void meshkernel::SampleAveragingInterpolator::SetDataSpan(const int propertyId,
     m_sampleData[propertyId].assign(sampleData.begin(), sampleData.end());
 }
 
-void meshkernel::SampleAveragingInterpolator::InterpolateSpan(const int propertyId, const std::span<const Point>& interpolationNodes, std::span<double>& result) const
+void meshkernel::SampleAveragingInterpolator::Interpolate(const int propertyId, const std::span<const Point> interpolationNodes, std::span<double> result) const
 {
     const std::vector<double>& propertyValues = m_sampleData.at(propertyId);
     std::vector<Sample> sampleCache;
@@ -293,7 +336,7 @@ void meshkernel::SampleAveragingInterpolator::InterpolateAtFaces(const int prope
     }
 }
 
-void meshkernel::SampleAveragingInterpolator::InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const
+void meshkernel::SampleAveragingInterpolator::Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double> result) const
 {
     std::ranges::fill(result, constants::missing::doubleValue);
 
diff --git a/libs/MeshKernel/src/SampleTriangulationInterpolator.cpp b/libs/MeshKernel/src/SampleTriangulationInterpolator.cpp
index 042d754476..6a2ce481d6 100644
--- a/libs/MeshKernel/src/SampleTriangulationInterpolator.cpp
+++ b/libs/MeshKernel/src/SampleTriangulationInterpolator.cpp
@@ -1,6 +1,6 @@
 #include "MeshKernel/SampleTriangulationInterpolator.hpp"
 
-void meshkernel::SampleTriangulationInterpolator::SetDataSpan(const int propertyId, const std::span<const double>& sampleData)
+void meshkernel::SampleTriangulationInterpolator::SetData(const int propertyId, const std::span<const double> sampleData)
 {
     if (m_triangulation.NumberOfNodes() != sampleData.size())
     {
@@ -11,7 +11,7 @@ void meshkernel::SampleTriangulationInterpolator::SetDataSpan(const int property
     m_sampleData[propertyId].assign(sampleData.begin(), sampleData.end());
 }
 
-void meshkernel::SampleTriangulationInterpolator::InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const
+void meshkernel::SampleTriangulationInterpolator::Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double> result) const
 {
     if (!Contains(propertyId))
     {
@@ -35,10 +35,10 @@ void meshkernel::SampleTriangulationInterpolator::InterpolateSpan(const int prop
         throw ConstraintError("Unknown location");
     }
 
-    InterpolateSpan(propertyId, meshNodes, result);
+    Interpolate(propertyId, meshNodes, result);
 }
 
-void meshkernel::SampleTriangulationInterpolator::InterpolateSpan(const int propertyId, const std::span<const Point>& interpolationNodes, std::span<double>& result) const
+void meshkernel::SampleTriangulationInterpolator::Interpolate(const int propertyId, const std::span<const Point> interpolationNodes, std::span<double> result) const
 {
     if (!Contains(propertyId))
     {
diff --git a/libs/MeshKernel/tests/src/MeshPropertyTests.cpp b/libs/MeshKernel/tests/src/MeshPropertyTests.cpp
index db5439039f..6ebb10e094 100644
--- a/libs/MeshKernel/tests/src/MeshPropertyTests.cpp
+++ b/libs/MeshKernel/tests/src/MeshPropertyTests.cpp
@@ -80,9 +80,12 @@ TEST(MeshPropertyTests, TriangulationTest)
 TEST(MeshPropertyTests, TriangulationFailureTest)
 {
     // Not enough points
-    EXPECT_THROW(mk::MeshTriangulation triangulation(std::vector<double>{0.0, 1.0}, std::vector<double>{0.0, 1.0}, mk::Projection::cartesian), mk::ConstraintError);
+    std::vector<double> twoPointsX{0.0, 1.0};
+    std::vector<double> twoPointsY{0.0, 1.0};
+    EXPECT_THROW(mk::MeshTriangulation triangulation(twoPointsX, twoPointsY, mk::Projection::cartesian), mk::ConstraintError);
     // x- and y-points not same size
-    EXPECT_THROW(mk::MeshTriangulation triangulation(std::vector<double>{0.0, 1.0, 2.0}, std::vector<double>{0.0, 1.0}, mk::Projection::cartesian), mk::ConstraintError);
+    std::vector<double> threePointsX{0.0, 1.0, 2.0};
+    EXPECT_THROW(mk::MeshTriangulation triangulation(threePointsX, twoPointsY, mk::Projection::cartesian), mk::ConstraintError);
 
     std::vector<double> xValues{0.0, 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 3.0, 0.0, 1.0, 2.0, 3.0};
     std::vector<double> yValues{0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0};