From ae001637f299b4036b008cd201279ed5dcf3e6e5 Mon Sep 17 00:00:00 2001 From: Joris Van den Bossche Date: Mon, 11 Nov 2024 10:21:54 +0100 Subject: [PATCH 1/3] ENH: add projection support in geoarrow output --- CMakeLists.txt | 2 +- src/geoarrow.cpp | 10 ++++++-- src/projections.cpp | 35 ++++++++++++++++++++++++++++ src/projections.hpp | 38 ++++++++++++++++++++++++++++++ src/spherely.cpp | 2 ++ tests/test_geoarrow.py | 52 ++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 136 insertions(+), 3 deletions(-) create mode 100644 src/projections.cpp create mode 100644 src/projections.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 326fa4f..a9e0bb1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -80,7 +80,7 @@ set(CPP_SOURCES src/spherely.cpp) if(${s2geography_VERSION} VERSION_GREATER_EQUAL "0.2.0") - set(CPP_SOURCES ${CPP_SOURCES} src/geoarrow.cpp) + set(CPP_SOURCES ${CPP_SOURCES} src/geoarrow.cpp src/projections.cpp) endif() add_library(spherely MODULE ${CPP_SOURCES}) diff --git a/src/geoarrow.cpp b/src/geoarrow.cpp index 1b7cc21..8d56b56 100644 --- a/src/geoarrow.cpp +++ b/src/geoarrow.cpp @@ -4,6 +4,7 @@ #include "constants.hpp" #include "creation.hpp" #include "geography.hpp" +#include "projections.hpp" #include "pybind11.hpp" namespace py = pybind11; @@ -197,6 +198,7 @@ class ArrowArrayHolder { ArrowArrayHolder to_geoarrow(py::array_t input, py::object output_schema, + Projection projection, bool planar, float tessellate_tolerance, int precision) { @@ -207,11 +209,10 @@ ArrowArrayHolder to_geoarrow(py::array_t input, s2geog::geoarrow::ExportOptions options; options.set_precision(precision); + options.set_projection(projection.s2_projection()); if (planar) { auto tol = S1Angle::Radians(tessellate_tolerance / EARTH_RADIUS_METERS); options.set_tessellate_tolerance(tol); - // TODO make this configurable - // options.set_projection(s2geog::geoarrow::mercator()); } if (!output_schema.is(py::none())) { @@ -315,6 +316,7 @@ void init_geoarrow(py::module& m) { py::pos_only(), py::kw_only(), py::arg("output_schema") = py::none(), + py::arg("projection") = Projection::lnglat(), py::arg("planar") = false, py::arg("tessellate_tolerance") = 100.0, py::arg("precision") = 6, @@ -338,6 +340,10 @@ void init_geoarrow(py::module& m) { can be a ``pyarrow.DataType`` or ``pyarrow.Field``, and you can use the ``geoarrow.pyarrow`` package to construct such geoarrow extension types. + projection : spherely.Projection, default Projection.lnglat() + The projection to use when converting the geographies to the output + encoding. By default, uses longitude/latitude coordinates ("plate + carree" projection). planar : bool, default False If set to True, the edges of linestrings and polygons in the output are assumed to be linear on the plane. In that case, additional diff --git a/src/projections.cpp b/src/projections.cpp new file mode 100644 index 0000000..9219ef5 --- /dev/null +++ b/src/projections.cpp @@ -0,0 +1,35 @@ +#include "projections.hpp" + +#include +#include + +#include "geography.hpp" +#include "pybind11.hpp" + +namespace py = pybind11; +namespace s2geog = s2geography; +using namespace spherely; + +// std::tuple project_mercator(py::object obj) { +// auto geog = (static_cast(obj)).as_geog_ptr(); +// if (geog->geog_type() != GeographyType::Point) { +// throw py::value_error("test"); +// } +// auto point = static_cast(geog); +// auto s2point = point->s2point(); + +// auto projection = S2::MercatorProjection(20037508.3427892); +// auto r2point = projection.Project(s2point); +// double x = r2point.x(); +// double y = r2point.y(); + +// return std::make_tuple(x, y); +// } + +void init_projections(py::module& m) { + // m.def("project_mercator", &project_mercator); + py::class_(m, "Projection") + .def("lnglat", &Projection::lnglat) + .def("pseudo_mercator", &Projection::pseudo_mercator) + .def("orthographic", &Projection::orthographic); +} diff --git a/src/projections.hpp b/src/projections.hpp new file mode 100644 index 0000000..7989053 --- /dev/null +++ b/src/projections.hpp @@ -0,0 +1,38 @@ +#ifndef SPHERELY_PROJECTIONS_H_ +#define SPHERELY_PROJECTIONS_H_ + +#include +#include +#include + +#include "pybind11.hpp" + +namespace py = pybind11; +namespace s2geog = s2geography; +using namespace spherely; + +class Projection { +public: + Projection(std::shared_ptr projection) + : m_s2_projection(std::move(projection)) {} + + std::shared_ptr s2_projection() { + return m_s2_projection; + } + + static Projection lnglat() { + return Projection(std::move(s2geog::lnglat())); + } + static Projection pseudo_mercator() { + return Projection(std::move(s2geog::pseudo_mercator())); + } + static Projection orthographic(double longitude, double latitude) { + return Projection( + std::move(s2geog::orthographic(S2LatLng::FromDegrees(latitude, longitude)))); + } + +private: + std::shared_ptr m_s2_projection; +}; + +#endif // SPHERELY_PROJECTIONS_H_ diff --git a/src/spherely.cpp b/src/spherely.cpp index 16c14ef..bd0bb66 100644 --- a/src/spherely.cpp +++ b/src/spherely.cpp @@ -13,6 +13,7 @@ void init_io(py::module&); #if defined(S2GEOGRAPHY_VERSION_MAJOR) && \ (S2GEOGRAPHY_VERSION_MAJOR >= 1 || S2GEOGRAPHY_VERSION_MINOR >= 2) void init_geoarrow(py::module&); +void init_projections(py::module&); #endif PYBIND11_MODULE(spherely, m) { @@ -31,6 +32,7 @@ PYBIND11_MODULE(spherely, m) { init_io(m); #if defined(S2GEOGRAPHY_VERSION_MAJOR) && \ (S2GEOGRAPHY_VERSION_MAJOR >= 1 || S2GEOGRAPHY_VERSION_MINOR >= 2) + init_projections(m); init_geoarrow(m); #endif diff --git a/tests/test_geoarrow.py b/tests/test_geoarrow.py index f3b8a39..31b0027 100644 --- a/tests/test_geoarrow.py +++ b/tests/test_geoarrow.py @@ -182,3 +182,55 @@ def test_to_geoarrow_invalid_output_schema(): with pytest.raises(ValueError, match="Did you pass a valid schema"): spherely.to_geoarrow(arr, output_schema=pa.schema([("test", pa.int64())])) + + +def test_to_geoarrow_projected(): + arr = spherely.points([1, 2, 3], [1, 2, 3]) + point_schema = ga.point().with_coord_type(ga.CoordType.INTERLEAVED) + result = pa.array( + spherely.to_geoarrow( + arr, output_schema=point_schema, projection=spherely.Projection.lnglat() + ) + ) + + coords = np.asarray(result.storage.values) + expected = np.array([1, 1, 2, 2, 3, 3], dtype="float64") + np.testing.assert_allclose(coords, expected) + + # import pyproj + # trans = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) + # trans.transform([1, 2, 3], [1, 2, 3]) + result = pa.array( + spherely.to_geoarrow( + arr, + output_schema=point_schema, + projection=spherely.Projection.pseudo_mercator(), + ) + ) + coords = np.asarray(result.storage.values) + expected = np.array( + [ + 111319.49079327357, + 111325.1428663851, + 222638.98158654713, + 222684.20850554405, + 333958.4723798207, + 334111.1714019596, + ], + dtype="float64", + ) + np.testing.assert_allclose(coords, expected) + + result = pa.array( + spherely.to_geoarrow( + arr, + output_schema=point_schema, + projection=spherely.Projection.orthographic(0.0, 0.0), + ) + ) + coords = np.asarray(result.storage.values) + expected = np.array( + [0.01744975, 0.01745241, 0.03487824, 0.0348995, 0.05226423, 0.05233596], + dtype="float64", + ) + np.testing.assert_allclose(coords, expected, rtol=1e-06) From 5d819fbeb7390a247d2e18ea8160e77d89fe3e75 Mon Sep 17 00:00:00 2001 From: Joris Van den Bossche Date: Fri, 15 Nov 2024 11:41:31 +0100 Subject: [PATCH 2/3] clean-up --- src/projections.cpp | 17 ----------------- tests/test_geoarrow.py | 8 +++++--- 2 files changed, 5 insertions(+), 20 deletions(-) diff --git a/src/projections.cpp b/src/projections.cpp index 9219ef5..9c2078b 100644 --- a/src/projections.cpp +++ b/src/projections.cpp @@ -10,24 +10,7 @@ namespace py = pybind11; namespace s2geog = s2geography; using namespace spherely; -// std::tuple project_mercator(py::object obj) { -// auto geog = (static_cast(obj)).as_geog_ptr(); -// if (geog->geog_type() != GeographyType::Point) { -// throw py::value_error("test"); -// } -// auto point = static_cast(geog); -// auto s2point = point->s2point(); - -// auto projection = S2::MercatorProjection(20037508.3427892); -// auto r2point = projection.Project(s2point); -// double x = r2point.x(); -// double y = r2point.y(); - -// return std::make_tuple(x, y); -// } - void init_projections(py::module& m) { - // m.def("project_mercator", &project_mercator); py::class_(m, "Projection") .def("lnglat", &Projection::lnglat) .def("pseudo_mercator", &Projection::pseudo_mercator) diff --git a/tests/test_geoarrow.py b/tests/test_geoarrow.py index 31b0027..d7706b4 100644 --- a/tests/test_geoarrow.py +++ b/tests/test_geoarrow.py @@ -197,9 +197,10 @@ def test_to_geoarrow_projected(): expected = np.array([1, 1, 2, 2, 3, 3], dtype="float64") np.testing.assert_allclose(coords, expected) - # import pyproj - # trans = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) - # trans.transform([1, 2, 3], [1, 2, 3]) + # Output to pseudo mercator - generation of expected result + # import pyproj + # trans = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) + # trans.transform([1, 2, 3], [1, 2, 3]) result = pa.array( spherely.to_geoarrow( arr, @@ -221,6 +222,7 @@ def test_to_geoarrow_projected(): ) np.testing.assert_allclose(coords, expected) + # Output to orthographic result = pa.array( spherely.to_geoarrow( arr, From 41ce10905de0d4210584dcb3525237be392881e3 Mon Sep 17 00:00:00 2001 From: Joris Van den Bossche Date: Fri, 15 Nov 2024 14:30:17 +0100 Subject: [PATCH 3/3] add to API docs --- docs/api.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/api.rst b/docs/api.rst index 83270db..a45da59 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -48,6 +48,7 @@ Input/Output to_wkt from_geoarrow to_geoarrow + Projection .. _api_measurement: