WGS84 <--> UTM projection header-only API

cpp/cuproj/README.md

@@ -0,0 +1,7 @@
+# <div align="left"><img src="https://rapids.ai/assets/images/rapids_logo.png" width="90px"/>&nbsp;cuProj: GPU-Accelerated Coordinate Projection</div>
+
+## Overview
+
+cuProj is a GPU-accelerated generic coordinate transformation library that transforms geospatial
+coordinates from one coordinate reference system (CRS) to another. cuProj is inspired by
+[PROJ](https://proj.org/en/9.2/), and aims to be compatible with it.

cpp/cuproj/include/cuproj/detail/pipeline.cuh

@@ -0,0 +1,120 @@
+/*
+ * Copyright (c) 2023, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <cuproj/operation/axis_swap.cuh>
+#include <cuproj/operation/clamp_angular_coordinates.cuh>
+#include <cuproj/operation/degrees_to_radians.cuh>
+#include <cuproj/operation/offset_scale_cartesian_coordinates.cuh>
+#include <cuproj/operation/operation.cuh>
+#include <cuproj/operation/transverse_mercator.cuh>
+
+#include <thrust/for_each.h>
+
+#include <iterator>
+
+namespace cuproj {
+
+namespace detail {
+
+/**
+ * @internal
+ * @brief A pipeline of projection operations applied in order to a coordinate
+ *
+ * @tparam Coordinate the coordinate type
+ * @tparam dir The direction of the pipeline, FORWARD or INVERSE
+ * @tparam T the coordinate value type
+ */
+template <typename Coordinate,
+          direction dir = direction::FORWARD,
+          typename T    = typename Coordinate::value_type>
+class pipeline {
+ public:
+  using iterator_type = std::conditional_t<dir == direction::FORWARD,
+                                           operation_type const*,
+                                           std::reverse_iterator<operation_type const*>>;
+
+  /**
+   * @brief Construct a new pipeline object with the given operations and parameters
+   *
+   * @param params The projection parameters
+   * @param ops The operations to apply
+   * @param num_stages The number of operations to apply
+   */
+  pipeline(projection_parameters<T> const& params,
+           operation_type const* ops,
+           std::size_t num_stages)
+    : params_(params), d_ops(ops), num_stages(num_stages)
+  {
+    if constexpr (dir == direction::FORWARD) {
+      first_ = d_ops;
+    } else {
+      first_ = std::reverse_iterator(d_ops + num_stages);
+    }
+  }
+
+  /**
+   * @brief Apply the pipeline to the given coordinate
+   *
+   * @param c The coordinate to transform
+   * @return The transformed coordinate
+   */
+  __device__ Coordinate operator()(Coordinate const& c) const
+  {
+    Coordinate c_out{c};
+    thrust::for_each_n(thrust::seq, first_, num_stages, [&](auto const& op) {
+      switch (op) {
+        case operation_type::AXIS_SWAP: {
+          auto op = axis_swap<Coordinate>{};
+          c_out   = op(c_out, dir);
+          break;
+        }
+        case operation_type::DEGREES_TO_RADIANS: {
+          auto op = degrees_to_radians<Coordinate>{};
+          c_out   = op(c_out, dir);
+          break;
+        }
+        case operation_type::CLAMP_ANGULAR_COORDINATES: {
+          auto op = clamp_angular_coordinates<Coordinate>{params_};
+          c_out   = op(c_out, dir);
+          break;
+        }
+        case operation_type::OFFSET_SCALE_CARTESIAN_COORDINATES: {
+          auto op = offset_scale_cartesian_coordinates<Coordinate>{params_};
+          c_out   = op(c_out, dir);
+          break;
+        }
+        case operation_type::TRANSVERSE_MERCATOR: {
+          auto op = transverse_mercator<Coordinate>{params_};
+          c_out   = op(c_out, dir);
+          break;
+        }
+      }
+    });
+    return c_out;
+  }
+
+ private:
+  projection_parameters<T> params_;
+  operation_type const* d_ops;
+  iterator_type first_;
+  std::size_t num_stages;
+};
+
+}  // namespace detail
+
+}  // namespace cuproj

cpp/cuproj/include/cuproj/ellipsoid.hpp

@@ -0,0 +1,71 @@
+/*
+ * Copyright (c) 2023, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <cassert>
+#include <cmath>
+
+namespace cuproj {
+
+/**
+ * @brief Ellipsoid parameters
+ *
+ * @tparam T Floating point type
+ */
+template <typename T>
+struct ellipsoid {
+  ellipsoid() = default;
+
+  /**
+   * @brief Construct an ellipsoid from semi-major axis and inverse flattening
+   *
+   * @param a Semi-major axis
+   * @param inverse_flattening Inverse flattening (a / (a - b), where b is the semi-minor axis)
+   */
+  constexpr ellipsoid(T a, T inverse_flattening) : a(a)
+  {
+    assert(inverse_flattening != 0.0);
+    b     = a * (1. - 1. / inverse_flattening);
+    f     = 1.0 / inverse_flattening;
+    es    = 2 * f - f * f;
+    e     = sqrt(es);
+    alpha = asin(e);
+    n     = pow(tan(alpha / 2), 2);
+  }
+
+  T a{};      // semi-major axis
+  T b{};      // semi-minor axis
+  T e{};      // first eccentricity
+  T es{};     // first eccentricity squared
+  T alpha{};  // angular eccentricity
+  T f{};      // flattening
+  T n{};      // third flattening
+};
+
+/**
+ * @brief Create the WGS84 ellipsoid
+ *
+ * @tparam T  Floating point type
+ * @return The WGS84 ellipsoid
+ */
+template <typename T>
+constexpr ellipsoid<T> make_ellipsoid_wgs84()
+{
+  return ellipsoid<T>{T{6378137.0}, T{298.257223563}};
+}
+
+}  // namespace cuproj

cpp/cuproj/include/cuproj/error.hpp

@@ -30,13 +30,13 @@ namespace cuproj {
  * @{
  */
 
-/**---------------------------------------------------------------------------*
+/**
  * @brief Exception thrown when logical precondition is violated.
  *
  * This exception should not be thrown directly and is instead thrown by the
  * CUPROJ_EXPECTS macro.
  *
- *---------------------------------------------------------------------------**/
+ */
 struct logic_error : public std::logic_error {
   logic_error(char const* const message) : std::logic_error(message) {}
   logic_error(std::string const& message) : std::logic_error(message) {}
@@ -58,7 +58,7 @@ struct cuda_error : public std::runtime_error {
 #define CUPROJ_STRINGIFY_DETAIL(x) #x
 #define CUPROJ_STRINGIFY(x)        CUPROJ_STRINGIFY_DETAIL(x)
 
-/**---------------------------------------------------------------------------*
+/**
  * @brief Macro for checking (pre-)conditions that throws an exception when
  * a condition is violated.
  *
@@ -72,13 +72,13 @@ struct cuda_error : public std::runtime_error {
  * @param[in] reason String literal description of the reason that cond is
  * expected to be true
  * @throw cuproj::logic_error if the condition evaluates to false.
- *---------------------------------------------------------------------------**/
-#define CUPROJ_EXPECTS(cond, reason)                                       \
-  (!!(cond)) ? static_cast<void>(0)                                        \
-             : throw cuspatial::logic_error("cuProj failure at: " __FILE__ \
-                                            ":" CUPROJ_STRINGIFY(__LINE__) ": " reason)
+ */
+#define CUPROJ_EXPECTS(cond, reason)                                    \
+  (!!(cond)) ? static_cast<void>(0)                                     \
+             : throw cuproj::logic_error("cuProj failure at: " __FILE__ \
+                                         ":" CUPROJ_STRINGIFY(__LINE__) ": " reason)
 
-/**---------------------------------------------------------------------------*
+/**
  * @brief Macro for checking (pre-)conditions that throws an exception when
  * a condition is violated.
  *
@@ -97,14 +97,14 @@ struct cuda_error : public std::runtime_error {
  * @throw cuproj::logic_error if the condition evaluates to false.
  * (if on device)
  * program terminates and assertion error message is printed to stderr.
- *---------------------------------------------------------------------------**/
+ */
 #ifndef __CUDA_ARCH__
 #define CUPROJ_HOST_DEVICE_EXPECTS(cond, reason) CUPROJ_EXPECTS(cond, reason)
 #else
 #define CUPROJ_HOST_DEVICE_EXPECTS(cond, reason) cuproj_assert(cond&& reason)
 #endif
 
-/**---------------------------------------------------------------------------*
+/**
  * @brief Indicates that an erroneous code path has been taken.
  *
  * In host code, throws a `cuproj::logic_error`.
@@ -116,7 +116,7 @@ struct cuda_error : public std::runtime_error {
  * ```
  *
  * @param[in] reason String literal description of the reason
- *---------------------------------------------------------------------------**/
+ */
 #define CUPROJ_FAIL(reason)                                                      \
   throw cuproj::logic_error("cuProj failure at: " __FILE__ ":" CUPROJ_STRINGIFY( \
     __LINE__) ":"                                                                \

cpp/cuproj/include/cuproj/operation/axis_swap.cuh

@@ -0,0 +1,43 @@
+/*
+ * Copyright (c) 2023, NVIDIA CORPORATION.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#pragma once
+
+#include <cuproj/operation/operation.cuh>
+
+namespace cuproj {
+
+/**
+ * @brief Axis swap operation: swap x and y coordinates
+ *
+ * @tparam Coordinate the coordinate type
+ */
+template <typename Coordinate>
+struct axis_swap : operation<Coordinate> {
+  /**
+   * @brief Swap x and y coordinates
+   *
+   * @param coord the coordinate to swap
+   * @param dir (unused) the direction of the operation
+   * @return the swapped coordinate
+   */
+  __host__ __device__ Coordinate operator()(Coordinate const& coord, direction) const
+  {
+    return Coordinate{coord.y, coord.x};
+  }
+};
+
+}  // namespace cuproj