feat: add projections

CMakeLists.txt

@@ -201,6 +201,7 @@ add_library(s2geography
     src/s2geography/accessors.cc
     src/s2geography/distance.cc
     src/s2geography/predicates.cc
+    src/s2geography/projections.cc
     src/s2geography/build.cc
     src/s2geography/geography.cc
     src/s2geography/geoarrow.cc

src/s2geography.h

@@ -15,5 +15,6 @@
 #include "s2geography/index.h"
 #include "s2geography/linear-referencing.h"
 #include "s2geography/predicates.h"
+#include "s2geography/projections.h"
 #include "s2geography/wkt-reader.h"
 #include "s2geography/wkt-writer.h"

src/s2geography/geoarrow.cc

@@ -6,7 +6,6 @@
 #include "geoarrow/geoarrow.h"
 #include "s2/s1angle.h"
 #include "s2/s2edge_tessellator.h"
-#include "s2/s2projections.h"
 #include "s2geography/geography.h"
 
 namespace s2geography {
@@ -15,10 +14,6 @@ namespace geoarrow {
 
 const char* version() { return GeoArrowVersion(); }
 
-S2::Projection* lnglat() {
-  static S2::PlateCarreeProjection projection(180);
-  return &projection;
-}
 
 // This should really be in nanoarrow or geoarrow
 // https://github.com/geoarrow/geoarrow-c-geos/blob/33ad0ba21c76c09e9d72fc4e4ae0b9ff9da61848/src/geoarrow_geos/geoarrow_geos.c#L323-L360

src/s2geography/geoarrow.h

@@ -8,6 +8,7 @@
 #include "s2/s2projections.h"
 #include "s2geography/arrow_abi.h"
 #include "s2geography/geography.h"
+#include "s2geography/projections.h"
 
 namespace s2geography {
 
@@ -16,24 +17,22 @@ namespace geoarrow {
 /// \brief Inspect the underlying GeoArrow implementation version
 const char* version();
 
-S2::Projection* lnglat();
-
 /// \brief Options used to build Geography objects from GeoArrow arrays
 
 class TessellationOptions {
  public:
   TessellationOptions()
-      : projection_(lnglat()),
+      : projection_(std::move(lnglat())),
         tessellate_tolerance_(S1Angle::Infinity()) {}
-  S2::Projection* projection() const { return projection_; }
-  void set_projection(S2::Projection* projection) { projection_ = projection; }
+  S2::Projection* projection() const { return projection_.get(); }
+  void set_projection(std::shared_ptr<S2::Projection> projection) { projection_ = std::move(projection); }
   S1Angle tessellate_tolerance() const { return tessellate_tolerance_; }
   void set_tessellate_tolerance(S1Angle tessellate_tolerance) {
     tessellate_tolerance_ = tessellate_tolerance;
   }
 
  protected:
-  S2::Projection* projection_;
+  std::shared_ptr<S2::Projection> projection_;
   S1Angle tessellate_tolerance_;
 };
 

src/s2geography/projections.cc

@@ -0,0 +1,70 @@
+
+#include "s2/s2projections.h"
+
+#include <s2/r2.h>
+#include <s2/s2latlng.h>
+#include <s2/s2point.h>
+#include <s2/s2pointutil.h>
+
+namespace s2geography {
+
+std::shared_ptr<S2::Projection> lnglat() {
+  std::shared_ptr<S2::Projection> projection = std::make_shared<S2::PlateCarreeProjection>(180);
+  return std::move(projection);
+}
+
+std::shared_ptr<S2::Projection> pseudo_mercator() {
+  // the semi-major axis of the WGS 84 ellipsoid is 6378137 meters
+  // -> half of the circumference of the sphere is PI * 6378137 = 20037508.3427892
+  std::shared_ptr<S2::Projection> projection = std::make_shared<S2::MercatorProjection>(20037508.3427892);
+  return std::move(projection);
+}
+
+class OrthographicProjection: public S2::Projection {
+public:
+  OrthographicProjection(const S2LatLng& centre):
+      centre_(centre) {
+    z_axis_ = S2Point(0, 0, 1);
+    y_axis_ = S2Point(0, 1, 0);
+  }
+
+  // Converts a point on the sphere to a projected 2D point.
+  R2Point Project(const S2Point& p) const {
+    S2Point out = S2::Rotate(p, z_axis_, -centre_.lng());
+    out = S2::Rotate(out, y_axis_, centre_.lat());
+    return R2Point(out.y(), out.z());
+  }
+
+  // Converts a projected 2D point to a point on the sphere.
+  S2Point Unproject(const R2Point& p) const {
+    double y = p.x();
+    double z = p.y();
+    double x = sqrt(1.0 - y * y - z * z);
+    S2Point pp(x, y, z);
+    S2Point out = S2::Rotate(pp, y_axis_, -centre_.lat());
+    out = S2::Rotate(out, z_axis_, centre_.lng());
+    return out;
+  }
+
+  R2Point FromLatLng(const S2LatLng& ll) const {
+    return Project(ll.ToPoint());
+  }
+
+  S2LatLng ToLatLng(const R2Point& p) const {
+    return S2LatLng(Unproject(p));
+  }
+
+  R2Point wrap_distance() const {return R2Point(0, 0); }
+
+private:
+  S2LatLng centre_;
+  S2Point z_axis_;
+  S2Point y_axis_;
+};
+
+std::shared_ptr<S2::Projection> orthographic(const S2LatLng& centre) {
+  std::shared_ptr<S2::Projection> projection = std::make_shared<OrthographicProjection>(centre);
+  return std::move(projection);
+}
+
+}  // namespace s2geography

src/s2geography/projections.h

@@ -0,0 +1,23 @@
+
+#pragma once
+
+#include "s2/s2projections.h"
+
+namespace s2geography {
+
+// Constructs the "plate carree" projection which maps coordinates on the sphere
+// to (longitude, latitude) pairs.
+// The x coordinates (longitude) span [-180, 180] and the y coordinates (latitude)
+// span [-90, 90].
+std::shared_ptr<S2::Projection> lnglat();
+
+// Constructs the spherical Mercator projection. When used together with WGS84
+// coordinates, known as the "Web Mercator" or "WGS84/Pseudo-Mercator" projection.
+std::shared_ptr<S2::Projection> pseudo_mercator();
+
+// Constructs an orthographic projection with the given centre point. The
+// resulting coordinates depict a single hemisphere of the globe as it appears
+// from outer space, centred on the given point.
+std::shared_ptr<S2::Projection> orthographic(const S2LatLng& centre);
+
+}  // namespace s2geography