feat: update geometry_py.h, module.cpp, and tests

Author: peng.li24

pycpp/geometry_py.h

@@ -22,6 +22,7 @@
 #include "../shapely/geometry/geometrycollection.h"
 
 #include <array>
+#include <cmath>
 #include <cstring>
 #include <tuple>
 #include <vector>
@@ -66,6 +67,20 @@ inline void ensure_c_contiguous(const py::buffer_info& buf) {
     }
 }
 
+/// Validate coordinate buffer — reject NaN/Inf before they reach GEOS.
+/// Throws std::invalid_argument with the index of the first offending value.
+template <typename T>
+inline void ensure_no_nan_inf(const T* data, size_t count, const char* label = "coords") {
+    for (size_t i = 0; i < count; ++i) {
+        if (std::isnan(data[i]))
+            throw std::invalid_argument(
+                std::string(label) + " contains NaN at index " + std::to_string(i));
+        if (std::isinf(data[i]))
+            throw std::invalid_argument(
+                std::string(label) + " contains Inf at index " + std::to_string(i));
+    }
+}
+
 /// Convert any py::array coordinate buffer to std::vector<double>.
 /// Extracts the first 2 columns per row — safe for (N,2) and (N,≥3) arrays.
 /// dtype double → direct copy; float32/other → cast & widen to double.
@@ -96,6 +111,8 @@ inline std::vector<double> array_to_double_vec(const py::array& arr) {
             tmp[r * 2 + 1] = static_cast<double>(src[base + 1]);
         }
     }
+    // Validate AFTER conversion: reject NaN/Inf before GEOS
+    ensure_no_nan_inf(tmp.data(), rows * 2, "coords");
     return tmp;
 }
 
@@ -104,45 +121,65 @@ inline std::vector<double> array_to_double_vec(const py::array& arr) {
 // ============================================================================
 
 // -- Point: scalar (x,y) + array ([x,y]) forms --
-inline Point<double> point(double x, double y) { return Point<double>(x, y); }
+inline Point<double> point(double x, double y) {
+    if (std::isnan(x) || std::isnan(y))
+        throw std::invalid_argument("point coords contains NaN");
+    if (std::isinf(x) || std::isinf(y))
+        throw std::invalid_argument("point coords contains Inf");
+    return Point<double>(x, y);
+}
 inline Point<float>  point(float x, float y)  { return Point<float>(x, y); }
 
 inline Point<double> point(const py::array_t<double>& arr) {
     auto buf = arr.request();
     ensure_c_contiguous(buf);
     const double* p = static_cast<const double*>(buf.ptr);
+    ensure_no_nan_inf(p, buf.size, "point coords");
     return Point<double>(buf.size > 0 ? p[0] : 0, buf.size > 1 ? p[1] : 0);
 }
 // float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
 inline Point<double> point(const py::array_t<float>& arr) {
     auto buf = arr.request();
     ensure_c_contiguous(buf);
     const float* p = static_cast<const float*>(buf.ptr);
-    return Point<double>(buf.size > 0 ? static_cast<double>(p[0]) : 0,
-                         buf.size > 1 ? static_cast<double>(p[1]) : 0);
+    double x = buf.size > 0 ? static_cast<double>(p[0]) : 0;
+    double y = buf.size > 1 ? static_cast<double>(p[1]) : 0;
+    if (std::isnan(x) || std::isnan(y))
+        throw std::invalid_argument("point coords contains NaN");
+    if (std::isinf(x) || std::isinf(y))
+        throw std::invalid_argument("point coords contains Inf");
+    return Point<double>(x, y);
 }
 // auto-double: accept any dtype (int, unknown), always produce Point<double>
 // Handles 1D [x,y] and 2D (N,≥2) arrays — takes first 2 elements.
 inline Point<double> point(const py::array& arr) {
     auto buf = arr.request();
     ensure_c_contiguous(buf);
     if (buf.size < 2) return Point<double>(0, 0);
+    double x, y;
     if (arr.dtype().is(py::dtype::of<double>())) {
         const double* p = static_cast<const double*>(buf.ptr);
-        return Point<double>(p[0], p[1]);
+        x = p[0]; y = p[1];
     } else {
         // float32 or int/other → cast through pybind11 to float32, then widen
         auto f32 = py::cast<py::array_t<float>>(arr);
         auto fbuf = f32.request();
         const float* p = static_cast<const float*>(fbuf.ptr);
-        return Point<double>(static_cast<double>(p[0]), static_cast<double>(p[1]));
+        x = static_cast<double>(p[0]);
+        y = static_cast<double>(p[1]);
     }
+    if (std::isnan(x) || std::isnan(y))
+        throw std::invalid_argument("point coords contains NaN");
+    if (std::isinf(x) || std::isinf(y))
+        throw std::invalid_argument("point coords contains Inf");
+    return Point<double>(x, y);
 }
 
 // -- LineString --
 inline LineString<double> linestring(const py::array_t<double>& arr) {
     auto buf = arr.request();
     ensure_c_contiguous(buf);
+    ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "linestring coords");
     return LineString<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 // float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
@@ -167,13 +204,15 @@ inline LineString<double> linestring(const std::vector<std::array<double, 2>>& p
         flat[i * 2] = pts[i][0];
         flat[i * 2 + 1] = pts[i][1];
     }
+    ensure_no_nan_inf(flat.data(), pts.size() * 2, "linestring coords");
     return LineString<double>(flat.data(), pts.size(), 2);
 }
 
 // -- Polygon --
 inline Polygon<double> polygon(const py::array_t<double>& arr) {
     auto buf = arr.request();
     ensure_c_contiguous(buf);
+    ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "polygon coords");
     return Polygon<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 // float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
@@ -192,6 +231,7 @@ inline Polygon<double> polygon(const py::array& arr) {
 inline LinearRing<double> linearring(const py::array_t<double>& arr) {
     auto buf = arr.request();
     ensure_c_contiguous(buf);
+    ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "linearring coords");
     return LinearRing<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline LinearRing<double> linearring(const py::array& arr) {
@@ -204,6 +244,7 @@ inline LinearRing<double> linearring(const py::array& arr) {
 inline MultiPoint<double> multipoint(const py::array_t<double>& arr) {
     auto buf = arr.request();
     ensure_c_contiguous(buf);
+    ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "multipoint coords");
     return MultiPoint<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 // float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
@@ -224,6 +265,7 @@ inline MultiLineString<double> multilinestring(const std::vector<py::array_t<dou
     for (auto& arr : arrays) {
         auto buf = arr.request();
         ensure_c_contiguous(buf);
+        ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "multilinestring coords");
         mls.add_line(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
     }
     return mls;
@@ -254,6 +296,7 @@ inline MultiPolygon<double> multipolygon(const std::vector<py::array_t<double>>&
     for (auto& arr : arrays) {
         auto buf = arr.request();
         ensure_c_contiguous(buf);
+        ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "multipolygon coords");
         mp.add_polygon(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
     }
     return mp;
@@ -443,9 +486,8 @@ inline std::tuple<double,double> centroid_linestring(const LineString<double>& l
 inline std::tuple<double,double> centroid_polygon(const Polygon<double>& p) { auto r=p.centroid(); return {r.x,r.y}; }
 
 inline double project_ls_pt(const LineString<double>& l, const Point<double>& p) { return l.project(p); }
-inline std::tuple<double,double> interpolate_ls(const LineString<double>& l, double d) { auto r=l.interpolate(d); return {r.x,r.y}; }
-inline std::tuple<double,double> interpolate_ls_normalized(const LineString<double>& l, double t) {
-    auto r = l.interpolate(t, true);
+inline std::tuple<double,double> interpolate_ls(const LineString<double>& l, double d, bool normalized = false) {
+    auto r = l.interpolate(d, normalized);
     return {r.x, r.y};
 }
 

tests/module.cpp

@@ -388,8 +388,8 @@ PYBIND11_MODULE(shapelycpp, m) {
     // Project / interpolate (from pycpp)
     // ======================================================================
     m.def("project_linestring_point",   &project_ls_pt);
-    m.def("interpolate_linestring",    &interpolate_ls);
-    m.def("interpolate_linestring_normalized", &interpolate_ls_normalized);
+    m.def("interpolate_linestring",    &interpolate_ls,
+          py::arg("ls"), py::arg("distance"), py::arg("normalized") = false);
 
     // ======================================================================
     // Intersection area (from pycpp)

tests/report.csv

@@ -230,11 +230,11 @@ ops,interpolate_linestring,unknown,f64,v_mid,PASS,0
 ops,interpolate_linestring,unknown,f64,h_start,PASS,0
 ops,interpolate_linestring,unknown,f64,h_end,PASS,0
 ops,interpolate_linestring,unknown,f64,v3pt,PASS,0
-ops,interpolate_linestring_normalized,unknown,f64,h_mid,PASS,0
-ops,interpolate_linestring_normalized,unknown,f64,h_start,PASS,0
-ops,interpolate_linestring_normalized,unknown,f64,h_end,PASS,0
-ops,interpolate_linestring_normalized,unknown,f64,v_mid,PASS,0
-ops,interpolate_linestring_normalized,unknown,f64,v3pt_end,PASS,0
+ops,interpolate_linestring,unknown,f64,norm_h_mid,PASS,0
+ops,interpolate_linestring,unknown,f64,norm_h_start,PASS,0
+ops,interpolate_linestring,unknown,f64,norm_h_end,PASS,0
+ops,interpolate_linestring,unknown,f64,norm_v_mid,PASS,0
+ops,interpolate_linestring,unknown,f64,norm_v3pt_end,PASS,0
 ops,centroid_point,unknown,f64,3-4,PASS,0
 ops,centroid_point,unknown,f64,neg,PASS,0
 ops,centroid_linestring,unknown,f64,d,PASS,0
@@ -249,3 +249,15 @@ ops,nearest_points,unknown,f64,poly_ls_touch,PASS,0
 ops,nearest_points,unknown,f64,poly_ls_near,PASS,0
 ops,nearest_points_ls_pt,unknown,f64,ls_pt_above,PASS,0
 ops,nearest_points_ls_pt,unknown,f64,ls_pt_endpoint,PASS,0
+errors,point,unknown,f64,nan,PASS,0
+errors,point,unknown,f64,inf,PASS,0
+errors,linestring,unknown,f64,nan,PASS,0
+errors,linestring,unknown,f64,inf,PASS,0
+errors,linestring,unknown,f64,ninf,PASS,0
+errors,polygon,unknown,f64,nan,PASS,0
+errors,polygon,unknown,f64,inf,PASS,0
+errors,linearring,unknown,f64,nan,PASS,0
+errors,linearring,unknown,f64,inf,PASS,0
+errors,multipoint,unknown,f64,nan,PASS,0
+errors,multilinestring,unknown,f64,nan,PASS,0
+errors,multipolygon,unknown,f64,nan,PASS,0

tests/test_all.py

@@ -147,7 +147,6 @@ def _compare_bit_exact(cpp_result, py_result, label=""):
 _OPS_NAMES = {
     "centroid_point", "centroid_linestring", "centroid_polygon",
     "project_linestring_point", "interpolate_linestring",
-    "interpolate_linestring_normalized",
     "intersection_area_polygon_polygon",
     "polygon_exterior",
     "nearest_points", "nearest_points_ls_pt",
@@ -353,16 +352,12 @@ def _call_ops(api_name, cpp, *args, **kwargs):
         c_ls = _build_cpp_geom(cpp, *geoms[0])
         py_ls = _build_py_geom(*geoms[0])
         dist = float(extras[0]) if extras else 5.0
-        c_result = cpp.interpolate_linestring(c_ls, dist)
-        py_pt = py_ls.interpolate(dist)
-        return np.array(c_result), np.array((py_pt.x, py_pt.y))
-
-    if api_name == "interpolate_linestring_normalized":
-        c_ls = _build_cpp_geom(cpp, *geoms[0])
-        py_ls = _build_py_geom(*geoms[0])
-        t = float(extras[0]) if extras else 0.5
-        c_result = cpp.interpolate_linestring_normalized(c_ls, t)
-        py_pt = py_ls.interpolate(t, normalized=True)
+        normalized = bool(kwargs.get("normalized", False))
+        if normalized:
+            c_result = cpp.interpolate_linestring(c_ls, dist, normalized=True)
+        else:
+            c_result = cpp.interpolate_linestring(c_ls, dist)
+        py_pt = py_ls.interpolate(dist, normalized=normalized)
         return np.array(c_result), np.array((py_pt.x, py_pt.y))
 
     if api_name == "intersection_area_polygon_polygon":
@@ -534,6 +529,19 @@ def _call_class_method(api_name, cpp, *args, **kwargs):
 _I_A   = [(0, 0), (5, 0), (5, 5), (0, 5)]
 _I_B   = [(3, 0), (8, 0), (8, 5), (3, 5)]
 
+# ---- NaN/Inf test data presets ------------------------------------------------
+
+_NAN_1 = [(float('nan'), 0), (1, 1)]
+_INF_1 = [(float('inf'), 0), (1, 1)]
+_NINF_1= [(float('-inf'), 0), (1, 1)]
+_NAN_PT = (float('nan'), 0)
+_INF_PT = (float('inf'), 0)
+
+_MNAN  = [(float('nan'), 0), (10, 10), (20, 5)]        # multipoint nan
+_MLNAN = [[(float('nan'), 0), (5, 0)], [(10, 0), (15, 0)]]  # multiline nan
+_MPNAN = [[(0, 0), (10, 0), (10, 10), (0, 10)],      # multipoly nan
+          [(float('nan'), 0), (20, 0), (20, 10), (10, 10)]]
+
 
 # ---- Group 1: factories (7 APIs) --------------------------------------------
 
@@ -877,7 +885,7 @@ def _catalog_ops():
         yield _TestCase("project_linestring_point", (ls, pt), {}, "f64", tag,
                        "scalar_eq", True, None, "ops")
 
-    # -- interpolate_linestring --
+    # -- interpolate_linestring (absolute distance) --
     for ls, dist, tag in [
         (("ls", _L_H), 5.0, "h_mid"),
         (("ls", _L_V), 5.0, "v_mid"),
@@ -888,15 +896,15 @@ def _catalog_ops():
         yield _TestCase("interpolate_linestring", (ls, dist), {}, "f64", tag,
                        "bit_exact", True, None, "ops")
 
-    # -- interpolate_linestring_normalized --
-    for ls, t, py_t_expected_tag, tag in [
-        (("ls", _L_H), 0.5, "h_mid", "h_mid"),       # t=0.5 → midpoint of 10-length horizontal
-        (("ls", _L_H), 0.0, "h_start", "h_start"),   # t=0.0 → start
-        (("ls", _L_H), 1.0, "h_end", "h_end"),        # t=1.0 → end
-        (("ls", _L_V), 0.5, "v_mid", "v_mid"),        # vertical 10-length, t=0.5
-        (("ls", _L_3), 1.0, "v3pt", "v3pt_end"),      # 3-point L, t=1.0
+    # -- interpolate_linestring (normalized) -- same API, kwargs={"normalized": True}
+    for ls, t, tag in [
+        (("ls", _L_H), 0.5, "norm_h_mid"),
+        (("ls", _L_H), 0.0, "norm_h_start"),
+        (("ls", _L_H), 1.0, "norm_h_end"),
+        (("ls", _L_V), 0.5, "norm_v_mid"),
+        (("ls", _L_3), 1.0, "norm_v3pt_end"),
     ]:
-        yield _TestCase("interpolate_linestring_normalized", (ls, t), {},
+        yield _TestCase("interpolate_linestring", (ls, t), {"normalized": True},
                        "f64", tag, "bit_exact", True, None, "ops")
 
     # -- Centroid --
@@ -955,6 +963,37 @@ def api_catalog():
         yield tc._replace(group="predicates")
     for tc in _catalog_ops():
         yield tc._replace(group="ops")
+    for tc in _catalog_errors():
+        yield tc._replace(group="errors")
+
+
+def _catalog_errors():
+    """NaN/Inf rejection: every geometry factory must reject NaN/Inf coords."""
+    # point
+    for coords, tag in [(_NAN_PT, "nan"), (_INF_PT, "inf")]:
+        yield _TestCase("point", (("pt", coords),), {}, "f64", tag,
+                       "expect_error", True, None, "errors")
+    # linestring
+    for coords, tag in [(_NAN_1, "nan"), (_INF_1, "inf"), (_NINF_1, "ninf")]:
+        yield _TestCase("linestring", (coords,), {}, "f64", tag,
+                       "expect_error", True, None, "errors")
+    # polygon
+    for coords, tag in [(_NAN_1, "nan"), (_INF_1, "inf")]:
+        yield _TestCase("polygon", (coords,), {}, "f64", tag,
+                       "expect_error", True, None, "errors")
+    # linearring
+    for coords, tag in [(_NAN_1, "nan"), (_INF_1, "inf")]:
+        yield _TestCase("linearring", (coords,), {}, "f64", tag,
+                       "expect_error", True, None, "errors")
+    # multipoint
+    yield _TestCase("multipoint", (_MNAN,), {}, "f64", "nan",
+                   "expect_error", True, None, "errors")
+    # multilinestring
+    yield _TestCase("multilinestring", (_MLNAN,), {}, "f64", "nan",
+                   "expect_error", True, None, "errors")
+    # multipolygon
+    yield _TestCase("multipolygon", (_MPNAN,), {}, "f64", "nan",
+                   "expect_error", True, None, "errors")
 
 
 # =============================================================================
@@ -1028,6 +1067,28 @@ def test_api(tc, cpp):
         tc.setup_fn(args, kwargs)
         return
 
+    # Error tests: expect exception from C++
+    if strategy == "expect_error":
+        try:
+            call_cpp_py(api_name, cpp, *args, **kwargs)
+            # Should have raised — didn't
+            raise AssertionError(
+                f"[{api_name}] expected exception for NaN/Inf input, but no error raised")
+        except Exception as e:
+            # Re-raise assertion errors, accept everything else
+            if isinstance(e, AssertionError):
+                _REPORT_ROWS.append(dict(
+                    category=tc.group, api=api_name, mode=_geos_info(cpp),
+                    dtype=tc.dtype_label, feature=tc.category,
+                    result="FAIL", ulp=-1))
+                raise
+            # Expected: any exception (invalid_argument, runtime_error, etc.)
+            _REPORT_ROWS.append(dict(
+                category=tc.group, api=api_name, mode=_geos_info(cpp),
+                dtype=tc.dtype_label, feature=tc.category,
+                result="PASS", ulp=0))
+        return
+
     cpp_r, py_r = call_cpp_py(api_name, cpp, *args, **kwargs)
     if check_py and py_r is None:
         pytest.skip(f"no shapely equivalent for {api_name}")