round rapidjson, round z in geobuf (#17)

* add round_z

* fix Encoder default constructor

* minor fix

* round rapidjson

* update

* fix windows

Author: district10

docs/about/release-notes.md

@@ -10,6 +10,11 @@ To upgrade `pybind11-geobuf` to the latest version, use pip:
 pip install -U pybind11-geobuf
 ```
 
+## Version 0.1.0 (2023-04-01)
+
+*   round z in geobuf encoding
+*   round rapidjson
+
 ## Version 0.0.9 (2023-03-30)
 
 *   Export to geobuf with only xy (no z), `only_xy=False`

setup.py

@@ -127,7 +127,7 @@ def build_extension(self, ext):
 # logic and declaration, and simpler if you include description/version in a file.
 setup(
     name="pybind11_geobuf",
-    version="0.0.9",
+    version="0.1.0",
     author="tzx",
     author_email="dvorak4tzx@gmail.com",
     url="https://geobuf-cpp.readthedocs.io",

src/geobuf/geobuf.cpp

@@ -221,7 +221,7 @@ std::string Encoder::encode(const mapbox::geojson::geojson &geojson)
     keys.clear();
     analyze(geojson);
     if (onlyXY) {
-        dim = 2u;
+        dim = dimXY;
     }
 
     {
@@ -357,7 +357,7 @@ inline double ROUND(double v, double s)
 
 void Encoder::analyzePoint(const mapbox::geojson::point &point)
 {
-    if (!onlyXY) {
+    if (!onlyXY && dim < dimXYZ) {
         dim = std::max(point.z == 0 ? dimXY : dimXYZ, dim);
     }
     if (e >= maxPrecision) {
@@ -509,13 +509,14 @@ void Encoder::writeValue(const mapbox::feature::value &value, Encoder::Pbf &pbf)
 
 void Encoder::writePoint(const mapbox::geojson::point &point, Encoder::Pbf &pbf)
 {
-    std::vector<int64_t> coords;
-    coords.reserve(dim);
-    const double *ptr = &point.x;
-    for (int i = 0; i < dim; ++i) {
-        coords.push_back(static_cast<int64_t>(std::floor(ptr[i] * e + 0.5)));
+    std::array<int64_t, 3> coords;
+    coords[0] = static_cast<int64_t>(std::floor(point.x * e + 0.5));
+    coords[1] = static_cast<int64_t>(std::floor(point.y * e + 0.5));
+    if (dim == 3) {
+        double z = zScale ? ROUND(point.z, *zScale) : point.z;
+        coords[2] = static_cast<int64_t>(std::floor(z * e + 0.5));
     }
-    pbf.add_packed_sint64(3, coords.begin(), coords.end());
+    pbf.add_packed_sint64(3, &coords[0], &coords[0] + dim);
 }
 
 void Encoder::writeLine(const PointsType &line, Encoder::Pbf &pbf)
@@ -579,12 +580,18 @@ void Encoder::populateLine(std::vector<int64_t> &coords, //
     int len = line.size() - (closed ? 1 : 0);
     auto sum = std::array<int64_t, 3>{0, 0, 0};
     for (int i = 0; i < len; ++i) {
-        const double *ptr = &line[i].x;
-        for (int j = 0; j < dim; ++j) {
-            auto n =
-                static_cast<int64_t>(std::floor(ptr[j] * e + 0.5)) - sum[j];
+        const auto &pt = line[i];
+        auto n = static_cast<int64_t>(std::floor(pt.x * e + 0.5)) - sum[0];
+        coords.push_back(n);
+        sum[0] += n;
+        n = static_cast<int64_t>(std::floor(pt.y * e + 0.5)) - sum[1];
+        coords.push_back(n);
+        sum[1] += n;
+        if (dim == 3) {
+            double z = zScale ? ROUND(pt.z, *zScale) : pt.z;
+            n = static_cast<int64_t>(std::floor(z * e + 0.5)) - sum[2];
             coords.push_back(n);
-            sum[j] += n;
+            sum[2] += n;
         }
     }
 }

src/geobuf/geobuf.hpp

@@ -7,6 +7,7 @@
 #include <protozero/pbf_reader.hpp>
 
 #include <map>
+#include <optional>
 #include <unordered_map>
 
 #define MAPBOX_GEOBUF_DEFAULT_PRECISION 6
@@ -81,8 +82,11 @@ struct Encoder
     using Pbf = protozero::pbf_writer;
     Encoder(uint32_t maxPrecision = std::pow(10,
                                              MAPBOX_GEOBUF_DEFAULT_PRECISION),
-            bool onlyXY = false)
-        : maxPrecision(maxPrecision), onlyXY(onlyXY)
+            bool onlyXY = false, //
+            std::optional<int> roundZ = std::nullopt)
+        : maxPrecision(maxPrecision), onlyXY(onlyXY),
+          zScale(roundZ ? std::make_optional(std::pow(10, *roundZ))
+                        : std::nullopt)
     {
     }
     std::string encode(const mapbox::geojson::geojson &geojson);
@@ -103,6 +107,14 @@ struct Encoder
     std::string encode(const RapidjsonValue &json);
     bool encode(const std::string &input_path, const std::string &output_path);
 
+    auto __maxPrecision() const { return maxPrecision; }
+    auto __onlyXY() const { return onlyXY; }
+    auto __roundZ() const
+    {
+        return zScale ? std::make_optional(std::log10(*zScale)) : std::nullopt;
+    }
+    auto __dim() const { return dim; }
+    auto __e() const { return e; }
     std::map<std::string, std::uint32_t> __keys() const
     {
         return std::map<std::string, std::uint32_t>(keys.begin(), keys.end());
@@ -140,6 +152,7 @@ struct Encoder
 
     const uint32_t maxPrecision;
     const bool onlyXY;
+    const std::optional<double> zScale;
     uint32_t dim = MAPBOX_GEOBUF_DEFAULT_DIM;
     uint32_t e = 1;
     std::unordered_map<std::string, std::uint32_t> keys;

src/geobuf/rapidjson_helpers.hpp

@@ -68,6 +68,36 @@ inline void sort_keys_inplace(RapidjsonValue &json)
     }
 }
 
+inline void round_rapidjson(RapidjsonValue &json, double scale, int depth = 1,
+                            const std::vector<std::string> &skip_keys = {})
+{
+    if (--depth < 0) {
+        return;
+    }
+    if (json.IsArray()) {
+        for (auto &e : json.GetArray()) {
+            round_rapidjson(e, scale, depth, skip_keys);
+        }
+    } else if (json.IsObject()) {
+        auto obj = json.GetObject();
+        for (auto &kv : obj) {
+            if (!skip_keys.empty() &&
+                std::find(skip_keys.begin(), skip_keys.end(),
+                          std::string(kv.name.GetString(),
+                                      kv.name.GetStringLength())) !=
+                    skip_keys.end()) {
+                continue;
+            }
+            round_rapidjson(kv.value, scale, depth, skip_keys);
+        }
+    } else if (json.IsDouble()) {
+        // see round_coords in geojson_helpers
+        json.SetDouble(std::floor(json.GetDouble() * scale + 0.5) / scale);
+    } else if (json.IsFloat()) {
+        json.SetFloat(std::floor(json.GetFloat() * scale + 0.5) / scale);
+    }
+}
+
 inline RapidjsonValue sort_keys(const RapidjsonValue &json)
 {
     RapidjsonAllocator allocator;

src/main.cpp

@@ -91,12 +91,20 @@ PYBIND11_MODULE(_pybind11_geobuf, m)
         py::kw_only(), //
         "indent"_a = "");
 
-    py::class_<Encoder>(m, "Encoder", py::module_local()) //
-        .def(py::init<uint32_t, bool>(),                  //
+    py::class_<Encoder>(m, "Encoder", py::module_local())    //
+        .def(py::init<uint32_t, bool, std::optional<int>>(), //
              py::kw_only(),
-             "max_precision"_a = std::pow(10, MAPBOX_GEOBUF_DEFAULT_PRECISION),
-             "only_xy"_a = false)
+             "max_precision"_a = static_cast<uint32_t>(
+                 std::pow(10, MAPBOX_GEOBUF_DEFAULT_PRECISION)),
+             "only_xy"_a = false, //
+             "round_z"_a = std::nullopt)
         //
+        .def("max_precision", &Encoder::__maxPrecision)
+        .def("only_xy", &Encoder::__onlyXY)
+        .def("round_z", &Encoder::__roundZ)
+        .def("dim", &Encoder::__dim)
+        .def("e", &Encoder::__e)
+        .def("keys", &Encoder::__keys)
         .def(
             "encode",
             [](Encoder &self, const mapbox::geojson::geojson &geojson) {

src/pybind11_geojson.cpp

@@ -134,14 +134,17 @@ void bind_geojson(py::module &geojson)
                 rvp::reference_internal)
             .def(
                 "to_geobuf",
-                [](const mapbox::geojson::geojson &self, int precision,
-                   bool only_xy) {
+                [](const mapbox::geojson::geojson &self, //
+                   int precision, bool only_xy, std::optional<int> round_z) {
                     auto bytes = mapbox::geobuf::Encoder(
-                                     std::pow(10, precision), only_xy)
+                                     std::pow(10, precision), only_xy, round_z)
                                      .encode(self);
                     return py::bytes(bytes);
                 },
-                py::kw_only(), "precision"_a = 8, "only_xy"_a = false)
+                py::kw_only(),       //
+                "precision"_a = 8,   //
+                "only_xy"_a = false, //
+                "round_z"_a = std::nullopt)
             //
             .def(
                 "load",
@@ -503,12 +506,16 @@ void bind_geojson(py::module &geojson)
             },
             rvp::reference_internal)
         .def("to_geobuf",
-            [](const mapbox::geojson::geometry &self, int precision, bool only_xy) {
+            [](const mapbox::geojson::geometry &self, //
+            int precision, bool only_xy, std::optional<int> round_z) {
             auto bytes =
-                mapbox::geobuf::Encoder(std::pow(10, precision), only_xy)
+                mapbox::geobuf::Encoder(std::pow(10, precision), only_xy, round_z)
                     .encode(self);
             return py::bytes(bytes);
-            }, py::kw_only(), "precision"_a = 8, "only_xy"_a = false)
+            }, py::kw_only(), //
+            "precision"_a = 8, //
+            "only_xy"_a = false, //
+            "round_z"_a = std::nullopt)
         .def("load",
             [](mapbox::geojson::geometry &self, const std::string &path) -> mapbox::geojson::geometry & {
                 if (endswith(path, ".pbf")) {
@@ -1763,14 +1770,17 @@ void bind_geojson(py::module &geojson)
             rvp::reference_internal)
         .def(
             "to_geobuf",
-            [](const mapbox::geojson::feature &self, int precision,
-               bool only_xy) {
-                auto bytes =
-                    mapbox::geobuf::Encoder(std::pow(10, precision), only_xy)
-                        .encode(self);
+            [](const mapbox::geojson::feature &self, //
+               int precision, bool only_xy, std::optional<int> round_z) {
+                auto bytes = mapbox::geobuf::Encoder(std::pow(10, precision),
+                                                     only_xy, round_z)
+                                 .encode(self);
                 return py::bytes(bytes);
             },
-            py::kw_only(), "precision"_a = 8, "only_xy"_a = false)
+            py::kw_only(),       //
+            "precision"_a = 8,   //
+            "only_xy"_a = false, //
+            "round_z"_a = std::nullopt)
         .def(
             "load",
             [](mapbox::geojson::feature &self,
@@ -1902,14 +1912,17 @@ void bind_geojson(py::module &geojson)
             rvp::reference_internal)
         .def(
             "to_geobuf",
-            [](const mapbox::geojson::feature_collection &self, int precision,
-               bool only_xy) {
-                auto bytes =
-                    mapbox::geobuf::Encoder(std::pow(10, precision), only_xy)
-                        .encode(self);
+            [](const mapbox::geojson::feature_collection &self, //
+               int precision, bool only_xy, std::optional<int> round_z) {
+                auto bytes = mapbox::geobuf::Encoder(std::pow(10, precision),
+                                                     only_xy, round_z)
+                                 .encode(self);
                 return py::bytes(bytes);
             },
-            py::kw_only(), "precision"_a = 8, "only_xy"_a = false)
+            py::kw_only(),       //
+            "precision"_a = 8,   //
+            "only_xy"_a = false, //
+            "round_z"_a = std::nullopt)
         .def(
             "load",
             [](mapbox::geojson::feature_collection &self,

src/pybind11_rapidjson.cpp

@@ -145,6 +145,14 @@ void bind_rapidjson(py::module &m)
                 sort_keys_inplace(self);
                 return self;
             }, rvp::reference_internal)
+            .def("round", [](RapidjsonValue &self, double precision, int depth, //
+                const std::vector<std::string> &skip_keys) -> RapidjsonValue & {
+                    round_rapidjson(self, std::pow(10, precision), depth, skip_keys);
+                return self;
+            }, rvp::reference_internal, py::kw_only(), //
+                "precision"_a = 3, //
+                "depth"_a = 32, //
+                "skip_keys"_a = std::vector<std::string>{})
             .def(
                 "get",
                 [](RapidjsonValue &self,

tests/test_geobuf.py

@@ -62,16 +62,44 @@ def test_geobuf():
     print(str2geojson2str(json.dumps(geojson), indent=True, sort_keys=True))
 
     # precision: 6(default), 7, 8(recommand), 9
+    encoder = Encoder()
+    assert encoder.max_precision() == 10**6
+    assert not encoder.only_xy()
+    assert encoder.round_z() is None
+
     encoder = Encoder(max_precision=int(10**8))
+    assert encoder.max_precision() == 10**8
     encoded = encoder.encode(geojson=json.dumps(geojson))
     print("encoded pbf bytes")
     print(pbf_decode(encoded))
 
+    encoder = Encoder(round_z=3)
+    assert encoder.round_z() == 3
+
+    with pytest.raises(Exception) as excinfo:
+        Encoder(max_precision=int(10**10))  # uint32_t overflow
+    assert "incompatible constructor arguments" in repr(excinfo)
+
     decoder = Decoder()
     geojson_text = decoder.decode(encoded, indent=True)
     print(geojson_text)
 
 
+def test_geobuf_roundz():
+    f = geojson.Feature().from_rapidjson(sample_geojson())
+    llas0 = geojson.Feature().from_geobuf(f.to_geobuf()).to_numpy()
+
+    llas = geojson.Feature().from_geobuf(f.to_geobuf(precision=3)).to_numpy()
+    assert np.all(np.round(llas, 3) == llas)
+
+    llas = geojson.Feature().from_geobuf(f.to_geobuf(only_xy=True)).to_numpy()
+    assert np.all(llas[:, 2] == np.zeros(len(llas)))
+
+    llas = geojson.Feature().from_geobuf(f.to_geobuf(round_z=3)).to_numpy()
+    assert np.all(np.round(llas0[:, :2], 8) == llas[:, :2])
+    assert np.all(np.round(llas0[:, 2], 3) == llas[:, 2])
+
+
 def test_rapidjson_empty():
     assert not bool(rapidjson())
     assert not bool(rapidjson([]))
@@ -247,6 +275,28 @@ def test_rapidjson_sort_dump():
     assert obj6 == obj5
 
 
+def test_rapidjson_round():
+    arr = rapidjson([1.23456, [3.2, 2.345]])
+    assert arr() == [1.23456, [3.2, 2.345]]
+    arr.round()
+    assert arr() == [1.235, [3.2, 2.345]]
+    arr.round(precision=2)
+    assert arr() == [1.24, [3.2, 2.35]]
+
+    obj = rapidjson(sample_geojson())
+    assert obj() == sample_geojson()
+    assert (
+        obj.clone().round(precision=1).dumps(sort_keys=True)
+        == '{"geometry":{"coordinates":[[120.4,31.4,1.1],[120.3,31.3,2.2],[120.4,31.2,3.3],[120.7,31.3,4.4]],"extra_key":"extra_value","type":"LineString"},"my_key":"my_value","properties":{"dict":{"key":42,"value":3.1},"double":3.1,"int":42,"int2":-101,"list":["a","list","is","a","list"],"string":"string"},"type":"Feature"}'  # noqa
+    )
+    assert (
+        obj.clone()
+        .round(precision=1, skip_keys=["geometry"])
+        .dumps(sort_keys=True)  # noqa
+        == '{"geometry":{"coordinates":[[120.40317479950272,31.416966084052177,1.111111],[120.28451900911591,31.30578266928819,2.22],[120.35592249359615,31.21781895672254,3.3333333333333],[120.67093786630113,31.299502266522722,4.4]],"extra_key":"extra_value","type":"LineString"},"my_key":"my_value","properties":{"dict":{"key":42,"value":3.1},"double":3.1,"int":42,"int2":-101,"list":["a","list","is","a","list"],"string":"string"},"type":"Feature"}'  # noqa
+    )
+
+
 def test_geojson_point():
     # as_numpy
     g1 = geojson.Point()