diff --git a/.github/workflows/ci-additional.yml b/.github/workflows/ci-additional.yml index 092ab8cc..e673ee2e 100644 --- a/.github/workflows/ci-additional.yml +++ b/.github/workflows/ci-additional.yml @@ -33,6 +33,10 @@ jobs: environment-file: environment-python-dev.yml cache-environment: true + - name: Install Healpix + run: | + micromamba install healpix_cxx + - name: Build and install Fastscapelib Python run: | python -m pip install . -v --no-build-isolation diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index f0aaaa69..7a2cc3f8 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -40,6 +40,11 @@ jobs: cache-environment: true cache-downloads: false + - name: Install Healpix (Linux and MacOS) + if: matrix.cfg.os != 'windows-latest' + run: | + micromamba install healpix_cxx + - name: Configure Fastscapelib (CMake) run: | cmake -S . -B build \ @@ -78,6 +83,11 @@ jobs: create-args: >- python=${{ matrix.python-version }} + - name: Install Healpix (Linux and MacOS) + if: matrix.os != 'windows-latest' + run: | + micromamba install healpix_cxx + - name: Build and install Fastscapelib Python run: | python -m pip install . -v --no-build-isolation @@ -100,7 +110,7 @@ jobs: - name: Build and install Fastscapelib Python (Numpy 2.x) run: | - python -m pip install . --config-settings=cmake.define.FS_DOWNLOAD_XTENSOR_PYTHON=ON + python -m pip install . --config-settings=cmake.args="-DFS_DOWNLOAD_XTENSOR_PYTHON=ON;-DFS_WITH_HEALPIX=OFF" - name: Run tests (Numpy 1.xx) run: | diff --git a/CMakeLists.txt b/CMakeLists.txt index 392c5d57..0c95f70d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,9 +17,14 @@ OPTION(FS_GTEST_SRC_DIR "Build gtest from local sources" OFF) OPTION(FS_BUILD_BENCHMARKS "Build fastscapelib benchmark suite" OFF) OPTION(FS_DOWNLOAD_GBENCHMARK "Build gbenchmark from downloaded sources" OFF) OPTION(FS_GBENCHMARK_SRC_DIR "Build gbenchmark from local sources" OFF) +OPTION(FS_WITH_HEALPIX "Turn on support of Healpix grid" ON) + +if(WIN32) + set(FS_WITH_HEALPIX OFF CACHE BOOL "Healpix grid not supported on Windows" FORCE) +endif() set(CMAKE_MODULE_PATH - ${CMAKE_SOURCE_DIR}/cmake + ${CMAKE_SOURCE_DIR}/cmake/modules ${CMAKE_MODULE_PATH}) # Dependencies @@ -76,6 +81,11 @@ endif() find_package(Threads REQUIRED) +if(FS_WITH_HEALPIX) + find_package(healpix REQUIRED) + message(STATUS "Found healpix: ${healpix_INCLUDE_DIRS}/healpix_cxx") +endif() + # Installation directories # ======================== @@ -131,6 +141,10 @@ set(FASTSCAPELIB_HEADERS version.hpp ) +if(FS_WITH_HEALPIX) + list(APPEND FASTSCAPELIB_HEADERS grid/healpix_grid.hpp) +endif() + set(FASTSCAPELIB_TARGET fastscapelib) add_library(${FASTSCAPELIB_TARGET} INTERFACE) @@ -143,6 +157,10 @@ target_include_directories(${FASTSCAPELIB_TARGET} target_link_libraries(${FASTSCAPELIB_TARGET} INTERFACE xtensor Threads::Threads) +if(FS_WITH_HEALPIX) + target_link_libraries(${FASTSCAPELIB_TARGET} INTERFACE healpix) +endif() + target_compile_features(${FASTSCAPELIB_TARGET} INTERFACE cxx_std_17) # -- optional subdirectories diff --git a/benchmark/benchmark_raster_grid.cpp b/benchmark/benchmark_raster_grid.cpp index 3cc9522c..6578f797 100644 --- a/benchmark/benchmark_raster_grid.cpp +++ b/benchmark/benchmark_raster_grid.cpp @@ -170,15 +170,18 @@ namespace fastscapelib auto n = static_cast(state.range(0)); std::array shape{ { n, n } }; + const short d8_row_offsets[9] = { 0, -1, -1, 0, 1, 1, 1, 0, -1 }; + const short d8_col_offsets[9] = { 0, 0, -1, -1, -1, 0, 1, 1, 1 }; + neighbors_offsets_type offsets(grid_type::n_neighbors_max()); - auto get_neighbors_indices - = [&shape, &offsets](auto& r, auto& c) -> neighbors_offsets_type + auto get_neighbors_indices = [&shape, &offsets, &d8_row_offsets, &d8_col_offsets]( + auto& r, auto& c) -> neighbors_offsets_type { for (std::size_t k = 1; k <= grid_type::n_neighbors_max(); ++k) { - const index_t kr = r + fs::consts::d8_row_offsets[k]; - const index_t kc = c + fs::consts::d8_col_offsets[k]; + const index_t kr = r + d8_row_offsets[k]; + const index_t kc = c + d8_col_offsets[k]; offsets[k - 1] = std::array({ kr, kc }); diff --git a/cmake/modules/Findhealpix.cmake b/cmake/modules/Findhealpix.cmake new file mode 100644 index 00000000..b9890d6b --- /dev/null +++ b/cmake/modules/Findhealpix.cmake @@ -0,0 +1,77 @@ +# The MIT License (MIT) +# +# Copyright (c) 2024 Benoit Bovy +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + +# +# FindHealpix +# ----------- +# +# Find Healpix include directories and libraries. +# +# This module will set the following variables: +# +# healpix_FOUND - System has Healpix +# healpix_INCLUDE_DIRS - The Healpix include directories +# healpix_LIBRARIES - The libraries needed to use Healpix +# + +include(FindPackageHandleStandardArgs) + +find_path(healpix_INCLUDE_DIR healpix_cxx + HINTS + ENV healpix_ROOT + ENV healpix_DIR + ENV CONDA_PREFIX + ${healpix_ROOT_DIR} + PATH_SUFFIXES + include + ) + +find_library(healpix_LIBRARY + NAMES healpix_cxx + HINTS + ENV healpix_ROOT + ENV healpix_DIR + ENV CONDA_PREFIX + ${healpix_ROOT_DIR} + PATH_SUFFIXES + lib + libs + Library + ) + +find_package_handle_standard_args(healpix + REQUIRED_VARS healpix_INCLUDE_DIR healpix_LIBRARY + ) + +if(healpix_FOUND) + set(healpix_INCLUDE_DIRS ${healpix_INCLUDE_DIR}) + set(healpix_LIBRARIES ${healpix_LIBRARY}) + + add_library(healpix SHARED IMPORTED) + set_target_properties(healpix PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES ${healpix_INCLUDE_DIRS} + IMPORTED_LOCATION ${healpix_LIBRARIES} + IMPORTED_IMPLIB ${healpix_LIBRARIES} + ) + + mark_as_advanced(healpix_INCLUDE_DIRS healpix_LIBRARIES) +endif() diff --git a/doc/environment.yml b/doc/environment.yml index fb81a31a..b002ea7f 100644 --- a/doc/environment.yml +++ b/doc/environment.yml @@ -5,6 +5,7 @@ dependencies: - python=3.11 - breathe - cmake + - healpix_cxx - xtensor-python - pybind11 - scikit-build-core diff --git a/environment-dev.yml b/environment-dev.yml index e28fe785..a716272c 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -4,5 +4,3 @@ channels: dependencies: - xtensor - cmake - - gtest - - benchmark diff --git a/include/fastscapelib/grid/base.hpp b/include/fastscapelib/grid/base.hpp index 79950422..5edf7227 100644 --- a/include/fastscapelib/grid/base.hpp +++ b/include/fastscapelib/grid/base.hpp @@ -44,7 +44,8 @@ namespace fastscapelib core = 0, /**< Inner grid node */ fixed_value = 1, /**< Dirichlet boundary condition */ fixed_gradient = 2, /**< Neumann boundary condition */ - looped = 3 /**< Reflective boundaries */ + looped = 3, /**< Reflective boundaries */ + ghost = 4 /**< Inactive grid node (outisde the domain) */ }; namespace detail diff --git a/include/fastscapelib/grid/healpix_grid.hpp b/include/fastscapelib/grid/healpix_grid.hpp new file mode 100644 index 00000000..7ce82788 --- /dev/null +++ b/include/fastscapelib/grid/healpix_grid.hpp @@ -0,0 +1,380 @@ +#ifndef FASTSCAPELIB_HEALPIX_GRID_H_ +#define FASTSCAPELIB_HEALPIX_GRID_H_ + +#include "healpix_cxx/healpix_base.h" +#include "healpix_cxx/healpix_tables.h" +#include "healpix_cxx/vec3.h" + +// conflict between healpix xcomplex macro and xtl xcomplex +#undef xcomplex +#include +#include + +#include "fastscapelib/grid/base.hpp" +#include "fastscapelib/utils/consts.hpp" +#include "fastscapelib/utils/xtensor_containers.hpp" + +#include +#include +#include +#include +#include + + +namespace fastscapelib +{ + namespace detail + { + inline double vec3_distance(vec3 a, vec3 b) + { + return std::sqrt(std::pow(a.x - b.x, 2) + std::pow(a.y - b.y, 2) + + std::pow(a.z - b.z, 2)); + } + + } + + template + class healpix_grid; + + /** + * Healpix grid specialized types + */ + template + struct grid_inner_types> + { + static constexpr bool is_structured = false; + static constexpr bool is_uniform = false; + + using grid_data_type = double; + + using container_selector = S; + static constexpr std::size_t container_ndims = 1; + + static constexpr uint8_t n_neighbors_max = 8; + using neighbors_cache_type = neighbors_no_cache; + }; + + /** + * @brief 2-dimensional grid on the sphere (HEALPix). + * + * Fastscapelib grid adapter for a HEALPix (Hierarchical Equal Area + * isoLatitude Pixelation of a sphere) grid. + * + * @tparam S The container selector for data array members. + * @tparam T The integer type used to store the HEALPix grid node indices. + */ + template + class healpix_grid : public grid> + { + public: + using self_type = healpix_grid; + using base_type = grid; + using inner_types = grid_inner_types; + + using grid_data_type = typename base_type::grid_data_type; + + using container_selector = typename base_type::container_selector; + using container_type = fixed_shape_container_t; + + using neighbors_type = typename base_type::neighbors_type; + using neighbors_indices_type = typename base_type::neighbors_indices_type; + using neighbors_distances_type = typename base_type::neighbors_distances_type; + + using nodes_status_type = typename base_type::nodes_status_type; + using nodes_status_array_type = fixed_shape_container_t; + + using size_type = typename base_type::size_type; + using shape_type = typename base_type::shape_type; + + healpix_grid(T nside, + const nodes_status_array_type& nodes_status, + double radius = numeric_constants<>::EARTH_RADIUS_METERS); + + // TODO: factory calculating nside from a given approx. cell area. + + void set_nodes_status(const nodes_status_array_type& nodes_status); + + std::pair nodes_lonlat() const; + std::pair nodes_lonlat(const size_type& idx) const; + std::tuple nodes_xyz() const; + std::tuple nodes_xyz(const size_type& idx) const; + + T nside() const; + double radius() const; + + protected: + using healpix_type = T_Healpix_Base; + std::unique_ptr m_healpix_obj_ptr; + + shape_type m_shape; + size_type m_size; + double m_radius; + double m_node_area; + + nodes_status_type m_nodes_status; + + using neighbors_distances_impl_type = typename base_type::neighbors_distances_impl_type; + using neighbors_indices_impl_type = typename base_type::neighbors_indices_impl_type; + + std::vector m_neighbors_count; + std::vector m_neighbors_indices; + std::vector m_neighbors_distances; + + void set_neighbors(); + + inline container_type nodes_areas_impl() const; + inline grid_data_type nodes_areas_impl(const size_type& idx) const noexcept; + + inline size_type neighbors_count_impl(const size_type& idx) const; + + void neighbors_indices_impl(neighbors_indices_impl_type& neighbors, + const size_type& idx) const; + + inline const neighbors_distances_impl_type& neighbors_distances_impl( + const size_type& idx) const; + + static constexpr std::size_t dimension_impl() noexcept; + + friend class grid; + }; + + + /** + * @name Constructors + */ + /** + * Creates a new HEALPix grid + * + * @param nside The number of divisions along the side of a base-resolution HEALPix pixel. + * @param radius The radius of the sphere (default: Earth radius in meters). + */ + template + healpix_grid::healpix_grid(T nside, + const nodes_status_array_type& nodes_status, + double radius) + : base_type(0) + , m_radius(radius) + { + m_healpix_obj_ptr + = std::make_unique(nside, Healpix_Ordering_Scheme::RING, SET_NSIDE); + + m_size = static_cast(m_healpix_obj_ptr->Npix()); + m_shape = { static_cast(m_size) }; + m_node_area = 4.0 * M_PI * m_radius * m_radius / m_size; + + // will also compute grid node neighbors + set_nodes_status(nodes_status); + } + //@} + + template + void healpix_grid::set_nodes_status(const nodes_status_array_type& nodes_status) + { + if (!xt::same_shape(nodes_status.shape(), m_shape)) + { + throw std::invalid_argument( + "invalid shape for nodes_status array (expects shape [N] where N is the total number of nodes)"); + } + m_nodes_status = nodes_status; + + // maybe invalidates the grid node neighbors so it must be (re)computed + set_neighbors(); + } + + /** + * @name HEALPix specific methods + */ + /** + * Returns the longitude and latitude of a given grid node (HEALPix cell centroid), in radians. + * + * @param idx The grid node indice. + */ + template + std::pair healpix_grid::nodes_lonlat(const size_type& idx) const + { + auto ang = m_healpix_obj_ptr->pix2ang(static_cast(idx)); + + return std::make_pair(double(ang.phi), + xt::numeric_constants::PI_2 - ang.theta); + } + + /** + * Returns the longitude and latitude of all grid nodes (HEALPix cell centroids), in radians. + */ + template + auto healpix_grid::nodes_lonlat() const -> std::pair + { + container_type lon(m_shape); + container_type lat(m_shape); + + for (size_type idx = 0; idx < m_size; ++idx) + { + auto ang = m_healpix_obj_ptr->pix2ang(static_cast(idx)); + lon(idx) = ang.phi; + lat(idx) = xt::numeric_constants::PI_2 - ang.theta; + } + + return std::make_pair(std::move(lon), std::move(lat)); + } + + /** + * Returns the x, y, z cartesian coordinates of a given grid node (HEALPix cell centroid). + * + * @param idx The grid node indice. + */ + template + std::tuple healpix_grid::nodes_xyz(const size_type& idx) const + { + auto vec = m_healpix_obj_ptr->pix2vec(static_cast(idx)); + + return std::make_tuple( + vec.x * m_radius, vec.y * m_radius, vec.z * m_radius); + } + + /** + * Returns the x, y, z cartesian coordinates of all grid nodes (HEALPix cell centroids). + */ + template + auto healpix_grid::nodes_xyz() const + -> std::tuple + { + container_type x(m_shape); + container_type y(m_shape); + container_type z(m_shape); + + for (size_type idx = 0; idx < m_size; ++idx) + { + auto vec = m_healpix_obj_ptr->pix2vec(static_cast(idx)); + x(idx) = vec.x * m_radius; + y(idx) = vec.y * m_radius; + z(idx) = vec.z * m_radius; + } + + return std::make_tuple( + std::move(x), std::move(y), std::move(z)); + } + //@} + + template + void healpix_grid::set_neighbors() + { + m_neighbors_count.resize(m_size); + m_neighbors_indices.resize(m_size); + m_neighbors_distances.resize(m_size); + + fix_arr temp_neighbors_indices; + + for (size_type inode = 0; inode < m_size; inode++) + { + size_type neighbors_count = 0; + + if (m_nodes_status[inode] == node_status::looped) + { + throw std::invalid_argument("node_status::looped is not allowed in " + "healpix grid"); + } + else if (m_nodes_status[inode] == node_status::ghost) + { + m_neighbors_count[inode] = neighbors_count; + continue; + } + + T inode_ = static_cast(inode); + auto inode_vec3 = m_healpix_obj_ptr->pix2vec(inode_); + + m_healpix_obj_ptr->neighbors(inode_, temp_neighbors_indices); + + for (size_type k = 1; k < 8; ++k) + { + auto ineighbor = temp_neighbors_indices[k]; + auto ineighbor_ = static_cast(ineighbor); + + if (ineighbor > -1 && m_nodes_status[ineighbor_] != node_status::ghost) + { + m_neighbors_indices[inode][neighbors_count] = ineighbor_; + + auto ineighbor_vec3 = m_healpix_obj_ptr->pix2vec(ineighbor); + m_neighbors_distances[inode][neighbors_count] + = detail::vec3_distance(inode_vec3, ineighbor_vec3); + + neighbors_count++; + } + } + + m_neighbors_count[inode] = neighbors_count; + } + } + + /** + * @name HEALPix specific properties + */ + /** + * HEALPix's Nside (number of divisions along the side of a HEALPix + * base-resolution pixel). + * + * The higher the value is, the finer the resolution of the grid is. + */ + template + auto healpix_grid::nside() const -> T + { + return m_healpix_obj_ptr->Nside(); + } + + /** + * Sphere radius + */ + template + double healpix_grid::radius() const + { + return m_radius; + } + //@} + + template + inline auto healpix_grid::nodes_areas_impl() const -> container_type + { + return xt::broadcast(m_node_area, m_shape); + } + + template + inline auto healpix_grid::nodes_areas_impl(const size_type& /*idx*/) const noexcept + -> grid_data_type + { + return m_node_area; + } + + template + inline auto healpix_grid::neighbors_count_impl(const size_type& idx) const -> size_type + { + return m_neighbors_count[idx]; + } + + template + void healpix_grid::neighbors_indices_impl(neighbors_indices_impl_type& neighbors, + const size_type& idx) const + { + const auto& size = m_neighbors_count[idx]; + + for (size_type i = 0; i < size; i++) + { + neighbors[i] = m_neighbors_indices[idx][i]; + } + } + + template + auto healpix_grid::neighbors_distances_impl(const size_type& idx) const + -> const neighbors_distances_impl_type& + { + return m_neighbors_distances[idx]; + } + + template + constexpr std::size_t healpix_grid::dimension_impl() noexcept + { + return 2; + } +} + +#endif // FASTSCAPELIB_HEALPIX_GRID_H_ diff --git a/include/fastscapelib/utils/consts.hpp b/include/fastscapelib/utils/consts.hpp index 0d0e0e9b..3ee9eaf4 100644 --- a/include/fastscapelib/utils/consts.hpp +++ b/include/fastscapelib/utils/consts.hpp @@ -8,26 +8,10 @@ namespace fastscapelib { - namespace consts + template + struct numeric_constants { - - /* - * @brief Row index offsets (from a central grid point) of D8 - * neighbour directions. - * - * The first element correspond to the central grid point itself. - */ - const short d8_row_offsets[9] = { 0, -1, -1, 0, 1, 1, 1, 0, -1 }; - - - /** - * @brief Column index offsets (from a central grid point) of D8 - * neighbour directions. - * - * The first element correspond to the central grid point itself. - */ - const short d8_col_offsets[9] = { 0, 0, -1, -1, -1, 0, 1, 1, 1 }; - - } // namespace consts + static constexpr T EARTH_RADIUS_METERS = 6.371e6; + }; } // namespace fastscapelib diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index d14f3672..3f1e3011 100644 --- a/python/CMakeLists.txt +++ b/python/CMakeLists.txt @@ -44,6 +44,10 @@ target_compile_definitions(${PY_FASTSCAPELIB} VERSION_INFO=${PROJECT_VERSION} ) +if(FS_WITH_HEALPIX) + target_compile_definitions(${PY_FASTSCAPELIB} PRIVATE WITH_HEALPIX) +endif() + # Installation # ============ diff --git a/python/fastscapelib/__init__.py b/python/fastscapelib/__init__.py index e4fd47e1..4af7ef0c 100644 --- a/python/fastscapelib/__init__.py +++ b/python/fastscapelib/__init__.py @@ -16,7 +16,8 @@ SingleFlowRouter, create_flow_kernel, ) -from fastscapelib.grid import ( # type: ignore[import] +from fastscapelib.grid import ( # type: ignore[import, attr-defined] + HealpixGrid, Neighbor, Node, NodeStatus, @@ -32,6 +33,7 @@ __all__ = [ "__version__", "DiffusionADIEroder", + "HealpixGrid", "SPLEroder", "Neighbor", "Node", diff --git a/python/fastscapelib/grid/__init__.py b/python/fastscapelib/grid/__init__.py index 8c421870..b306776e 100644 --- a/python/fastscapelib/grid/__init__.py +++ b/python/fastscapelib/grid/__init__.py @@ -1,4 +1,4 @@ -from _fastscapelib_py.grid import ( # type: ignore[import] +from _fastscapelib_py.grid import ( Neighbor, Node, NodeStatus, @@ -11,7 +11,14 @@ TriMesh, ) +try: + from _fastscapelib_py.grid import HealpixGrid +except ImportError: + # Python module built with no healpix grid support (e.g., Windows) + HealpixGrid = None + __all__ = [ + "HealpixGrid", "Neighbor", "Node", "NodeStatus", diff --git a/python/src/flow_graph.cpp b/python/src/flow_graph.cpp index 36e58cfb..98ecbf78 100644 --- a/python/src/flow_graph.cpp +++ b/python/src/flow_graph.cpp @@ -421,6 +421,9 @@ add_flow_graph_bindings(py::module& m) fs::register_py_flow_graph_init(pyfgraph, true); fs::register_py_flow_graph_init(pyfgraph, false); fs::register_py_flow_graph_init(pyfgraph, false); +#ifdef WITH_HEALPIX + fs::register_py_flow_graph_init(pyfgraph, false); +#endif pyfgraph.def_property_readonly( "operators", diff --git a/python/src/grid.cpp b/python/src/grid.cpp index db67a5e8..bab1c2a1 100644 --- a/python/src/grid.cpp +++ b/python/src/grid.cpp @@ -4,6 +4,7 @@ #include "fastscapelib/grid/base.hpp" #include "fastscapelib/grid/profile_grid.hpp" #include "fastscapelib/grid/raster_grid.hpp" +#include "fastscapelib/utils/consts.hpp" #include "xtensor-python/pytensor.hpp" #include "xtensor-python/pyarray.hpp" @@ -134,7 +135,7 @@ add_grid_bindings(py::module& m) Create a new mesh with a dictionary of node status (optional). - 2. ``__init__(self, points: numpy.ndarray, triangles: numpy.ndarray, nodes_status: np.ndarray) -> None`` + 2. ``__init__(self, points: numpy.ndarray, triangles: numpy.ndarray, nodes_status: numpy.ndarray) -> None`` Create a new mesh with an array of node status. @@ -165,6 +166,106 @@ add_grid_bindings(py::module& m) fs::register_base_grid_properties(tmesh); fs::register_grid_methods(tmesh); +#ifdef WITH_HEALPIX + /* + ** Healpix grid + */ + + py::class_ hgrid(m, "HealpixGrid", "A HEALPix grid."); + + hgrid.def( + py::init(), + py::arg("nside"), + py::arg("nodes_status"), + py::arg("radius") = fs::numeric_constants<>::EARTH_RADIUS_METERS, + R"doc(__init__(nside: int, nodes_status: numpy.ndarray, radius: float = 6.371e6) -> None + + HEALPix grid initializer. + + Parameters + ---------- + nside : int + Number of divisions along the side of a HEALPix base-resolution pixel. + A higher value sets a grid with a finer resolution and a greater size N + (i.e., total number of grid nodes). + nodes_status : array-like + Array of shape [N] setting the status of all grid nodes at once. + radius : float, optional + Sphere radius (default, Earth radius in meters). + + )doc"); + + fs::register_grid_static_properties(hgrid); + fs::register_base_grid_properties(hgrid); + fs::register_grid_methods(hgrid); + + hgrid + .def_property_readonly( + "nside", + &fs::py_healpix_grid::nside, + "HEALPix's Nside (number of division along the side of a HEALPix base-resolution pixel).") + .def_property_readonly("radius", &fs::py_healpix_grid::radius, "Sphere radius"); + + hgrid + .def("nodes_lonlat", + py::overload_cast( + &fs::py_healpix_grid::nodes_lonlat, py::const_), + py::arg("idx"), + R"doc(nodes_lonlat(*args) -> tuple + + Return the longitude and latitude coordinates of one or all grid nodes + (HEALPix cell centroids), in radians. + + Overloaded method that supports the following signatures: + + 1. ``nodes_lonlat(idx: int) -> tuple[double, double]`` + + 2. ``nodes_lonlat() -> tuple[numpy.ndarray, numpy.ndarray]`` + + Parameters + ---------- + idx : int + Grid node indice. + + Returns + ------- + lonlat : tuple + Longitude and latitude coordinates (scalars or arrays) in radians. + + )doc") + .def("nodes_lonlat", py::overload_cast<>(&fs::py_healpix_grid::nodes_lonlat, py::const_)); + + hgrid + .def("nodes_xyz", + py::overload_cast( + &fs::py_healpix_grid::nodes_xyz, py::const_), + py::arg("idx"), + R"doc(nodes_xyz(*args) -> tuple + + Return the x,y,z cartesian coordinates of one or all grid nodes + (HEALPix cell centroids). + + Overloaded method that supports the following signatures: + + 1. ``nodes_xyz(idx: int) -> tuple[double, double, double]`` + + 2. ``nodes_xyz() -> tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]`` + + Parameters + ---------- + idx : int + Grid node indice. + + Returns + ------- + x, y, z : floats (scalar or arrays) + X, Y and Z coordinates in the same scale and units than the grid's + sphere radius. + + )doc") + .def("nodes_xyz", py::overload_cast<>(&fs::py_healpix_grid::nodes_xyz, py::const_)); +#endif + /* ** Profile Grid */ diff --git a/python/src/grid.hpp b/python/src/grid.hpp index 2b43b759..4cce985c 100644 --- a/python/src/grid.hpp +++ b/python/src/grid.hpp @@ -12,6 +12,9 @@ #include "xtensor-python/pytensor.hpp" #include "fastscapelib/grid/base.hpp" +#ifdef WITH_HEALPIX +#include "fastscapelib/grid/healpix_grid.hpp" +#endif #include "fastscapelib/grid/profile_grid.hpp" #include "fastscapelib/grid/raster_grid.hpp" #include "fastscapelib/grid/trimesh.hpp" @@ -29,6 +32,9 @@ namespace fastscapelib using py_profile_grid = fs::profile_grid; using py_raster_grid = fs::raster_grid; using py_trimesh = fs::trimesh_xt; +#ifdef WITH_HEALPIX + using py_healpix_grid = fs::healpix_grid; +#endif template void register_grid_static_properties(py::class_& pyg) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0a9eb3e0..454fb7f5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -86,6 +86,10 @@ set(FASTSCAPELIB_TEST_SRC test_trimesh.cpp ) +if(FS_WITH_HEALPIX) + list(APPEND FASTSCAPELIB_TEST_SRC test_healpix_grid.cpp) +endif() + enable_testing() # -- build a target for each test diff --git a/test/test_healpix_grid.cpp b/test/test_healpix_grid.cpp new file mode 100644 index 00000000..f0ff168d --- /dev/null +++ b/test/test_healpix_grid.cpp @@ -0,0 +1,146 @@ +#include "fastscapelib/grid/base.hpp" +#include "fastscapelib/grid/healpix_grid.hpp" +#include "fastscapelib/utils/consts.hpp" + +#include "xtensor/xtensor.hpp" +#include "xtensor/xmath.hpp" + +#include "gtest/gtest.h" + + +namespace fs = fastscapelib; + + +namespace fastscapelib +{ + namespace testing + { + + class healpix_grid : public ::testing::Test + { + protected: + int nside = 32; + std::size_t size = 12288; + + xt::xtensor nodes_status = xt::zeros({ size }); + }; + + TEST_F(healpix_grid, static_expr) + { + EXPECT_EQ(fs::healpix_grid<>::is_structured(), false); + EXPECT_EQ(fs::healpix_grid<>::is_uniform(), false); + EXPECT_EQ(fs::healpix_grid<>::n_neighbors_max(), 8u); + EXPECT_EQ(fs::healpix_grid<>::container_ndims(), 1); + } + + TEST_F(healpix_grid, ctor) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + EXPECT_EQ(grid.nside(), nside); + EXPECT_EQ(grid.size(), size); + EXPECT_EQ(grid.radius(), fs::numeric_constants<>::EARTH_RADIUS_METERS); + } + + TEST_F(healpix_grid, nodes_status) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + EXPECT_EQ(grid.nodes_status(), nodes_status); + + { + SCOPED_TRACE("setting ghost node and updating node status should update neighbors"); + + auto nodes_status2 = nodes_status; + nodes_status2(1) = fs::node_status::ghost; + grid.set_nodes_status(nodes_status2); + EXPECT_EQ(grid.neighbors_count(0), 6); + + auto nodes_status3 = nodes_status; + nodes_status3(0) = fs::node_status::ghost; + grid.set_nodes_status(nodes_status3); + EXPECT_EQ(grid.neighbors_count(0), 0); + } + + { + SCOPED_TRACE("looped boundary conditions not supported"); + + auto nodes_status2 = nodes_status; + nodes_status2(0) = fs::node_status::looped; + EXPECT_THROW(grid.set_nodes_status(nodes_status2), std::invalid_argument); + } + } + + TEST_F(healpix_grid, nodes_areas) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + double expected = 41509152987.45021; + EXPECT_NEAR(grid.nodes_areas(0), expected, 1e-5); + + xt::xtensor expected_arr({ grid.size() }, expected); + EXPECT_TRUE(xt::allclose(grid.nodes_areas(), expected_arr)); + } + + TEST_F(healpix_grid, nodes_lonlat) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + auto actual = grid.nodes_lonlat(0); + auto actual_arr = grid.nodes_lonlat(); + auto actual_arr0 + = std::make_pair(actual_arr.first(0), actual_arr.second(0)); + + auto expected = std::make_pair(0.7853981633974483, 1.5452801164374776); + + EXPECT_EQ(actual, actual_arr0); + EXPECT_EQ(actual, expected); + } + + TEST_F(healpix_grid, nodes_xyz) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + auto actual = grid.nodes_xyz(0); + auto actual_arr = grid.nodes_xyz(); + auto actual_arr0 = std::make_tuple( + std::get<0>(actual_arr)(0), std::get<1>(actual_arr)(0), std::get<2>(actual_arr)(0)); + + auto expected = std::make_tuple( + 114937.4753788243, 114937.4753788243, 6368926.106770833); + + EXPECT_EQ(actual, actual_arr0); + EXPECT_EQ(actual, expected); + } + + TEST_F(healpix_grid, neighbors_count) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + EXPECT_EQ(grid.neighbors_count(0), 7); + EXPECT_EQ(grid.neighbors_count(1984), 6); + } + + TEST_F(healpix_grid, neighbors_indices) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + EXPECT_EQ(grid.neighbors_indices(0), + (xt::xtensor{ 11, 3, 2, 1, 6, 5, 13 })); + } + + TEST_F(healpix_grid, nodes_distances) + { + auto grid = fs::healpix_grid<>(nside, nodes_status); + + xt::xtensor expected({ 0.04752047, + 0.03608146, + 0.05102688, + 0.03608146, + 0.04752047, + 0.02914453, + 0.0510435 }); + EXPECT_TRUE(xt::allclose(grid.neighbors_distances(0), expected)); + } + } +}