Merge pull request #35 from Quetzal-framework/develop

doc: update examples and include system of units in geography

.github/workflows/cmake.yml

@@ -18,7 +18,7 @@ jobs:
     runs-on: ubuntu-22.04
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
       with:
         fetch-depth: 0
         submodules: recursive

.gitmodules

@@ -1,3 +1,4 @@
 [submodule "docs/doxygen-awesome-css"]
 	path = docs/doxygen-awesome-css
 	url = https://github.com/jothepro/doxygen-awesome-css.git
+	branch = main

docs/4-tutorials.md

@@ -31,7 +31,7 @@
   - @subpage spatial_graphs
     - @subpage spatial_graph_construction
     - @subpage dispersal_kernels
-    - @subpage Setting Vertex/Edge information at construction time
+    - @subpage spatial_graphs_information
     - @subpage Growth Expressions
     - @subpage Memory Management
 
@@ -958,6 +958,26 @@ dynamically adding vertices and edges based on user-defined conditions, resultin
 
 ---
 
+[//]: # (----------------------------------------------------------------------)
+@page spatial_graphs_information Embedding information along Vertices and Edges
+@tableofcontents
+
+The focal point of the provided code snippet revolves around the process of constructing a spatial graph from raster data. This spatial graph serves as a powerful tool for visualizing and analyzing spatial relationships and patterns present within the underlying geographic data.
+
+In addition to forming the basic structure of vertices and edges, the code snippet emphasizes the integration of supplementary information. By utilizing custom structures, developers can enrich vertices and edges with additional details specific to their application domain. For instance, vertices can be decorated with labels or population data, while edges can incorporate distance metrics or other spatial attributes.
+
+This approach enables the creation of spatial graphs that not only capture the topological relationships between locations but also incorporate contextual information essential for comprehensive spatial analysis and modeling. 
+
+**Input**
+
+@include{lineno} geography_graph_7.cpp
+
+**Output**
+
+@include{lineno} geography_graph_7.txt
+
+---
+
 [//]: # (----------------------------------------------------------------------)
 @page dispersal_kernels Distance-based dispersal Kernels
 @tableofcontents
@@ -984,16 +1004,3 @@ Quetzal incorporates various kernel types available in the `quetzal::demography:
 
 ---
 
-[//]: # (----------------------------------------------------------------------)
-@page spatial_graph_information Embedding Vertex and/or Edge information
-@tableofcontents
-
-**Input**
-
-@include{lineno} geography_graph_7.cpp
-
-**Output**
-
-@include{lineno} geography_graph_7.txt
-
----

example/demography_dispersal_kernel.cpp

@@ -14,7 +14,7 @@ int main()
     // Dispersal kernels operate on distances between geographic coordinates
     constexpr auto source = lonlat(2.25, 48.9176);
     constexpr auto target = lonlat(2.25, 48.95);
-    constexpr auto distance = source.great_circle_distance_to(target) * km;
+    constexpr auto distance = source.great_circle_distance_to(target);
 
     // Used to parametrize the kernel
     constexpr auto param = gaussian<>::param_type{.a = 1000. * m};

example/geography_graph_7.cpp

@@ -2,17 +2,21 @@
 
 namespace geo = quetzal::geography;
 
+// User-defined data structure to store arbitrary vertex information
 struct MyVertexInfo {
-    std::string label = "NA";
-    std::vector<double> pop_data_chunk;
+    std::string population_label = "NA";
+    std::vector<double> population_data_chunk;
+    // ... anything else required by the user context
 };
 
+// User-defined data structure to store arbitrary edge information
 struct MyEdgeInfo {
-    double distance;
+    mp_units::quantity<mp_units::si::metre> distance;
+    // needs to be default-constructible
     MyEdgeInfo() = default;
-    // required by from_grid
-    MyEdgeInfo(auto u, auto v, auto const& land){ 
-        distance = land.to_lonlat(u).great_circle_distance_to( land.to_lonlat(v) );
+    // constructor required by from_grid method: takes a source and a target vertex descriptor, and a spatial grid
+    MyEdgeInfo(auto s, auto t, auto const& raster){ 
+        distance = raster.to_lonlat(s).great_circle_distance_to( raster.to_lonlat(t) );
     }
 };
 
@@ -24,8 +28,11 @@ int main()
     std::vector<int> times(10);
     std::iota(times.begin(), times.end(), 2001);
 
+    // Landscape construction
     using landscape_type = geo::raster<>;
     auto land = geo::raster<>::from_file(file, times);
+
+    // Graph construction with vertex and edge info
     auto graph = geo::from_grid(land, MyVertexInfo(), MyEdgeInfo(), geo::connect_fully(), geo::isotropy(), geo::mirror());
 
     std::cout << "Graph has " << graph.num_vertices() << " vertices, " << graph.num_edges() << " edges." << std::endl;

example/geography_graph_complete_1.cpp

@@ -1,5 +1,4 @@
 #include "quetzal/quetzal.hpp"
-#include <mp-units/ostream.h>
 
 #include <cassert>
 #include <filesystem>
@@ -41,7 +40,7 @@ int main()
     // Transform the edges of the graph to a vector of probability to disperse from source to target
     auto probabilities = graph.edges()
     	| std::views::transform( [&](auto const& e){ return sphere_distance( graph.source( e ), graph.target( e ) );} )
-    	| std::views::transform( [&](auto distance){ return kernel::pdf( distance*km, { .a=200.*km , .b=5.5 } ) ;} );
+    	| std::views::transform( [&](auto distance){ return kernel::pdf( distance, { .a=200.*km , .b=5.5 } ) ;} );
 
     // Print the result
     for (auto const& i : probabilities)

src/include/quetzal/geography/coordinates.hpp

@@ -12,10 +12,17 @@
 #include <cmath> // trigonometry
 #include <ostream>
 #include <tuple>
+#include <numbers>
+
+#include <mp-units/systems/isq/isq.h> // QuantityOf<isq::length>
+#include <mp-units/systems/si/si.h>
+#include <mp-units/ostream.h>
 
 namespace quetzal::geography
 {
 
+using std::numbers::pi;
+
 class latlon;
 class lonlat;
 class rowcol;
@@ -44,7 +51,7 @@ template <GridCoordinate C1, GridCoordinate C2> auto operator<=>(const C1 &c1, c
 /// @brief GridCoordinates are streamable
 template <GridCoordinate C> std::ostream &operator<<(std::ostream &stream, const C &c)
 {
-    stream << "(Column: " << c.col << " , Row: " << c.row << ")";
+    stream << "(Column: " << c.col << ", Row: " << c.row << ")";
     return stream;
 }
 
@@ -113,7 +120,7 @@ template <Coordinate C1, Coordinate C2> auto operator<=>(const C1 &c1, const C2
 /// @brief Coordinates are streamable
 template <Coordinate C> std::ostream &operator<<(std::ostream &stream, const C &c)
 {
-    stream << "(Lat: " << c.lat << " , Lon: " << c.lon << ")";
+    stream << "(Lat: " << c.lat << ", Lon: " << c.lon << ")";
     return stream;
 }
 
@@ -130,15 +137,13 @@ void check(double lat, double lon)
 /// @brief Converts decimal degrees to radians
 constexpr double deg2rad(double deg) noexcept
 {
-    const double pi = 3.14159265358979323846;
-    return (deg * pi / 180);
+    return (deg * pi / 180.0);
 }
 
 /// @brief Converts radians to decimal degrees
 constexpr double rad2deg(double rad) noexcept
 {
-    const double pi = 3.14159265358979323846;
-    return (rad * 180 / pi);
+    return (rad * 180.0 / pi);
 }
 
 /// @brief Returns the distance between two points on the Earth in kilometers
@@ -176,9 +181,6 @@ struct latlon
     /// @typedef Latitude and Longitude value type
     using decimal_degree = double;
 
-    /// @typedef The great circle distance value type
-    using km = double;
-
     /// @brief Latitude
     decimal_degree lat;
 
@@ -199,14 +201,15 @@ struct latlon
     /// @param other the other coordinate
     /// @return the distance in kilometers
     /// @invariant The returned distance is not negative.
-    template <class T> constexpr km great_circle_distance_to(const T &other) const noexcept
+    template <class T> constexpr 
+    mp_units::quantity<mp_units::si::metre> great_circle_distance_to(const T &other) const noexcept
     {
-        const km d = details::distanceEarth(this->lat, this->lon, other.lat, other.lon);
+        const double d = details::distanceEarth(this->lat, this->lon, other.lat, other.lon);
         assert(d >= 0.);
-        return d;
+        return d * 1000.0 * mp_units::si::metre;
     }
+};
 
-}; // end struct latlon
 static_assert(Coordinate<latlon>);
 
 /// @brief Geographic coordinates.
@@ -217,9 +220,6 @@ struct lonlat
     /// @typedef Latitude and Longitude value type
     using decimal_degree = double;
 
-    /// @typedef The great circle distance value type
-    using km = double;
-
     /// @brief Longitude
     decimal_degree lon;
 
@@ -240,11 +240,13 @@ struct lonlat
     /// @param other the other coordinate
     /// @return the distance in kilometers
     /// @invariant The returned distance is not negative.
-    template <class T> constexpr km great_circle_distance_to(const T &other) const noexcept
+    template <class T> 
+    constexpr 
+    mp_units::quantity<mp_units::si::metre> great_circle_distance_to(const T &other) const noexcept
     {
-        const km d = details::distanceEarth(this->lat, this->lon, other.lat, other.lon);
-        assert(d >= 0.);
-        return d;
+        const double d = details::distanceEarth(this->lat, this->lon, other.lat, other.lon);
+        assert(d >= 0.0);
+        return d * 1000.0 * mp_units::si::metre;
     }
 
 }; // end struct latlon

test/unit_test/geography_test.cpp

@@ -13,6 +13,7 @@
 
 #include <filesystem>
 #include <iostream>
+#include <mp-units/ostream.h>
 
 namespace utf = boost::unit_test;
 
@@ -49,15 +50,17 @@ BOOST_AUTO_TEST_CASE(resolution)
 
 BOOST_AUTO_TEST_CASE(coordinates)
 {
+    using namespace mp_units;
+    using namespace mp_units::si::unit_symbols;
     using coord_type = quetzal::geography::latlon;
     coord_type paris(48.856614, 2.3522219000000177);
     coord_type moscow(55.7522200, 37.6155600);
     bool b = paris < moscow;
     BOOST_TEST(b);
-    coord_type::km computed = paris.great_circle_distance_to(moscow);
-    coord_type::km expected = 2486.22;
-    coord_type::km EPSILON = 1.; // a 1-km error is acceptable.
-    coord_type::km diff = expected - computed;
+    auto computed = paris.great_circle_distance_to(moscow);
+    auto expected = 2486.22 * km ;
+    auto EPSILON = 1.0 * km; // a 1-km error is acceptable.
+    auto diff = expected - computed;
     BOOST_TEST(diff < EPSILON);
     BOOST_TEST(-diff < EPSILON);
 }