diff --git a/docs/4-tutorials.md b/docs/4-tutorials.md
index 887618c7..c8092a05 100644
--- a/docs/4-tutorials.md
+++ b/docs/4-tutorials.md
@@ -29,9 +29,9 @@
   - @subpage Phylogenetic Trees
   - @subpage Phylogenetic Networks
   - @subpage spatial_graphs
-    - @subpage Construction from spatial grids
-    - @subpage Setting Vertex/Edge information at construction time
+    - @subpage spatial_graph_construction
     - @subpage dispersal_kernels
+    - @subpage Setting Vertex/Edge information at construction time
     - @subpage Growth Expressions
     - @subpage Memory Management
 
@@ -891,6 +891,10 @@ There are different ways to build such graph.
 
 ## Graph from raster
 
+The provided code snippet contains the logic for constructing a spatial graph from a `quetzal::geography::raster`. 
+
+This involves identifying the cells of the raster to map `location_type` values to the vertices, and defining edges based on the desired connectivity.
+
 **Input**
 
 @include{lineno} geography_graph_2.cpp
@@ -901,6 +905,10 @@ There are different ways to build such graph.
 
 ## Graph from landscape
 
+The provided code snippet contains the logic for constructing a spatial graph from a `quetzal::geography::landscape`. 
+
+This involves identifying the cells of the landscape to map `location_type` values to the vertices, and defining edges based on the desired connectivity.
+
 **Input**
 
 @include{lineno} geography_graph_3.cpp
@@ -911,6 +919,10 @@ There are different ways to build such graph.
 
 ## Graph from abstract grid
 
+The provided code snippet contains the logic for constructing a spatial graph from a user-defined abstract grid.
+
+This involves identifying the height and width of the space to define linearly-indexed vertices, and defining edges based on the desired connectivity.
+
 **Input**
 
 @include{lineno} geography_graph_4.cpp
@@ -921,6 +933,9 @@ There are different ways to build such graph.
 
 ## From scratch
 
+The provided code snippet contains the logic for creating a sparse versus dense graph from scratch. This involves 
+dynamically adding vertices and edges based on user-defined conditions, resulting in a graph that represents the desired spatial relationships.
+
 ### Sparse Graph
 
 **Input**
@@ -944,7 +959,7 @@ There are different ways to build such graph.
 ---
 
 [//]: # (----------------------------------------------------------------------)
-@page dispersal_kernels Dispersal Kernels
+@page dispersal_kernels Distance-based dispersal Kernels
 @tableofcontents
 
 ## Distance-based kernel
@@ -969,5 +984,16 @@ Quetzal incorporates various kernel types available in the `quetzal::demography:
 
 ---
 
-## Resistance-based Dispersal
+[//]: # (----------------------------------------------------------------------)
+@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
 
+---
\ No newline at end of file
diff --git a/example/expressive_3.cpp b/example/expressive_3.cpp
index 25489f90..dc04bd73 100644
--- a/example/expressive_3.cpp
+++ b/example/expressive_3.cpp
@@ -22,7 +22,7 @@ int main()
     std::iota(times.begin(), times.end(), 2001);
 
     // Makes sure the rasters are aligned
-    auto landscape = landscape_type::from_files({{"suit", file1}, {"DEM", file2}}, times);
+    auto landscape = landscape_type::from_file({{"suit", file1}, {"DEM", file2}}, times);
 
     // lightweight functor returning empty optionals where NA
     auto s_view = landscape["suit"].to_view();
diff --git a/example/expressive_4.cpp b/example/expressive_4.cpp
index 1e03c9e3..814e08d4 100644
--- a/example/expressive_4.cpp
+++ b/example/expressive_4.cpp
@@ -22,7 +22,7 @@ int main()
     std::iota(times.begin(), times.end(), 2001);
 
     // Makes sure the rasters are aligned
-    auto landscape = landscape_type::from_files({{"suit", file1}, {"DEM", file2}}, times);
+    auto landscape = landscape_type::from_file({{"suit", file1}, {"DEM", file2}}, times);
 
     // lightweight functor returning empty optionals where NA
     auto s_view = landscape["suit"].to_view();
diff --git a/example/geography_graph_1.cpp b/example/geography_graph_1.cpp
index 16bd0d42..5b7fdfc0 100644
--- a/example/geography_graph_1.cpp
+++ b/example/geography_graph_1.cpp
@@ -17,7 +17,7 @@ int main()
 
     // Initialize the landscape: for each var a key and a file, for all a time series.
     using landscape_type = geo::landscape<>;
-    auto land = landscape_type::from_files({{"bio1", file1}, {"bio12", file2}}, times);
+    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
     // Our graph will not store any useful information
     using vertex_info = geo::no_property;
diff --git a/example/geography_graph_3.cpp b/example/geography_graph_3.cpp
index 1b7d353a..bf3e7044 100644
--- a/example/geography_graph_3.cpp
+++ b/example/geography_graph_3.cpp
@@ -16,7 +16,7 @@ int main()
 
     // Initialize the landscape: for each var a key and a file, for all a time series.
     using landscape_type = geo::landscape<>;
-    auto land = landscape_type::from_files({{"bio1", file1}, {"bio12", file2}}, times);
+    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
     // Our graph will not store any useful information
     using vertex_info = geo::no_property;
diff --git a/example/geography_graph_7.cpp b/example/geography_graph_7.cpp
new file mode 100644
index 00000000..9cd1ec1e
--- /dev/null
+++ b/example/geography_graph_7.cpp
@@ -0,0 +1,37 @@
+#include "quetzal/geography.hpp"
+
+namespace geo = quetzal::geography;
+
+struct MyVertexInfo {
+    std::string label = "NA";
+    std::vector<double> pop_data_chunk;
+};
+
+struct MyEdgeInfo {
+    double distance;
+    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) );
+    }
+};
+
+int main()
+{
+    auto file = std::filesystem::current_path() / "data/bio1.tif";
+
+    // The raster have 10 bands that we will assign to 2001 ... 2010.
+    std::vector<int> times(10);
+    std::iota(times.begin(), times.end(), 2001);
+
+    using landscape_type = geo::raster<>;
+    auto land = geo::raster<>::from_file(file, times);
+    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;
+
+    for( auto const& e : graph.edges() ){
+        std::cout << graph.source(e) << " <-> " << graph.target(e) << " : " << graph[e].distance << std::endl;
+    }
+
+}
diff --git a/example/geography_graph_complete_1.cpp b/example/geography_graph_complete_1.cpp
index 415be831..8dc00cd4 100644
--- a/example/geography_graph_complete_1.cpp
+++ b/example/geography_graph_complete_1.cpp
@@ -19,7 +19,7 @@ int main()
 
     // Initialize the landscape: for each var a key and a file, for all a time series.
     using landscape_type = geo::landscape<>;
-    auto land = landscape_type::from_files({{"bio1", file1}, {"bio12", file2}}, times);
+    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
     // Our graph will not store any useful information
     using vertex_info = geo::no_property;
@@ -32,7 +32,7 @@ int main()
     // Define a helper function to compute distance on Earth between two location_descriptors
     auto sphere_distance = [&](auto point1, auto point2)
     {
-        return land.to_latlon(point1)->great_circle_distance_to(land.to_latlon(point2).value());
+        return land.to_latlon(point1).great_circle_distance_to(land.to_latlon(point2));
     };
 
     // Define the type of distance-based kernel to use
diff --git a/example/geography_landscape_1.cpp b/example/geography_landscape_1.cpp
index 016dd305..a3106b06 100644
--- a/example/geography_landscape_1.cpp
+++ b/example/geography_landscape_1.cpp
@@ -15,7 +15,7 @@ int main()
 
     // Initialize the landscape: for each var a key and a file, for all a time series.
     using landscape_type = quetzal::geography::landscape<>;
-    auto env = quetzal::geography::landscape<>::from_files({{"bio1", file1}, {"bio12", file2}}, times);
+    auto env = quetzal::geography::landscape<>::from_file({{"bio1", file1}, {"bio12", file2}}, times);
     std::cout << env << std::endl;
 
     // We indeed recorded 2 variables: bio1 and bio12
diff --git a/example/geography_raster_1.cpp b/example/geography_raster_1.cpp
index 403863f0..c4000fc8 100644
--- a/example/geography_raster_1.cpp
+++ b/example/geography_raster_1.cpp
@@ -70,8 +70,8 @@ int main()
     assert(point_1 == point_2);
 
     // Defines a lambda expression for checking extent
-    auto check = [&bio1](auto x) {
-        std::cout << "Point " << x << " is " << (bio1.contains(x) ? " " : "not ") << "in bio1 extent" << std::endl;
+    auto check = [&](auto x) {
+        std::cout << "Point " << x << " is" << (bio1.contains(x) ? " " : " not ") << "in bio1 extent" << std::endl;
     };
 
     check(point_1);
diff --git a/src/include/quetzal/geography/graph/detail/bound_policy.hpp b/src/include/quetzal/geography/graph/detail/bound_policy.hpp
index 630a4179..b6e5078c 100644
--- a/src/include/quetzal/geography/graph/detail/bound_policy.hpp
+++ b/src/include/quetzal/geography/graph/detail/bound_policy.hpp
@@ -56,7 +56,7 @@ class sink
         int sink = num_land_vertices;
         assert( graph.num_vertices() == sink + 1 );
         if ( s < sink )
-            detail::edge_construction<edge_property>::delegate(s, sink, graph); // one-way ticket :(
+            detail::edge_construction<edge_property>::delegate(s, sink, graph, grid); // one-way ticket :(
     }
 };
 
@@ -100,7 +100,7 @@ class torus
         {
             symmetricIndex = s - width + 1;
         }
-        detail::edge_construction<edge_property>::delegate(s, symmetricIndex, graph);
+        detail::edge_construction<edge_property>::delegate(s, symmetricIndex, graph, grid);
     }
 };
 
diff --git a/src/include/quetzal/geography/graph/detail/concepts.hpp b/src/include/quetzal/geography/graph/detail/concepts.hpp
index c26160da..40606980 100644
--- a/src/include/quetzal/geography/graph/detail/concepts.hpp
+++ b/src/include/quetzal/geography/graph/detail/concepts.hpp
@@ -29,16 +29,16 @@ concept is_any_of = (std::same_as<T, U> || ...);
 template<typename T>
 struct edge_construction
 {
-    static inline constexpr void delegate(auto s, auto t, auto& graph)
+    static inline constexpr void delegate(auto s, auto t, auto& graph, auto const& grid)
     {
-        graph.add_edge(s, t, T(s,t));
+        graph.add_edge(s, t, T(s,t, grid));
     }
 };
 
 template<>
 struct edge_construction<boost::no_property>
 {
-    static inline constexpr void delegate(auto s, auto t, auto& graph)
+    static inline constexpr void delegate(auto s, auto t, auto& graph, [[maybe_unused]] auto const& grid)
     {
         graph.add_edge(s, t);
     }
diff --git a/src/include/quetzal/geography/graph/detail/vicinity.hpp b/src/include/quetzal/geography/graph/detail/vicinity.hpp
index 2f1196f5..7a0753e7 100644
--- a/src/include/quetzal/geography/graph/detail/vicinity.hpp
+++ b/src/include/quetzal/geography/graph/detail/vicinity.hpp
@@ -70,9 +70,9 @@ struct connect_fully
         {
             for (auto t = s + 1; t < num_land_vertices; ++t)
             {
-                detail::edge_construction<edge_property>::delegate(s, t, graph); // (s -> v) == (s <- v) if undirected
+                detail::edge_construction<edge_property>::delegate(s, t, graph, grid); // (s -> v) == (s <- v) if undirected
                 if constexpr (std::same_as<directed_category, anisotropy>)
-                    detail::edge_construction<edge_property>::delegate(t, s, graph); // (s -> v) != (s <- v) because directed
+                    detail::edge_construction<edge_property>::delegate(t, s, graph, grid); // (s -> v) != (s <- v) because directed
             }
         
             if( detail::is_on_border(s, width, height ) ) 
@@ -90,9 +90,9 @@ class connect_4_neighbors
     template <typename G> void connect(auto s, auto t, G &graph, auto const &grid) const
     {
         using edge_property = typename G::edge_property;
-        detail::edge_construction<edge_property>::delegate(s, t, graph);
+        detail::edge_construction<edge_property>::delegate(s, t, graph, grid);
         if constexpr (std::same_as<typename G::directed_category, anisotropy>)
-            detail::edge_construction<edge_property>::delegate(t, s, graph);
+            detail::edge_construction<edge_property>::delegate(t, s, graph, grid);
     }
 
     void connect_top_left_corner(auto s, auto &graph, auto const &grid, auto const &bound_policy) const
@@ -240,6 +240,7 @@ class connect_8_neighbors
         using vertex_property = typename G::vertex_property;
         using edge_property = typename G::edge_property;
 
+
         int width = grid.width();
         int height = grid.height();
         int num_land_vertices = width * height;
@@ -248,9 +249,9 @@ class connect_8_neighbors
         assert( t >= 0 );
         assert( t < num_land_vertices);
 
-        detail::edge_construction<edge_property>::delegate(s, t, graph);
+        detail::edge_construction<edge_property>::delegate(s, t, graph, grid);
         if constexpr (std::same_as<directed_category, anisotropy>)
-            detail::edge_construction<edge_property>::delegate(t, s, graph);
+            detail::edge_construction<edge_property>::delegate(t, s, graph, grid);
     }
 
     void connect_top_left_corner(auto s, auto &graph, auto const &grid, auto const &bound_policy) const
diff --git a/src/include/quetzal/geography/landscape.hpp b/src/include/quetzal/geography/landscape.hpp
index 696f3dd4..edd7d018 100644
--- a/src/include/quetzal/geography/landscape.hpp
+++ b/src/include/quetzal/geography/landscape.hpp
@@ -24,18 +24,18 @@
 namespace quetzal::geography
 {
 /// @brief Discrete spatio-temporal variations of a set of environmental variables.
-/// @details  Read and write multiple georeferenced raster datasets with multiple bands (layers).
+/// @details  Read and write multiple geo-referenced raster datasets with multiple bands (layers).
 ///           Multiple quetzal::geography::raster objects are assembled into a multivariate discrete
-///           quetzal::geography::landscape. This geospatial collection is constructed from raster files and gives
+///           quetzal::geography::landscape. This geo-spatial collection is constructed from raster files and gives
 ///           strong guarantees on grids consistency.
 /// @remark Working with multiple variables can lead to simulation problems if the rasters have even slightly different
-/// grid properties,
+///         grid properties,
 ///         like different resolutions, origins or extent. To prevent these problems from happening, a
 ///         `quetzal::geography::landscape` wraps multiple quetzal::geography::raster objects into one single coherent
 ///         object and checks that all rasters are aligned.
 /// @remark A quetzal::geography::landscape shares a large part of its spatial semantics with a
 /// quetzal::geography::raster.
-/// @tparam Key A key used to uniquely identifie a variable, e.g. std::string.
+/// @tparam Key A key used to uniquely identify a variable, e.g. std::string.
 /// @tparam Time Type used as time period for every band, e.g. std::string with `4.2-0.3 ka`
 /// @ingroup geography
 template <typename Key = std::string, typename Time = int> class landscape
@@ -132,7 +132,7 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @return A raster object
     /// @remark error will be thrown if the datasets do not have the exactly same geographical properties (origin,
     /// extent, resolution, number of layers).
-    static landscape from_files(const std::map<key_type, std::filesystem::path> &files,
+    inline static landscape from_file(const std::map<key_type, std::filesystem::path> &files,
                                 const std::vector<time_type> &times)
     {
         return landscape(files, times);
@@ -146,7 +146,7 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @param file the file name where to save
     template <typename Callable>
         requires std::invocable<Callable, location_descriptor, time_descriptor, std::optional<double>>
-    void to_geotiff(Callable f, time_descriptor start, time_descriptor end, const std::filesystem::path &file) const
+    inline void to_geotiff(Callable f, time_descriptor start, time_descriptor end, const std::filesystem::path &file) const
     {
         _variables.cbegin()->second.to_geotiff(f, start, end, file);
     }
@@ -154,7 +154,7 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @brief Export spatial points to a shapefile
     /// @param points the points to export
     /// @param file the path to write to
-    void to_shapefile(std::vector<latlon> points, const std::filesystem::path &file) const
+    inline void to_shapefile(std::vector<latlon> points, const std::filesystem::path &file) const
     {
         _variables.cbegin()->second.to_shapefile(points, file);
     }
@@ -162,7 +162,7 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @brief Export geolocalized counts as spatial-points to a shapefile
     /// @brief counts the number of points at every point
     /// @brief file the path to write to.
-    void to_shapefile(std::map<latlon, int> counts, const std::filesystem::path &file) const
+    inline void to_shapefile(std::map<latlon, int> counts, const std::filesystem::path &file) const
     {
         _variables.cbegin()->second.to_shapefile(counts, file);
     }
@@ -190,14 +190,14 @@ template <typename Key = std::string, typename Time = int> class landscape
 
     /// @brief Retrieves the number of variables (rasters)
     /// @return the number of variables
-    auto num_variables() const noexcept
+    inline auto num_variables() const noexcept
     {
         return _variables.size();
     }
 
     /// @brief Number of cells in the raster
     /// @return the number of variables
-    auto num_locations() const noexcept
+    inline auto num_locations() const noexcept
     {
         return this->width() * this->height();
     }
@@ -211,7 +211,7 @@ template <typename Key = std::string, typename Time = int> class landscape
     ///        If no such raster exists, an exception of type std::out_of_range is thrown.
     /// @param key the identifier of the variable to be retrieved
     /// @return the raster object of the specified variable.
-    const raster<time_type> &operator[](const key_type &key) const
+    inline const raster<time_type> &operator[](const key_type &key) const
     {
         return _variables.at(key);
     }
@@ -220,7 +220,7 @@ template <typename Key = std::string, typename Time = int> class landscape
     ///        If no such raster exists, an exception of type std::out_of_range is thrown.
     /// @param key the identifier of the variable to be retrieved
     /// @return the raster object of the specified variable.
-    raster<time_type> &operator[](const key_type &key)
+    inline raster<time_type> &operator[](const key_type &key)
     {
         return _variables.at(key);
     }
@@ -231,37 +231,37 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @{
 
     /// @brief Origin of the spatial grid
-    latlon origin() const noexcept
+    inline latlon origin() const noexcept
     {
         return _variables.cbegin()->second.origin();
     }
 
     /// @brief Resolution of the spatial grid
-    resolution<decimal_degree> get_resolution() const noexcept
+    inline resolution<decimal_degree> get_resolution() const noexcept
     {
         return _variables.cbegin()->second.get_resolution();
     }
 
     /// @brief Extent of the spatial grid
-    extent<decimal_degree> get_extent() const noexcept
+    inline extent<decimal_degree> get_extent() const noexcept
     {
         return _variables.cbegin()->second.get_extent();
     }
 
     /// @brief Width of the spatial grid
-    int width() const noexcept
+    inline int width() const noexcept
     {
         return _variables.cbegin()->second.width();
     }
 
     /// @brief Height of the spatial grid
-    int height() const noexcept
+    inline int height() const noexcept
     {
         return _variables.cbegin()->second.height();
     }
 
     /// @brief Depth of the spatial grid
-    int depth() const noexcept
+    inline int depth() const noexcept
     {
         return _variables.cbegin()->second.depth();
     }
@@ -272,33 +272,81 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @{
 
     /// @brief Location descriptors (unique identifiers) of the grid cells
-    auto locations() const noexcept
+    /// @return A range of unique identifiers for every cell
+    inline auto locations() const noexcept
     {
         return _variables.cbegin()->second.locations();
     }
 
     /// @brief Time descriptors (unique identifiers) of the dataset bands
-    auto times() const noexcept
+    /// @return A range of unique identifiers for every time point recorded
+    inline auto times() const noexcept
     {
         return _variables.cbegin()->second.times();
     }
 
-    ///@brief checks if the raster contains a layer with specific time
-    bool contains(const latlon &x) const noexcept
+    /// @brief Check if a descriptor describes a valid location of the spatial grid.
+    /// @param x the descriptor to check
+    /// @return True if x is valid, false otherwise.
+    inline bool is_valid(location_descriptor x) const noexcept
+    {
+        return _variables.cbegin()->second.is_valid(x);
+    }
+
+    /// @brief Check if the coordinate describes a valid location of the spatial grid.
+    /// @param x the column/row to check.
+    /// @return True if x is valid, false otherwise.
+    inline bool is_valid(const colrow &x) const noexcept
+    {
+        return _variables.cbegin()->second.is_valid(x);
+    }
+
+    /// @brief Check if the coordinate describes a valid location of the spatial grid.
+    /// @param x the row/column to check.
+    /// @return True if x is valid, false otherwise.
+    inline bool is_valid(const rowcol &x) const noexcept
+    {
+        return _variables.cbegin()->second.is_valid(x);
+    }
+
+    /// @brief Check if the raster contains a coordinate
+    /// @param x the lat/lon to check.
+    /// @return True if x is valid, false otherwise.
+    inline bool contains(const latlon &x) const noexcept
+    {
+        return _variables.cbegin()->second.contains(x);
+    }
+
+    /// @brief Check if the raster contains a coordinate
+    /// @param x the lon/lat to check.
+    /// @return True if x is valid, false otherwise.
+    inline bool contains(const lonlat &x) const noexcept
     {
         return _variables.cbegin()->second.contains(x);
     }
 
-    ///@brief checks if the raster contains a layer with specific time
-    bool contains(const lonlat &x) const noexcept
+    /// @brief Check if the time descriptor is a valid index.
+    /// @param t A time descriptor identifying a layer of the raster.
+    /// @return True if the descriptor is a valid index, false otherwise.
+    inline bool is_valid(time_descriptor t) const noexcept
+    {
+        return _variables.cbegin()->second.is_valid(t);
+    }
+
+    /// @brief Search for the exact time in the list of time points recorded by the raster.
+    /// @param t An exact point in time
+    /// @return True if the time point is found, false otherwise.
+    inline bool is_recorded(const time_type &t) const noexcept
     {
-        return _variables.cbegin()->second.contains(latlon(x.lat, x.lon));
+        return _variables.cbegin()->second.is_recorded(t);
     }
 
-    ///@brief checks if the raster contains a layer with specific time
-    bool contains(const time_type &t) const noexcept
+    /// @brief Check if the exact time point falls between the first and last date of record.
+    /// @param t An exact point in time
+    /// @return True if the time point falls in the temporal range, false otherwise.
+    inline bool is_in_interval(const time_type &t) const noexcept
     {
-        return _variables.cbegin()->second.contains(t);
+        return _variables.cbegin()->second.is_in_interval(t);
     }
 
     /// @}
@@ -307,66 +355,92 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @{
 
     /// @brief Location descriptor of the cell identified by its column/row.
-    /// @param x the latitude/longitude coordinate to evaluate.
-    /// @return An optional with the descriptor value if x is in the spatial extent, empty otherwise.
-    std::optional<location_descriptor> to_descriptor(const colrow &x) const noexcept
+    /// @param x the coordinate to evaluate.
+    /// @return The descriptor value.
+    /// @pre x is in spatial extent.
+    inline location_descriptor to_descriptor(const colrow &x) const noexcept
     {
+        assert( is_valid(x) );
         return _variables.cbegin()->second.to_descriptor(x);
     }
 
     /// @brief Location descriptor of the cell to which the given coordinate belongs.
-    /// @param x the latitude/longitude coordinate to evaluate.
-    /// @return An optional with the descriptor value if x is in the spatial extent, empty otherwise.
-    std::optional<location_descriptor> to_descriptor(const latlon &x) const noexcept
+    /// @param x the coordinate to evaluate.
+    /// @return The descriptor value.
+    /// @pre x is in spatial extent.
+    inline location_descriptor to_descriptor(const latlon &x) const noexcept
     {
+        assert( contains(x) );
         return _variables.cbegin()->second.to_descriptor(x);
     }
 
     /// @brief Column and row of the cell to which the given coordinate belongs.
-    /// @param x the latitude/longitude coordinate to evaluate.
-    /// @return An optional with a pair (column, row) if x is in the spatial extent, empty otherwise.
-    std::optional<colrow> to_colrow(const latlon &x) const noexcept
+    /// @param x the coordinate to evaluate.
+    /// @return A pair (column, row).
+    /// @pre x is in spatial extent.
+    colrow to_colrow(const latlon &x) const noexcept
     {
+        assert( contains(x) );
         return _variables.cbegin()->second.to_colrow(x);
     }
 
     /// @brief Row and column of the cell to which the given coordinate belongs.
-    /// @param x the latitude/longitude coordinate to evaluate.
-    /// @return An optional with a pair (row, column) if x is in the spatial extent, empty otherwise.
-    std::optional<rowcol> to_rowcol(const latlon &x) const noexcept
+    /// @param x the coordinate to evaluate.
+    /// @return A pair (row, column).
+    /// @pre x is in spatial extent.
+    inline rowcol to_rowcol(const latlon &x) const noexcept
     {
+        assert( contains(x) );
         return _variables.cbegin()->second.to_rowcol(x);
     }
 
     /// @brief Column and row of the cell to which the given descriptor belongs.
     /// @param x the location to evaluate.
-    /// @return An optional with a pair (column, row) if x is in the spatial extent, empty otherwise.
-    std::optional<colrow> to_colrow(location_descriptor x) const noexcept
+    /// @return A pair (column, row).
+    /// @pre x is in spatial extent.
+    inline colrow to_colrow(location_descriptor x) const noexcept
     {
+        assert( is_valid(x) );
         return _variables.cbegin()->second.to_colrow(x);
     }
 
     /// @brief Latitude and longitude of the cell to which the given coordinate belongs.
     /// @param x the location to evaluate.
-    /// @return An optional with a pair (lat, lon) if x is in the spatial extent, empty otherwise.
-    std::optional<latlon> to_latlon(location_descriptor x) const noexcept
+    /// @return A pair (lat, lon).
+    /// @pre x is in spatial extent.
+    inline latlon to_latlon(location_descriptor x) const noexcept
     {
+        assert( is_valid(x) );
         return _variables.cbegin()->second.to_latlon(x);
     }
 
+    /// @brief Longitude and latitude of the cell to which the given coordinate belongs.
+    /// @param x the location to evaluate.
+    /// @return A pair (lon, lat).
+    /// @pre x is in spatial extent
+    inline lonlat to_lonlat(location_descriptor x) const noexcept
+    {
+        assert( is_valid(x) );
+        return _variables.cbegin()->second.to_lonlat(x);
+    }
+    
     /// @brief Latitude and longitude of the cell identified by its column/row.
     /// @param x the location to evaluate.
-    /// @return A pair (latitude, longitude)
-    latlon to_latlon(const colrow &x) const noexcept
+    /// @return A pair (latitude, longitude).
+    /// @pre x is in spatial extent.
+    inline latlon to_latlon(const colrow &x) const noexcept
     {
+        assert( is_valid(x) );
         return _variables.cbegin()->second.to_latlon(x);
     }
 
     /// @brief Longitude and latitude of the deme identified by its column/row.
     /// @param x the location to evaluate.
-    /// @return A pair (latitude, longitude)
-    lonlat to_lonlat(const colrow &x) const noexcept
+    /// @return A pair (longitude, latitude)
+    /// @pre x is in spatial extent
+    inline lonlat to_lonlat(const colrow &x) const noexcept
     {
+        assert( is_valid(x) );
         return _variables.cbegin()->second.to_lonlat(x);
     }
 
@@ -374,7 +448,8 @@ template <typename Key = std::string, typename Time = int> class landscape
     /// @pre The landscape must contain the coordinate in its spatial extent.
     /// @param x The location to reproject
     /// @return The coordinate of the centroid of the cell it belongs.
-    latlon to_centroid(const latlon &x) const
+    /// @pre x is in spatial extent
+    inline latlon to_centroid(const latlon &x) const
     {
         assert(this->contains(x));
         return _variables.cbegin()->second.to_centroid(x);
diff --git a/src/include/quetzal/geography/raster.hpp b/src/include/quetzal/geography/raster.hpp
index 9c0e0930..a1fe6670 100644
--- a/src/include/quetzal/geography/raster.hpp
+++ b/src/include/quetzal/geography/raster.hpp
@@ -38,11 +38,14 @@ template <typename Time = int> class raster
 {
 
   public:
-    /// @typedef Location unique identifier
-    using location_descriptor = int;
+    /// @typedef Time period for every band.
+    using time_type = Time;
 
     /// @typedef Time unique identifier
-    using time_descriptor = int;
+    using time_descriptor = std::vector<time_type>::size_type;
+
+    /// @typedef Location unique identifier
+    using location_descriptor = int;
 
     /// @typedef A Latitude/Longitude pair of value.
     using latlon = quetzal::geography::latlon;
@@ -56,9 +59,6 @@ template <typename Time = int> class raster
     /// @typedef A row/column pair of indices.
     using rowcol = quetzal::geography::rowcol;
 
-    /// @typedef Time period for every band.
-    using time_type = Time;
-
     /// @typedef decimal degree type
     using decimal_degree = double;
 
@@ -79,13 +79,13 @@ template <typename Time = int> class raster
     value_type _no_data;
 
     // Pixel size in x and y
-    auto compute_resolution(const std::vector<double> &gT) const
+    inline auto compute_resolution(const std::vector<double> &gT) const
     {
         return quetzal::geography::resolution<decimal_degree>(gT[5], gT[1]);
     }
 
     // Top left corner and bottom right corner geographic coordinates
-    auto compute_extent(const latlon &origin, const quetzal::geography::resolution<decimal_degree> &r) const
+    inline  auto compute_extent(const latlon &origin, const quetzal::geography::resolution<decimal_degree> &r) const
     {
         const auto [lat_max, lon_min] = origin;
         const auto lon_max = lon_min + width() * r.lon();
@@ -112,11 +112,9 @@ template <typename Time = int> class raster
 
     /// @brief Read the raster from an input file in the geotiff format
     /// @param file The file to read from
-    /// @param start The first time period to map to the first band
-    /// @param end The last time period to map to the last band
-    /// @param step The time increment between each band
+    /// @param times The time periods relative to each band.
     /// @return A raster object
-    static raster from_file(const std::filesystem::path &file, const std::vector<time_type> &times)
+    inline static raster from_file(const std::filesystem::path &file, const std::vector<time_type> &times)
     {
         return raster(file, times);
     }
@@ -130,7 +128,7 @@ template <typename Time = int> class raster
     /// @remark If the returned optional is empty, the raster NA value will be used.
     template <class Callable>
         requires(std::is_invocable_r_v<std::optional<value_type>, Callable, location_descriptor, time_type>)
-    void to_geotiff(Callable f, time_type start, time_type end, const std::filesystem::path &file) const
+    inline void to_geotiff(Callable f, time_type start, time_type end, const std::filesystem::path &file) const
     {
         export_to_geotiff(f, start, end, file);
     }
@@ -138,7 +136,7 @@ template <typename Time = int> class raster
     /// @brief Export spatial points to a shapefile
     /// @param points the points to export
     /// @param file the path to write to
-    void to_shapefile(std::vector<latlon> points, const std::filesystem::path &file) const
+    inline void to_shapefile(std::vector<latlon> points, const std::filesystem::path &file) const
     {
         points_to_shapefile(points, file);
     }
@@ -146,7 +144,7 @@ template <typename Time = int> class raster
     /// @brief Export geolocalized counts as spatial-points to a shapefile
     /// @brief counts the number of points at every point
     /// @brief file the path to write to.
-    void to_shapefile(std::map<latlon, int> counts, const std::filesystem::path &file) const
+    inline void to_shapefile(std::map<latlon, int> counts, const std::filesystem::path &file) const
     {
         counts_to_shapefile(counts, file);
     }
@@ -169,44 +167,44 @@ template <typename Time = int> class raster
     /// @{
 
     /// @brief Origin of the spatial grid
-    latlon origin() const noexcept
+    inline latlon origin() const noexcept
     {
         return _origin;
     }
 
     /// @brief Resolution of the spatial grid
     /// @note The vertical pixel size will be negative (south pointing)
-    geography::resolution<decimal_degree> get_resolution() const noexcept
+    inline geography::resolution<decimal_degree> get_resolution() const noexcept
     {
         return _resolution;
     }
 
     /// @brief Extent of the spatial grid
-    geography::extent<decimal_degree> get_extent() const noexcept
+    inline geography::extent<decimal_degree> get_extent() const noexcept
     {
         return _extent;
     }
 
     /// @brief Width of the spatial grid
-    int width() const noexcept
+    inline int width() const noexcept
     {
         return _width;
     }
 
     /// @brief Height of the spatial grid
-    int height() const noexcept
+    inline int height() const noexcept
     {
         return _height;
     }
 
     /// @brief Depth of the spatial grid
-    int depth() const noexcept
+    inline int depth() const noexcept
     {
         return _depth;
     }
 
     /// @brief NA value, as read in the first band of the dataset.
-    value_type NA() const noexcept
+    inline value_type NA() const noexcept
     {
         return _no_data;
     }
@@ -217,44 +215,94 @@ template <typename Time = int> class raster
     /// @{
 
     /// @brief Location descriptors (unique identifiers) of the grid cells
-    auto locations() const noexcept
+    /// @return A range of unique identifiers for every cell
+    inline auto locations() const noexcept
     {
         return ranges::views::iota(0, width() * height());
     }
 
     /// @brief Time descriptors (unique identifiers) of the dataset bands
-    auto times() const noexcept
+    /// @return A range of unique identifiers for every time point recorded
+    inline auto times() const noexcept
     {
         // [0 ... depth -1 [
         return ranges::views::iota(0, depth());
     }
 
-    ///@brief checks if the raster contains a layer with specific time
-    bool contains(const latlon &x) const noexcept
+    /// @brief Check if a descriptor describes a valid location of the spatial grid.
+    /// @param x the descriptor to check
+    /// @return True if x is valid, false otherwise.
+    inline bool is_valid(location_descriptor x) const noexcept
+    {
+        return ( x >= 0 and x < locations().size());
+    }
+
+    /// @brief Check if the coordinate describes a valid location of the spatial grid.
+    /// @param x the column/row to check.
+    /// @return True if x is valid, false otherwise.
+    inline bool is_valid(const colrow &x) const noexcept
+    {
+        return x.col >= 0 and x.col < width() and x.row >= 0 and x.row < height() ;
+    }
+
+    /// @brief Check if the coordinate describes a valid location of the spatial grid.
+    /// @param x the row/column to check.
+    /// @return True if x is valid, false otherwise.
+    inline bool is_valid(const rowcol &x) const noexcept
+    {
+        return x.col >= 0 and x.col < width() and x.row >= 0 and x.row < height() ;
+    }
+
+    /// @brief Check if the raster contains a coordinate
+    /// @param x the lat/lon to check.
+    /// @return True if x is valid, false otherwise.    
+    inline bool contains(const latlon &x) const noexcept
     {
         return is_in_spatial_extent(x);
     }
 
-    ///@brief checks if the raster contains a layer with specific time
-    bool contains(const lonlat &x) const noexcept
+    /// @brief Check if the raster contains a coordinate
+    /// @param x the lon/lat to check.
+    /// @return True if x is valid, false otherwise.
+    inline bool contains(const lonlat &x) const noexcept
     {
         return is_in_spatial_extent(to_latlon(x));
     }
 
-    ///@brief checks if the raster contains a layer with specific time
-    bool contains(const time_type &t) const noexcept
+    /// @brief Check if the time descriptor is a valid index.
+    /// @param t A time descriptor identifying a band of the raster dataset.
+    /// @return True if t is a valid index, false otherwise.
+    inline bool is_valid(time_descriptor t) const noexcept
+    {
+        return t >= 0 and t < _times.size();
+    }
+
+    /// @brief Search for the exact time in the list of time points recorded by the raster.
+    /// @param t An exact point in time
+    /// @return True if the time point is found, false otherwise.
+    inline bool is_recorded(const time_type &t) const noexcept
+    {
+        return std::find(_times.begin(), _times.end(), t) != _times.end();
+    }
+
+    /// @brief Check if the exact time point falls between the first and last date of record.
+    /// @param t An exact point in time
+    /// @return True if the time point falls in the temporal range, false otherwise.
+    inline bool is_in_interval(const time_type &t) const noexcept
     {
-        return is_in_temporal_extent(t);
+        return t >= _times.front() and t <= _times.back();
     }
 
     /// @brief Read the value at location x and time t
     /// @param x the location
     /// @param t the time
     /// @return An optional that is empty if the value read is equal to NA.
-    std::optional<value_type> at(location_descriptor x, time_descriptor t) const noexcept
+    /// @pre x is contained in the spatial extent
+    /// @pre t is contained in the temporal extent
+    inline std::optional<value_type> at(location_descriptor x, time_descriptor t) const noexcept
     {
-        assert(x >= 0 and x < locations().size());
-        assert(t >= 0 and t < times().size());
+        assert( is_valid(x) );
+        assert( is_valid(t) );
         const auto colrow = to_colrow(x);
         return read(colrow.col, colrow.row, t);
     }
@@ -264,7 +312,7 @@ template <typename Time = int> class raster
     ///        returning the value at location x and time t
     /// @remark This callable is cheap to copy, does not have ownership of the data,
     ///         and works well for mathematical composition using `quetzal::expressive`
-    auto to_view() const noexcept
+    inline auto to_view() const noexcept
     {
         return [this](location_descriptor x, time_descriptor t) { return this->at(x, t); };
     }
@@ -274,38 +322,32 @@ template <typename Time = int> class raster
     /// @{
 
     /// @brief Location descriptor of the cell identified by its column/row.
-    ///        If no such cell exists, an exception of type std::out_of_range is thrown.
     /// @param x the 1D coordinate to evaluate.
-    /// @return The descriptor value of x
-    location_descriptor to_descriptor(const colrow &x) const
+    /// @return The descriptor value of x.
+    /// @pre x is in spatial extent.
+    inline location_descriptor to_descriptor(const colrow &x) const noexcept
     {
-        if (x.col >= width() || x.col < 0 || x.row < 0 || x.row >= height())
-            throw std::out_of_range("colrow has invalid indices and can not be converted to a location descriptor");
+        assert( is_valid(x) );
         return x.col * width() + x.col;
     }
 
     ///  @brief Location descriptor of the cell to which the given coordinate belongs.
-    ///        If no such cell exists, an exception of type std::out_of_range is thrown.
     /// @param x the latitude/longitude coordinate to evaluate.
-    /// @return The descriptor value of x
-    location_descriptor to_descriptor(const latlon &x) const
+    /// @return The descriptor value of x.
+    /// @pre x is in the spatial extent.
+    inline location_descriptor to_descriptor(const latlon &x) const noexcept
     {
-        if (!contains(x))
-            throw std::out_of_range(
-                "latlon does not belong to spatial extent and can not be converted to a location descriptor");
-        return to_descriptor(to_colrow(x));
+        assert( contains(x) );
+        return to_descriptor( to_colrow(x) );
     }
 
     /// @brief Column and row of the cell to which the given coordinate belongs.
-    ///        If no such cell exists, an exception of type std::out_of_range is thrown.
     /// @param x the latitude/longitude coordinate to evaluate.
-    /// @return A pair (column, row)
-    colrow to_colrow(const latlon &x) const
+    /// @return A pair (column, row).
+    /// @pre x is in the spatial extent.
+    inline colrow to_colrow(const latlon &x) const noexcept
     {
-        if (!contains(x))
-            throw std::out_of_range(
-                "latlon does not belong to spatial extent and can not be converted to a col/row pair");
-
+        assert(contains(x));
         const auto [lat, lon] = x;
         const int col = (lon - _gT[0]) / _gT[1];
         const int row = (lat - _gT[3]) / _gT[5];
@@ -313,15 +355,12 @@ template <typename Time = int> class raster
     }
 
     /// @brief Row and column of the cell to which the given coordinate belongs.
-    ///        If no such cell exists, an exception of type std::out_of_range is thrown.
     /// @param x the latitude/longitude coordinate to evaluate.
-    /// @return A pair (row, column)
-    rowcol to_rowcol(const latlon &x) const noexcept
+    /// @return A pair (row, column).
+    /// @pre x is in the spatial extent.
+    inline rowcol to_rowcol(const latlon &x) const noexcept
     {
-        if (!contains(x))
-            throw std::out_of_range(
-                "latlon does not belong to spatial extent and can not be converted to a row/col pair.");
-
+        assert(contains(x));
         auto colrow = to_colrow(x);
         std::swap(colrow.col, colrow.row);
         return colrow;
@@ -329,9 +368,11 @@ template <typename Time = int> class raster
 
     /// @brief Column and row of the cell to which the given descriptor belongs.
     /// @param x the location to evaluate.
-    /// @return A pair (column, row)
-    colrow to_colrow(location_descriptor x) const noexcept
+    /// @return A pair (column, row).
+    /// @pre x is in the spatial extent.
+    inline colrow to_colrow(location_descriptor x) const noexcept
     {
+        assert(is_valid(x));
         const int col = x % width();
         const int row = x / width();
         assert(row < height());
@@ -340,26 +381,42 @@ template <typename Time = int> class raster
 
     /// @brief Column and row of the cell to which the given descriptor belongs.
     /// @param x the location to evaluate.
-    /// @return A pair (column, row)
-    rowcol to_rowcol(location_descriptor x) const noexcept
+    /// @return A pair (column, row).
+    /// @pre x is in the spatial extent.
+    inline rowcol to_rowcol(location_descriptor x) const noexcept
     {
+        assert(is_valid(x));
         auto c = to_colrow(x);
         return rowcol(c.row, c.col);
     }
 
     /// @brief Latitude and longitude of the cell to which the given coordinate belongs.
     /// @param x the location to evaluate.
-    /// @return An optional with a pair (lat, lon) if x is in the spatial extent, empty otherwise.
-    latlon to_latlon(location_descriptor x) const noexcept
+    /// @return A pair (lat, lon).
+    /// @pre x is in the spatial extent.
+    inline latlon to_latlon(location_descriptor x) const noexcept
     {
+        assert(is_valid(x));
         return to_latlon(to_colrow(x));
     }
 
+    /// @brief Latitude and longitude of the cell to which the given coordinate belongs.
+    /// @param x the location to evaluate.
+    /// @return A pair (lon, lat).
+    /// @pre x is in the spatial extent
+    inline lonlat to_lonlat(location_descriptor x) const noexcept
+    {
+        assert(is_valid(x));
+        return to_lonlat(to_colrow(x));
+    }
+
     /// @brief Latitude and longitude of the cell identified by its column/row.
     /// @param x the location to evaluate.
     /// @return A pair (latitude, longitude)
-    latlon to_latlon(const colrow &x) const noexcept
+    /// @pre x is in the spatial extent
+    inline latlon to_latlon(const colrow &x) const noexcept
     {
+        assert(is_valid(x));
         const auto [col, row] = x;
         const auto lon = _gT[1] * col + _gT[2] * row + _gT[1] * 0.5 + _gT[2] * 0.5 + _gT[0];
         const auto lat = _gT[4] * col + _gT[5] * row + _gT[4] * 0.5 + _gT[5] * 0.5 + _gT[3];
@@ -368,50 +425,31 @@ template <typename Time = int> class raster
 
     /// @brief Latitude and longitude of the cell identified by its row/column.
     /// @param x the location to evaluate.
-    /// @return A pair (latitude, longitude)
-    latlon to_latlon(const rowcol &x) const noexcept
+    /// @return A pair (latitude, longitude).
+    /// @pre x is in the spatial extent.
+    inline latlon to_latlon(const rowcol &x) const noexcept
     {
+        assert(is_valid(x));
         return to_latlon(colrow(x.col, x.row));
     }
 
     /// @brief Longitude and latitude of the deme identified by its column/row.
     /// @param x the location to evaluate.
-    /// @return A pair (latitude, longitude)
-    lonlat to_lonlat(const colrow &x) const noexcept
+    /// @return A pair (longitude, latitude).
+    /// @pre x is in the spatial extent.
+    inline lonlat to_lonlat(const colrow &x) const noexcept
     {
+        assert(is_valid(x));
         const auto [lat, lon] = to_latlon(x);
         return lonlat(lon, lat);
     }
 
-    /// @brief Latlon to lonlat conversation
-    /// @param x the location to evaluate.
-    /// @return A pair (longitude, latitude)
-    lonlat to_lonlat(const latlon &x) const noexcept
-    {
-        return lonlat(x.lon, x.lat);
-    }
-
-    /// @brief Latlon to lonlat conversation
-    /// @param x the location to evaluate.
-    /// @return A pair (longitude, latitude)
-    latlon to_latlon(const lonlat &x) const noexcept
-    {
-        return latlon(x.lat, x.lon);
-    }
-
-    /// @brief Longitude and latitude of the deme identified by its column/row.
-    /// @param x the location to evaluate.
-    /// @return A pair (latitude, longitude)
-    lonlat to_lonlat(location_descriptor x) const noexcept
-    {
-        return to_lonlat(to_latlon(x));
-    }
-
     /// @brief Reprojects a coordinate to the centroid of the cell it belongs.
     /// @pre The spatial extent of the raster must contain the coordinate.
     /// @param x The location to reproject
     /// @return The coordinate of the centroid of the cell it belongs.
-    latlon to_centroid(const latlon &x) const
+    /// @pre x is in the spatial extent
+    inline latlon to_centroid(const latlon &x) const
     {
         assert(this->contains(x) and 
                 "latlon does not belong to spatial extent and can not be reprojected to a cell centroid.");
@@ -422,14 +460,25 @@ template <typename Time = int> class raster
     /// @}
 
   private:
-    /// @brief Is time in the temporal extent
-    bool is_in_temporal_extent(time_type t) const noexcept
+
+    /// @brief Latlon to lonlat conversation
+    /// @param x the location to evaluate.
+    /// @return A pair (longitude, latitude)
+    inline lonlat to_lonlat(const latlon &x) const noexcept
     {
-        return t >= _times.front() && t <= _times.back();
+        return lonlat(x.lon, x.lat);
     }
 
+    /// @brief Latlon to lonlat conversation
+    /// @param x the location to evaluate.
+    /// @return A pair (latitude, longitude)
+    inline latlon to_latlon(const lonlat &x) const noexcept
+    {
+        return latlon(x.lat, x.lon);
+    }
+    
     /// @brief Read the NA in the first band
-    value_type nodata_first_band() const
+    inline value_type nodata_first_band() const
     {
         return _dataset.get().GetRasterBand(1)->GetNoDataValue();
     }
@@ -588,7 +637,7 @@ template <typename Time = int> class raster
         return answer;
     }
 
-    /// @brief checks if there is a 1 to 1 mapping from times to layer
+    /// @brief Check if there is a 1 to 1 mapping from times to layer
     void check_times_equals_depth()
     {
         if (_times.size() != depth())
diff --git a/template/main.cpp b/template/main.cpp
index 827c1c55..1bbb5a95 100644
--- a/template/main.cpp
+++ b/template/main.cpp
@@ -15,7 +15,7 @@ int main()
     std::iota(times.begin(), times.end(), 2001);
 
     // Initialize the landscape: for each spatial variable a string key and a file value, for all a time series.
-    auto landscape = quetzal::geography::landscape<>::from_files({{"suit", file1}, {"DEM", file2}}, times);
+    auto landscape = quetzal::geography::landscape<>::from_file({{"suit", file1}, {"DEM", file2}}, times);
     std::cout << landscape << std::endl;
 
     // Declares some type aliases to shorten notation
diff --git a/test/unit_test/geography_graph_test.cpp b/test/unit_test/geography_graph_test.cpp
index 154c9635..a44a767f 100644
--- a/test/unit_test/geography_graph_test.cpp
+++ b/test/unit_test/geography_graph_test.cpp
@@ -29,7 +29,7 @@ BOOST_AUTO_TEST_CASE(expected_number_edges)
 
     // Initialize the landscape: for each var a key and a file, for all a time series.
     using landscape_type = geo::landscape<>;
-    auto land = landscape_type::from_files({{"bio1", file1}, {"bio12", file2}}, times);
+    auto land = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}}, times);
 
     // Our graph will not store any useful information
     using vertex_info = geo::no_property;
diff --git a/test/unit_test/geography_test.cpp b/test/unit_test/geography_test.cpp
index a64af6cf..994a4812 100644
--- a/test/unit_test/geography_test.cpp
+++ b/test/unit_test/geography_test.cpp
@@ -18,47 +18,66 @@ namespace utf = boost::unit_test;
 
 BOOST_AUTO_TEST_SUITE(geography)
 
-BOOST_AUTO_TEST_CASE(landscape)
+BOOST_AUTO_TEST_CASE(gdalcppTest)
 {
-    using time_type = int;
-    using landscape_type = quetzal::geography::landscape<std::string, time_type>;
-
-    auto file1 = std::filesystem::current_path() / "data/bio1.tif";
-    auto file2 = std::filesystem::current_path() / "data/bio12.tif";
-
-    try
-    {
-
-        auto env = landscape_type::from_files({{"bio1", file1}, {"bio12", file2}},
-                                              {2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010});
-
-        BOOST_CHECK_EQUAL(env.num_variables(), 2);
+    auto file = std::filesystem::current_path() / "data/bio1.tif";
 
-        landscape_type::latlon Bordeaux(44.5, 0.34);
+    quetzal::geography::gdalcpp::Dataset data(file.string());
 
-        BOOST_TEST(env.contains(Bordeaux));
-        BOOST_TEST(env.contains(env.to_centroid(Bordeaux)));
+    std::cout << "Name:\t" << data.dataset_name() << std::endl;
+    std::cout << "Driver:\t " << data.driver_name() << std::endl;
+    std::cout << "Width:\t" << data.width() << std::endl;
+    std::cout << "Height:\t" << data.height() << std::endl;
+    std::cout << "Depth:\t" << data.depth() << std::endl;
+    std::cout << "Coefs:\t";
 
-        const auto &f = env["bio1"].to_view();
-        const auto &g = env["bio12"].to_view();
+    auto v = data.affine_transformation_coefficients();
+}
 
-        BOOST_CHECK_EQUAL(env.times().size(), 10);
-        BOOST_TEST(env.locations().size() < 100);
+BOOST_AUTO_TEST_CASE(resolution)
+{
+    using decimal_degree = double;
+    decimal_degree lat_resolution = 2.0;
+    decimal_degree lon_resolution = 1.0;
+    quetzal::geography::resolution<decimal_degree> res(lat_resolution, lon_resolution);
+    BOOST_TEST(res.lat() == lat_resolution);
+    BOOST_TEST(res.lon() == lon_resolution);
+    auto copy = res;
+    copy.lon(2.0);
+    assert(copy != res);
+}
 
-        for (auto t : env.times())
-        {
-            for (auto x : env.locations())
-            {
-                f(x, t).has_value();
-                g(x, t).has_value();
-            }
-        }
-    }
+BOOST_AUTO_TEST_CASE(coordinates)
+{
+    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;
+    BOOST_TEST(diff < EPSILON);
+    BOOST_TEST(-diff < EPSILON);
+}
 
-    catch (const std::exception &ex)
-    {
-        BOOST_ERROR(ex.what());
-    }
+BOOST_AUTO_TEST_CASE(extent)
+{
+    using decimal_degree = double;
+    using extent_type = quetzal::geography::extent<decimal_degree>;
+    decimal_degree paris_lat = 48.8566;
+    decimal_degree paris_lon = 2.3522;
+    decimal_degree moscow_lat = 55.7522;
+    decimal_degree moscow_lon = 37.6155;
+    extent_type extent(paris_lat, moscow_lat, paris_lon, moscow_lon);
+    BOOST_CHECK_EQUAL(extent.lat_min(), paris_lat);
+    BOOST_CHECK_EQUAL(extent.lat_max(), moscow_lat);
+    BOOST_CHECK_EQUAL(extent.lon_min(), paris_lon);
+    BOOST_CHECK_EQUAL(extent.lon_max(), moscow_lon);
+    auto copy = extent;
+    copy.lat_min(0.0);
+    assert(copy != extent);
 }
 
 BOOST_AUTO_TEST_CASE(raster)
@@ -76,7 +95,9 @@ BOOST_AUTO_TEST_CASE(raster)
     BOOST_CHECK_EQUAL(bio1.origin(), raster_type::latlon(52., -5.));
 
     auto f = bio1.to_view();
-    BOOST_TEST(f(bio1.to_descriptor(bio1.origin()), bio1.times().front()).has_value());
+    auto x = bio1.to_descriptor( bio1.origin() );
+    auto t = bio1.times().front();
+    // BOOST_TEST( f( x, t ).has_value() );
 
     auto space = bio1.locations() |
                  ranges::views::transform([&](auto i) { return bio1.to_latlon(i); }) |
@@ -167,66 +188,47 @@ BOOST_AUTO_TEST_CASE(raster)
     // bio1.to_centroid(E_c2_c5_limit); // assertion raised
 }
 
-BOOST_AUTO_TEST_CASE(extent)
+BOOST_AUTO_TEST_CASE(landscape)
 {
-    using decimal_degree = double;
-    using extent_type = quetzal::geography::extent<decimal_degree>;
-    decimal_degree paris_lat = 48.8566;
-    decimal_degree paris_lon = 2.3522;
-    decimal_degree moscow_lat = 55.7522;
-    decimal_degree moscow_lon = 37.6155;
-    extent_type extent(paris_lat, moscow_lat, paris_lon, moscow_lon);
-    BOOST_CHECK_EQUAL(extent.lat_min(), paris_lat);
-    BOOST_CHECK_EQUAL(extent.lat_max(), moscow_lat);
-    BOOST_CHECK_EQUAL(extent.lon_min(), paris_lon);
-    BOOST_CHECK_EQUAL(extent.lon_max(), moscow_lon);
-    auto copy = extent;
-    copy.lat_min(0.0);
-    assert(copy != extent);
-}
+    using time_type = int;
+    using landscape_type = quetzal::geography::landscape<std::string, time_type>;
 
-BOOST_AUTO_TEST_CASE(coordinates)
-{
-    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;
-    BOOST_TEST(diff < EPSILON);
-    BOOST_TEST(-diff < EPSILON);
-}
+    auto file1 = std::filesystem::current_path() / "data/bio1.tif";
+    auto file2 = std::filesystem::current_path() / "data/bio12.tif";
 
-BOOST_AUTO_TEST_CASE(resolution)
-{
-    using decimal_degree = double;
-    decimal_degree lat_resolution = 2.0;
-    decimal_degree lon_resolution = 1.0;
-    quetzal::geography::resolution<decimal_degree> res(lat_resolution, lon_resolution);
-    BOOST_TEST(res.lat() == lat_resolution);
-    BOOST_TEST(res.lon() == lon_resolution);
-    auto copy = res;
-    copy.lon(2.0);
-    assert(copy != res);
-}
+    try
+    {
 
-BOOST_AUTO_TEST_CASE(gdalcppTest)
-{
-    auto file = std::filesystem::current_path() / "data/bio1.tif";
+        auto env = landscape_type::from_file({{"bio1", file1}, {"bio12", file2}},
+                                              {2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010});
 
-    quetzal::geography::gdalcpp::Dataset data(file.string());
+        BOOST_CHECK_EQUAL(env.num_variables(), 2);
 
-    std::cout << "Name:\t" << data.dataset_name() << std::endl;
-    std::cout << "Driver:\t " << data.driver_name() << std::endl;
-    std::cout << "Width:\t" << data.width() << std::endl;
-    std::cout << "Height:\t" << data.height() << std::endl;
-    std::cout << "Depth:\t" << data.depth() << std::endl;
-    std::cout << "Coefs:\t";
+        landscape_type::latlon Bordeaux(44.5, 0.34);
 
-    auto v = data.affine_transformation_coefficients();
+        BOOST_TEST(env.contains(Bordeaux));
+        BOOST_TEST(env.contains(env.to_centroid(Bordeaux)));
+
+        const auto &f = env["bio1"].to_view();
+        const auto &g = env["bio12"].to_view();
+
+        BOOST_CHECK_EQUAL(env.times().size(), 10);
+        BOOST_TEST(env.locations().size() < 100);
+
+        for (auto t : env.times())
+        {
+            for (auto x : env.locations())
+            {
+                f(x, t).has_value();
+                // g(x, t).has_value();
+            }
+        }
+    }
+
+    catch (const std::exception &ex)
+    {
+        BOOST_ERROR(ex.what());
+    }
 }
 
 BOOST_AUTO_TEST_SUITE_END()