Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for getting projected coordinates #23

Draft
wants to merge 1 commit into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@ endif()
add_library(spherely MODULE
src/geography.cpp
src/predicates.cpp
src/projections.cpp
src/spherely.cpp
)

Expand Down
2 changes: 1 addition & 1 deletion src/geography.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -244,7 +244,7 @@ void init_geography(py::module &m) {
)pbdoc");

pygeography.def("__repr__", [](const Geography &geog) {
s2geog::WKTWriter writer;
s2geog::WKTWriter writer(6);
return writer.write_feature(geog.geog());
});

Expand Down
49 changes: 49 additions & 0 deletions src/projections.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
#include <stdexcept>
#include <iostream>

#include <s2/s2projections.h>
#include <s2/s2latlng.h>
#include <s2geography.h>

#include "geography.hpp"
#include "pybind11.hpp"

namespace py = pybind11;
namespace s2geog = s2geography;
using namespace spherely;


std::tuple<double, double> project_mercator(py::object obj) {
auto geog = (static_cast<PyObjectGeography&>(obj)).as_geog_ptr();
if (geog->geog_type() != GeographyType::Point) {
throw py::value_error("test");
}
auto point = static_cast<Point*>(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("intersects", py::vectorize(&intersects), py::arg("a"), py::arg("b"),
// R"pbdoc(
// Returns True if A and B share any portion of space.

// Intersects implies that overlaps, touches and within are True.

// Parameters
// ----------
// a, b : :py:class:`Geography` or array_like
// Geography object(s)

// )pbdoc");

m.def("project_mercator", &project_mercator);
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A problem with vectorizing a function like this is that it returns multiple values for each input value (so in numpy ufunc terms, it would be a generalize ufunc with a signature of "()->(n)" (adding one dimension)). But this is not something that py::vectorize supports AFAIK. So we might need to write this more manually (eg writing the loop ourselves, only supporting 1D input, which is what we do in shapely as well).

And for non-point geographies, it doesn't fit in the ufunc paradigm anyway, since you have a variable number of coordinate pairs that gets returned (like shapely.get_coordinates)

Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we follow shapely.get_coordinates API then? I.e., return a (N, 2) array with a flat index. I guess with tessellation enabled we'll almost always want to get the index.


}
3 changes: 3 additions & 0 deletions src/spherely.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,8 @@ namespace py = pybind11;

void init_geography(py::module&);
void init_predicates(py::module&);
void init_projections(py::module&);


PYBIND11_MODULE(spherely, m) {
m.doc() = R"pbdoc(
Expand All @@ -19,6 +21,7 @@ PYBIND11_MODULE(spherely, m) {

init_geography(m);
init_predicates(m);
init_projections(m);

#ifdef VERSION_INFO
m.attr("__version__") = MACRO_STRINGIFY(VERSION_INFO);
Expand Down