diff --git a/src/operation/BoundaryOp.cpp b/src/operation/BoundaryOp.cpp
index b407a2a41c..c0c452c455 100644
--- a/src/operation/BoundaryOp.cpp
+++ b/src/operation/BoundaryOp.cpp
@@ -139,7 +139,9 @@ BoundaryOp::boundaryMultiLineString(const geom::MultiLineString& mLine)
 
     // return Point or MultiPoint
     if (bdyPts->size() == 1) {
-        return std::unique_ptr<Geometry>(m_geomFact.createPoint(bdyPts->getAt(0)));
+        return bdyPts->applyAt(0, [this](const auto& c) {
+            return m_geomFact.createPoint(c);
+        });
     }
     // this handles 0 points case as well
     return std::unique_ptr<Geometry>(m_geomFact.createMultiPoint(*bdyPts));
@@ -148,8 +150,8 @@ BoundaryOp::boundaryMultiLineString(const geom::MultiLineString& mLine)
 std::unique_ptr<CoordinateSequence>
 BoundaryOp::computeBoundaryCoordinates(const geom::MultiLineString& mLine)
 {
-    auto bdyPts = detail::make_unique<CoordinateSequence>();
-    std::map<Coordinate, int> endpointMap;
+    auto bdyPts = detail::make_unique<CoordinateSequence>(0, mLine.hasZ(), mLine.hasM());
+    std::map<geom::CoordinateXYZM, int> endpointMap;
 
     for (std::size_t i = 0; i < mLine.getNumGeometries(); i++) {
       const LineString* line = mLine.getGeometryN(i);
@@ -158,14 +160,20 @@ BoundaryOp::computeBoundaryCoordinates(const geom::MultiLineString& mLine)
         continue;
       }
 
-      endpointMap[line->getCoordinateN(0)]++;
-      endpointMap[line->getCoordinateN(line->getNumPoints() - 1)]++;
+      const CoordinateSequence& pts = *line->getCoordinatesRO();
+
+      geom::CoordinateXYZM start;
+      geom::CoordinateXYZM end;
+      pts.getAt(0, start);
+      pts.getAt(pts.size() - 1, end);
+
+      endpointMap[start]++;
+      endpointMap[end]++;
     }
 
-    for (const auto& entry: endpointMap) {
-        auto valence = entry.second;
+    for (const auto& [coord, valence] : endpointMap) {
         if (m_bnRule.isInBoundary(valence)) {
-            bdyPts->add(entry.first);
+            bdyPts->add(coord);
         }
     }
 
diff --git a/tests/unit/capi/GEOSBoundaryTest.cpp b/tests/unit/capi/GEOSBoundaryTest.cpp
index d24bd558f3..17f4af46e3 100644
--- a/tests/unit/capi/GEOSBoundaryTest.cpp
+++ b/tests/unit/capi/GEOSBoundaryTest.cpp
@@ -41,5 +41,27 @@ void object::test<2>()
     ensure(result_ == nullptr);
 }
 
+template<>
+template<>
+void object::test<3>()
+{
+    input_ = fromWKT("MULTILINESTRING M ((10 11 3, 20 21 4), (20 21 4, 32 21 3))");
+    result_ = GEOSBoundary(input_);
+    expected_ = fromWKT("MULTIPOINT M ((10 11 3), (32 21 3))");
+
+    ensure_geometry_equals_identical(result_, expected_);
+}
+
+template<>
+template<>
+void object::test<4>()
+{
+    input_ = fromWKT("POLYGON M ((0 0 0, 1 0 1, 1 1 2, 0 1 3, 0 0 4))");
+    result_ = GEOSBoundary(input_);
+    expected_ = fromWKT("LINESTRING M (0 0 0, 1 0 1, 1 1 2, 0 1 3, 0 0 4)");
+
+    ensure_geometry_equals_identical(result_, expected_);
+}
+
 } // namespace tut
 
diff --git a/tests/unit/capi/capi_test_utils.h b/tests/unit/capi/capi_test_utils.h
index 803c51cf21..aef573565b 100644
--- a/tests/unit/capi/capi_test_utils.h
+++ b/tests/unit/capi/capi_test_utils.h
@@ -72,7 +72,10 @@ namespace capitest {
         fromWKT(const char* wkt)
         {
             GEOSGeometry* g = GEOSGeomFromWKT(wkt);
-            tut::ensure(g != nullptr);
+            if (g == nullptr) {
+                std::string message = "WKT is invalid: " + std::string(wkt);
+                tut::ensure(message, g != nullptr);
+            }
             return g;
         }