Repository: computationalgeography/lue (language: python)

Issue:
Update LUE to support running the pcr-globwb model
- [x] #659 
- [x] #668 
- [x] #696 
- [x] #705 
- [x] #710 
- [x] #712 


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/linux.yml b/.github/workflows/linux.yml
index ec509bbf3..a1d04fe1c 100644
--- a/.github/workflows/linux.yml
+++ b/.github/workflows/linux.yml
@@ -12,7 +12,7 @@ jobs:
         strategy:
             matrix:
                 os: ["ubuntu-22.04"]
-                python_version: ["3.9", "3.12"]
+                python_version: ["3.10", "3.12"]
                 compiler:
                     # To be recognized by Conan, name must be one of:
                     #     'Visual Studio', 'apple-clang', 'clang', 'gcc', 'intel', 'intel-cc', 'mcst-lcc',
diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml
index 43a868948..daa92180a 100644
--- a/.github/workflows/windows.yml
+++ b/.github/workflows/windows.yml
@@ -16,7 +16,7 @@ jobs:
                   - { name: "visual studio", c: cl, cxx: cl }
                 vs:
                   - { name: "Visual Studio Enterprise 2022", version: "17", msvc_toolset: "14" }
-                python_version: ["3.9"]
+                python_version: ["3.10"]
                 architecture: ["amd64"]
             fail-fast: false
 
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 587ab52fc..03ba518a7 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,7 +1,6 @@
 # See https://pre-commit.com for more information
 # See https://pre-commit.com/hooks.html for more hooks
 exclude: |
-  # Skip the files that must be kept as they are
   (?x)^(
     environment/cmake/.*\.patch|
     environment/cmake/UseLATEX\.cmake|
@@ -11,18 +10,18 @@ exclude: |
   )$
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.6.0
+  rev: v5.0.0
   hooks:
   - id: trailing-whitespace
   - id: end-of-file-fixer
   # - id: check-yaml
   - id: check-added-large-files
 - repo: https://github.com/psf/black
-  rev: 24.8.0
+  rev: 24.10.0
   hooks:
   - id: black
 - repo: https://github.com/pre-commit/mirrors-mypy
-  rev: v1.11.1
+  rev: v1.11.2
   hooks:
   - id: mypy
     # args: [--verbose]
@@ -46,7 +45,7 @@ repos:
         .*\.py.in|
       )$
 - repo: https://github.com/markdownlint/markdownlint
-  rev: v0.13.0
+  rev: v0.12.0
   hooks:
   - id: markdownlint
 - repo: https://github.com/codespell-project/codespell
@@ -56,7 +55,7 @@ repos:
     additional_dependencies:
     - tomli
 - repo: https://github.com/crate-ci/typos
-  rev: v1.23.6
+  rev: v1.26.0
   hooks:
   - id: typos
     args: []  # Relevant: don't update files, just report
diff --git a/environment/cmake/LueConfiguration.cmake b/environment/cmake/LueConfiguration.cmake
index 39c5f2b1c..5822ca70f 100644
--- a/environment/cmake/LueConfiguration.cmake
+++ b/environment/cmake/LueConfiguration.cmake
@@ -250,12 +250,12 @@ endif()
 
 # ------------------------------------------------------------------------------
 if(LUE_PYTHON_REQUIRED)
-    # This is the first numpy version supporting Python 3.9
-    set(LUE_MIN_NUMPY_VERSION 1.19)
+    # This is the first numpy version supporting Python 3.10
+    set(LUE_MIN_NUMPY_VERSION 1.21)
     string(REPLACE "." "_" LUE_NPY_NO_DEPRECATED_API "NPY_${LUE_MIN_NUMPY_VERSION}_API_VERSION")
 
     # Order matters: Pybind11 must be searched for after Python has been found.
-    find_package(Python 3.9 REQUIRED COMPONENTS Interpreter Development NumPy)
+    find_package(Python 3.10 REQUIRED COMPONENTS Interpreter Development NumPy)
 
     if((Python_INTERPRETER_ID STREQUAL "Anaconda") OR (Python_EXECUTABLE MATCHES "^.*conda.*$"))
         set(LUE_PYTHON_FROM_CONDA TRUE)
@@ -542,7 +542,8 @@ if(LUE_DOCOPT_REQUIRED)
         GIT_REPOSITORY https://github.com/docopt/docopt.cpp.git
         GIT_TAG 42ebcec9dc2c99a1b3a4542787572045763ad196  # 0.6.3
         SYSTEM
-        FIND_PACKAGE_ARGS
+        # Using the system version is inconvenient since headers are then in a docopt/ subdirectory...
+        # FIND_PACKAGE_ARGS
     )
     FetchContent_MakeAvailable(docopt)
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/abs.hpp b/source/framework/algorithm/include/lue/framework/algorithm/abs.hpp
index b06ea64a5..622c447c3 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/abs.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/abs.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_signed_v<InputElement> || std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"abs"};
+
                 using OutputElement = InputElement;
 
                 auto operator()(InputElement const input_element) const noexcept -> OutputElement
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/acos.hpp b/source/framework/algorithm/include/lue/framework/algorithm/acos.hpp
index 5c501650f..7ffc0d07e 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/acos.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/acos.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"acos"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/add.hpp b/source/framework/algorithm/include/lue/framework/algorithm/add.hpp
index cc5ddb537..b0ae21b48 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/add.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/add.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"add"};
+
                 using OutputElement = OutputElement_;
 
                 constexpr auto operator()(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/asin.hpp b/source/framework/algorithm/include/lue/framework/algorithm/asin.hpp
index 06dc86f05..2a89310b7 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/asin.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/asin.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"asin"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/atan.hpp b/source/framework/algorithm/include/lue/framework/algorithm/atan.hpp
index 13d17ef40..ee4957247 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/atan.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/atan.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"atan"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/cast.hpp b/source/framework/algorithm/include/lue/framework/algorithm/cast.hpp
index 2d2f22f4d..4b8501c26 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/cast.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/cast.hpp
@@ -12,6 +12,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"cast"};
+
                 using OutputElement = OutputElement_;
 
                 static_assert(std::is_arithmetic_v<InputElement>);
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/ceil.hpp b/source/framework/algorithm/include/lue/framework/algorithm/ceil.hpp
index dc2bea834..7738d7e09 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/ceil.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/ceil.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"ceil"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/close_to.hpp b/source/framework/algorithm/include/lue/framework/algorithm/close_to.hpp
index b8574b41c..c9fcc9e28 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/close_to.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/close_to.hpp
@@ -17,6 +17,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"close_to"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/cos.hpp b/source/framework/algorithm/include/lue/framework/algorithm/cos.hpp
index 8de1f6997..2d7166074 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/cos.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/cos.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"cos"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/create_partitioned_array.hpp b/source/framework/algorithm/include/lue/framework/algorithm/create_partitioned_array.hpp
index 2687328d0..4c50f6685 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/create_partitioned_array.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/create_partitioned_array.hpp
@@ -434,14 +434,14 @@ namespace lue {
 
 
         template<typename Policies, typename Shape, typename Functor>
-        std::tuple<
-            Localities<rank<Shape>>,
-            PartitionsT<PartitionedArray<OutputElementT<Functor>, rank<Shape>>>>
-        instantiate_partitions(
+        auto instantiate_partitions(
             Policies const& policies,
             Shape const& array_shape,
             Shape const& partition_shape,
             Functor const& partition_creator)
+            -> std::tuple<
+                Localities<rank<Shape>>,
+                PartitionsT<PartitionedArray<OutputElementT<Functor>, rank<Shape>>>>
         {
             // Create array containing partitions. Each of these partitions
             // will be a component client instance referring to a, possibly
@@ -528,14 +528,14 @@ namespace lue {
 
 
         template<typename Policies, typename InputPartitions, typename Functor>
-        std::tuple<
-            Localities<rank<InputPartitions>>,
-            PartitionsT<PartitionedArray<OutputElementT<Functor>, rank<InputPartitions>>>>
-        instantiate_partitions(
+        auto instantiate_partitions(
             Policies const& policies,
             ShapeT<InputPartitions> const& array_shape,
             InputPartitions const& input_partitions,
             Functor const& partition_creator)
+            -> std::tuple<
+                Localities<rank<InputPartitions>>,
+                PartitionsT<PartitionedArray<OutputElementT<Functor>, rank<InputPartitions>>>>
         {
             // Create partitions. Each of these partitions will be a component client instance
             // referring to a, possibly remote, component server instance.
@@ -611,19 +611,18 @@ namespace lue {
 
 
             template<typename Policies>
-            Partition instantiate(
+            auto instantiate(
                 hpx::id_type const locality_id,
                 [[maybe_unused]] Policies const& policies,
                 [[maybe_unused]] Shape const& array_shape,
                 [[maybe_unused]] lue::Index const partition_idx,
                 Offset const& offset,
-                Shape const& partition_shape)
+                Shape const& partition_shape) -> Partition
             {
                 return hpx::async(
 
-                    [locality_id, offset, partition_shape]() {
-                        return Partition{locality_id, offset, partition_shape};
-                    }
+                    [locality_id, offset, partition_shape]()
+                    { return Partition{locality_id, offset, partition_shape}; }
 
                 );
             }
@@ -665,9 +664,8 @@ namespace lue {
                 Shape const& partition_shape) const -> Partition
             {
                 return _fill_value.then(
-                    [locality_id, offset, partition_shape](hpx::shared_future<Element> const& fill_value) {
-                        return Partition{locality_id, offset, partition_shape, fill_value.get()};
-                    });
+                    [locality_id, offset, partition_shape](hpx::shared_future<Element> const& fill_value)
+                    { return Partition{locality_id, offset, partition_shape, fill_value.get()}; });
             }
 
 
@@ -740,13 +738,13 @@ namespace lue {
 
 
             template<typename Policies>
-            Partition instantiate(
+            auto instantiate(
                 hpx::id_type const locality_id,
                 Policies const& policies,
                 Shape const& array_shape,
                 [[maybe_unused]] lue::Index const partition_idx,
                 Offset const& offset,
-                Shape const& partition_shape)
+                Shape const& partition_shape) -> Partition
             {
                 using Data = lue::ArrayPartitionData<OutputElement, rank>;
 
@@ -824,36 +822,6 @@ namespace lue {
     };
 
 
-    template<typename Element, lue::Rank rank>
-    class FunctorTraits<InstantiateFilled<Element, rank>>
-    {
-
-        public:
-
-            static constexpr bool is_functor{true};
-    };
-
-
-    template<typename Element, lue::Rank rank>
-    class FunctorTraits<InstantiateDefaultInitialized<Element, rank>>
-    {
-
-        public:
-
-            static constexpr bool is_functor{true};
-    };
-
-
-    template<typename BufferHandle, Rank rank>
-    class FunctorTraits<InstantiateFromBuffer<BufferHandle, rank>>
-    {
-
-        public:
-
-            static constexpr bool is_functor{true};
-    };
-
-
     // TODO, as part of range algorithm
     // template<
     //     typename Element,
@@ -894,6 +862,13 @@ namespace lue {
     }  // namespace policy::create_partitioned_array
 
 
+    template<typename Functor>
+    concept PartitionInstantiator = requires(Functor functor)
+    {
+        functor.instantiate_per_locality;
+    };
+
+
     /*!
         @brief      Create a partitioned array
         @param      policies Operation policies
@@ -912,11 +887,7 @@ namespace lue {
         - Fill all elements with random values.
         - Fill all elements with values read from a file on disk.
     */
-    template<
-        typename Policies,
-        Rank rank,
-        typename Functor,
-        std::enable_if_t<is_functor_v<Functor>>* = nullptr>
+    template<typename Policies, Rank rank, PartitionInstantiator Functor>
     auto create_partitioned_array(
         Policies const& policies,
         Shape<Count, rank> const& array_shape,
@@ -939,15 +910,11 @@ namespace lue {
     }
 
 
-    template<
-        typename Policies,
-        typename RouteID,
-        Rank rank,
-        typename Functor,
-        std::enable_if_t<is_functor_v<Functor>>* = nullptr>
+    template<typename Policies, typename RouteID, Rank rank, PartitionInstantiator Functor>
     auto create_partitioned_array(
-        Policies const& policies, SerialRoute<RouteID, rank> const& route, Functor const& partition_creator)
-        -> PartitionedArray<OutputElementT<Functor>, rank>
+        Policies const& policies,
+        SerialRoute<RouteID, rank> const& route,
+        Functor const& partition_creator) -> PartitionedArray<OutputElementT<Functor>, rank>
     {
         // Use the route passed in as a clone / template for creating the new array. The shape
         // of the output array and the individual partitions, and the distribution of the
@@ -968,21 +935,18 @@ namespace lue {
 
         The partition shape is based on @a default_partition_shape.
     */
-    template<
-        typename Policies,
-        Rank rank,
-        typename Functor,
-        std::enable_if_t<is_functor_v<Functor>>* = nullptr>
+    template<typename Policies, Rank rank, PartitionInstantiator Functor>
     auto create_partitioned_array(
-        Policies const& policies, Shape<Count, rank> const& array_shape, Functor const& partition_creator)
-        -> PartitionedArray<OutputElementT<Functor>, rank>
+        Policies const& policies,
+        Shape<Count, rank> const& array_shape,
+        Functor const& partition_creator) -> PartitionedArray<OutputElementT<Functor>, rank>
     {
         return create_partitioned_array<Policies, rank, Functor>(
             policies, array_shape, default_partition_shape(array_shape), partition_creator);
     }
 
 
-    template<Rank rank, typename Functor, std::enable_if_t<is_functor_v<Functor>>* = nullptr>
+    template<Rank rank, PartitionInstantiator Functor>
     auto create_partitioned_array(
         Shape<Count, rank> const& array_shape,
         Shape<Count, rank> const& partition_shape,
@@ -996,7 +960,7 @@ namespace lue {
     }
 
 
-    template<Rank rank, typename Functor, std::enable_if_t<is_functor_v<Functor>>* = nullptr>
+    template<Rank rank, PartitionInstantiator Functor>
     auto create_partitioned_array(Shape<Count, rank> const& array_shape, Functor const& partition_creator)
         -> PartitionedArray<OutputElementT<Functor>, rank>
     {
@@ -1005,11 +969,7 @@ namespace lue {
     }
 
 
-    template<
-        typename RouteID,
-        Rank rank,
-        typename Functor,
-        std::enable_if_t<is_functor_v<Functor>>* = nullptr>
+    template<typename RouteID, Rank rank, PartitionInstantiator Functor>
     auto create_partitioned_array(SerialRoute<RouteID, rank> const& route, Functor const& partition_creator)
         -> PartitionedArray<OutputElementT<Functor>, rank>
     {
@@ -1023,8 +983,8 @@ namespace lue {
 
     template<typename Element, Rank rank, std::enable_if_t<std::is_arithmetic_v<Element>>* = nullptr>
     auto create_partitioned_array(
-        Shape<Count, rank> const& array_shape, Shape<Count, rank> const& partition_shape)
-        -> PartitionedArray<Element, rank>
+        Shape<Count, rank> const& array_shape,
+        Shape<Count, rank> const& partition_shape) -> PartitionedArray<Element, rank>
     {
         using Policies = policy::create_partitioned_array::DefaultPolicies<Element>;
         using Functor = InstantiateDefaultInitialized<Element, rank>;
@@ -1069,8 +1029,8 @@ namespace lue {
 
     template<typename Element, Rank rank>
     requires std::is_arithmetic_v<Element> auto create_partitioned_array(
-        Shape<Count, rank> const& array_shape, Scalar<Element> const& fill_value)
-        -> PartitionedArray<Element, rank>
+        Shape<Count, rank> const& array_shape,
+        Scalar<Element> const& fill_value) -> PartitionedArray<Element, rank>
     {
         return create_partitioned_array(array_shape, default_partition_shape(array_shape), fill_value);
     }
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/accu_threshold3.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/accu_threshold3.hpp
index 499627d92..dc26a0d15 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/accu_threshold3.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/accu_threshold3.hpp
@@ -239,7 +239,8 @@ namespace lue {
                             hpx::future<std::array<CellsIdxs, nr_neighbours<rank>()>>&&
                                 output_cells_idxs_f) mutable
                         {
-                            AnnotateFunction annotation{"intra_partition_stream"};
+                            AnnotateFunction const annotation{
+                                "accu_threshold: partition: intra_partition_stream"};
                             auto const flow_direction_partition_ptr{
                                 ready_component_ptr(flow_direction_partition)};
                             FlowDirectionData const& flow_direction_data{
@@ -343,7 +344,8 @@ namespace lue {
                         hpx::future<MaterialData>&& outflow_data_f,
                         hpx::future<MaterialData>&& remainder_data_f) mutable
                     {
-                        AnnotateFunction annotation{"inter_partition_stream"};
+                        AnnotateFunction const annotation{
+                            "accu_threshold: partition: inter_partition_stream"};
 
                         std::vector<accu::Direction> const directions{
                             accu::Direction::north,
@@ -626,7 +628,7 @@ namespace lue {
                    hpx::future<MaterialData>&& outflow_data_f,
                    hpx::future<MaterialData>&& remainder_data_f)
                 {
-                    AnnotateFunction annotation{"create_result_partitions"};
+                    AnnotateFunction const annotation{"accu_threshold: partition: create_result_partitions"};
                     using Server = typename MaterialPartition::Server;
 
                     Offset const partition_offset{ready_component_ptr(flow_direction_partition)->offset()};
@@ -668,7 +670,7 @@ namespace lue {
         PartitionedArray<MaterialElement, rank> const& external_inflow,
         PartitionedArray<MaterialElement, rank> const& threshold)
     {
-        AnnotateFunction annotation{"accu_threshold_meh"};
+        AnnotateFunction const annotation{"accu_threshold: array"};
 
         detail::verify_compatible(flow_direction, external_inflow, threshold);
 
@@ -707,7 +709,8 @@ namespace lue {
         for (Index p = 0; p < nr_partitions; ++p)
         {
             hpx::tie(outflow_partitions[p], remainder_partitions[p]) = hpx::split_future(hpx::async(
-                hpx::annotated_function(action, "accu_threshold"),
+                // hpx::annotated_function(action, "accu_threshold"),
+                action,
                 localities[p],
                 policies,
                 flow_direction.partitions()[p],
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/aspect.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/aspect.hpp
index 6700381fd..95f550685 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/aspect.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/aspect.hpp
@@ -17,12 +17,14 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"aspect"};
+
                 using OutputElement = InputElement;
 
                 OutputElement operator()(
                     InputElement const gradient_x, InputElement const gradient_y) const noexcept
                 {
-                    OutputElement result;
+                    OutputElement result{};
 
                     if (gradient_x == 0)
                     {
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/atan2.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/atan2.hpp
index 532a68a05..ae54c52f6 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/atan2.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/atan2.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"atan2"};
+
                 using OutputElement = InputElement;
 
                 OutputElement operator()(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/binary_local_operation.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/binary_local_operation.hpp
index c8d0ec426..b372a6fc3 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/binary_local_operation.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/binary_local_operation.hpp
@@ -1,6 +1,7 @@
 #pragma once
 #include "lue/framework/algorithm/binary_local_operation.hpp"
 #include "lue/framework/algorithm/detail/verify_compatible.hpp"
+#include "lue/framework/algorithm/functor_traits.hpp"
 #include "lue/framework/algorithm/local_operation_export.hpp"
 #include "lue/framework/core/annotate.hpp"
 #include "lue/framework/core/component.hpp"
@@ -67,7 +68,8 @@ namespace lue {
                                     InputData1 const& input_partition_data1,
                                     InputData2 const& input_partition_data2)
                                 {
-                                    AnnotateFunction annotation{"binary_local_operation_partition"};
+                                    AnnotateFunction const annotation{
+                                        fmt::format("{}: partition", functor_name<Functor>)};
 
                                     HPX_UNUSED(input_partition1);
                                     HPX_UNUSED(input_partition2);
@@ -193,7 +195,8 @@ namespace lue {
                                 [policies, input_partition, functor, input_scalar](
                                     Offset const& offset, InputData const& input_partition_data)
                                 {
-                                    AnnotateFunction annotation{"binary_local_operation_partition"};
+                                    AnnotateFunction const annotation{
+                                        fmt::format("{}: partition", functor_name<Functor>)};
 
                                     HPX_UNUSED(input_partition);
 
@@ -324,7 +327,8 @@ namespace lue {
                                 [policies, input_partition, functor, input_scalar](
                                     Offset const& offset, InputData const& input_partition_data)
                                 {
-                                    AnnotateFunction annotation{"binary_local_operation_partition"};
+                                    AnnotateFunction const annotation{
+                                        fmt::format("{}: partition", functor_name<Functor>)};
 
                                     HPX_UNUSED(input_partition);
 
@@ -438,6 +442,8 @@ namespace lue {
         using InputPartition = ArrayPartition<policy::InputElementT<Policies, 0>, rank>;
         using OutputPartition = ArrayPartition<policy::OutputElementT<Policies, 0>, rank>;
 
+        AnnotateFunction const annotation{fmt::format("{}: partition", functor_name<Functor>)};
+
         lue_hpx_assert(input_partition.valid());
         lue_hpx_assert(input_scalar.valid());
 
@@ -449,13 +455,7 @@ namespace lue {
             Functor>
             action;
 
-        return hpx::async(
-            hpx::annotated_function(action, "binary_local_operation"),
-            locality_id,
-            policies,
-            input_partition,
-            input_scalar,
-            functor);
+        return hpx::async(action, locality_id, policies, input_partition, input_scalar, functor);
     }
 
 
@@ -482,6 +482,8 @@ namespace lue {
 
         using Shape = ShapeT<OutputArray>;
 
+        AnnotateFunction const annotation{fmt::format("{}: array", functor_name<Functor>)};
+
         detail::verify_compatible(input_array1, input_array2);
 
         lue_hpx_assert(all_are_valid(input_array1.partitions()));
@@ -507,12 +509,7 @@ namespace lue {
         for (Index p = 0; p < nr_partitions; ++p)
         {
             output_partitions[p] = hpx::async(
-                hpx::annotated_function(action, "binary_local_operation"),
-                localities[p],
-                policies,
-                input_partitions1[p],
-                input_partitions2[p],
-                functor);
+                action, localities[p], policies, input_partitions1[p], input_partitions2[p], functor);
         }
 
         return OutputArray{shape(input_array1), localities, std::move(output_partitions)};
@@ -537,6 +534,8 @@ namespace lue {
 
         using Shape = ShapeT<OutputArray>;
 
+        AnnotateFunction const annotation{fmt::format("{}: array", functor_name<Functor>)};
+
         lue_hpx_assert(all_are_valid(input_array.partitions()));
         lue_hpx_assert(input_scalar.valid());
 
@@ -558,13 +557,8 @@ namespace lue {
 
         for (Index p = 0; p < nr_partitions; ++p)
         {
-            output_partitions[p] = hpx::async(
-                hpx::annotated_function(action, "binary_local_operation"),
-                localities[p],
-                policies,
-                input_partitions[p],
-                input_scalar,
-                functor);
+            output_partitions[p] =
+                hpx::async(action, localities[p], policies, input_partitions[p], input_scalar, functor);
         }
 
         return OutputArray{shape(input_array), localities, std::move(output_partitions)};
@@ -589,6 +583,8 @@ namespace lue {
 
         using Shape = ShapeT<OutputArray>;
 
+        AnnotateFunction const annotation{fmt::format("{}: array", functor_name<Functor>)};
+
         lue_hpx_assert(input_scalar.valid());
         lue_hpx_assert(all_are_valid(input_array.partitions()));
 
@@ -610,13 +606,8 @@ namespace lue {
 
         for (Index p = 0; p < nr_partitions; ++p)
         {
-            output_partitions[p] = hpx::async(
-                hpx::annotated_function(action, "binary_local_operation"),
-                localities[p],
-                policies,
-                input_scalar,
-                input_partitions[p],
-                functor);
+            output_partitions[p] =
+                hpx::async(action, localities[p], policies, input_scalar, input_partitions[p], functor);
         }
 
         return OutputArray{shape(input_array), localities, std::move(output_partitions)};
@@ -633,6 +624,8 @@ namespace lue {
     {
         using OutputElement = policy::OutputElementT<Policies, 0>;
 
+        AnnotateFunction const annotation{fmt::format("{}: scalar", functor_name<Functor>)};
+
         return hpx::dataflow(
             hpx::launch::async,
             hpx::unwrapping(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/clump.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/clump.hpp
index d0650fe49..2c315635a 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/clump.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/clump.hpp
@@ -825,7 +825,7 @@ namespace lue {
                     [policies, zone_partition, connectivity](
                         auto const& offset, auto const& zone_data) -> LocalResult
                     {
-                        AnnotateFunction const annotate{"clump_partition"};
+                        AnnotateFunction const annotate{"clump: partition"};
 
                         HPX_UNUSED(zone_partition);
 
@@ -1220,7 +1220,7 @@ namespace lue {
                     [policies, shape_in_partitions = local_result_fs.shape(), connectivity](
                         auto&& local_result_fs) -> LocalResults<Partition>
                     {
-                        AnnotateFunction const annotate{"determine_reclass_tables"};
+                        AnnotateFunction const annotation{"clump: array: determine_reclass_tables"};
                         auto const& indp = std::get<0>(policies.inputs_policies()).input_no_data_policy();
                         auto const& ondp = std::get<0>(policies.outputs_policies()).output_no_data_policy();
 
@@ -1835,7 +1835,7 @@ namespace lue {
         auto reclass_partition(NoDataPolicy const& ndp, Partition const& clump_partition, LUT const& lut)
             -> Partition
         {
-            AnnotateFunction const annotate{"reclass_partition"};
+            AnnotateFunction const annotation{"clump: partition: reclass"};
 
             lue_hpx_assert(clump_partition.is_ready());
 
@@ -1934,7 +1934,7 @@ namespace lue {
             std::is_same_v<policy::InputElementT<Policies, 0>, policy::OutputElementT<Policies, 0>>);
         static_assert(rank == 2);
 
-        AnnotateFunction const annotate{"clump"};
+        AnnotateFunction const annotation{"clump: array"};
 
         using ClumpArray = PartitionedArray<policy::OutputElementT<Policies, 0>, rank>;
         using ClumpPartitions = PartitionsT<ClumpArray>;
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/convolve.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/convolve.hpp
index ae6bf9275..891a61f6d 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/convolve.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/convolve.hpp
@@ -17,6 +17,8 @@ namespace lue {
 
                 static_assert(std::is_convertible_v<InputElement, OutputElement>);
 
+                static constexpr char const* name{"convolve"};
+
 
                 Convolve() = default;
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/d8_flow_direction.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/d8_flow_direction.hpp
index eac2b2775..b168439bd 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/d8_flow_direction.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/d8_flow_direction.hpp
@@ -16,6 +16,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"d8_flow_direction"};
+
                 using OutputElement = FlowDirectionElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream.hpp
index 3bc3b6e9b..8d773425e 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"downstream"};
+
                 using OutputElement = MaterialElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream_distance.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream_distance.hpp
index da92a5ac4..61834a5ae 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream_distance.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/downstream_distance.hpp
@@ -19,6 +19,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<DistanceElement>);
 
+                static constexpr char const* name{"downstream_distance"};
+
                 using OutputElement = DistanceElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_diversity.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_diversity.hpp
index 5d50f8610..20e7011fb 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_diversity.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_diversity.hpp
@@ -21,6 +21,8 @@ namespace lue {
                 static_assert(std::is_integral_v<InputElement>);
                 static_assert(std::is_integral_v<OutputElement>);
 
+                static constexpr char const* name{"focal_diversity"};
+
 
                 template<typename Kernel, typename OutputPolicies, typename InputPolicies, typename Subspan>
                 Count operator()(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_high_pass.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_high_pass.hpp
index a9570a762..6bfb4f9b5 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_high_pass.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_high_pass.hpp
@@ -16,6 +16,8 @@ namespace lue {
                 // If not, an OutputElement type must be passed in
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"focal_high_pass"};
+
                 using OutputElement = InputElement;
 
 
@@ -40,7 +42,7 @@ namespace lue {
                     auto const& indp = input_policies.input_no_data_policy();
                     auto const& ondp = output_policies.output_no_data_policy();
 
-                    OutputElement focal_value;
+                    OutputElement focal_value{};
                     OutputElement sum_of_neighbours{0};
                     Count nr_neighbours{0};
 
@@ -79,7 +81,7 @@ namespace lue {
                         }
                     }
 
-                    OutputElement high_pass;
+                    OutputElement high_pass{};
 
                     {
                         if (nr_neighbours == 0)
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_majority.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_majority.hpp
index ead06c3f6..94c938540 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_majority.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_majority.hpp
@@ -16,6 +16,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<InputElement>);
 
+                static constexpr char const* name{"focal_majority"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_maximum.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_maximum.hpp
index 54ae748db..6b71f4191 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_maximum.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_maximum.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"focal_maximum"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_mean.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_mean.hpp
index e41f06435..01914070c 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_mean.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_mean.hpp
@@ -16,6 +16,8 @@ namespace lue {
                 // If not, an OutputElement type must be passed in
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"focal_mean"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_minimum.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_minimum.hpp
index 2a670a55b..ec594779a 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_minimum.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_minimum.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"focal_minimum"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_operation.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_operation.hpp
index 5a9733ce7..18844ae94 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_operation.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_operation.hpp
@@ -3,6 +3,7 @@
 #include "lue/framework/algorithm/detail/verify_compatible.hpp"
 #include "lue/framework/algorithm/detail/when_all_get.hpp"
 #include "lue/framework/algorithm/export.hpp"
+#include "lue/framework/algorithm/functor_traits.hpp"
 #include "lue/framework/algorithm/policy.hpp"
 #include "lue/framework/core/annotate.hpp"
 #include "lue/framework/core/array.hpp"
@@ -851,7 +852,8 @@ namespace lue {
                         [policies, input_partitions = focal_input_partitions, kernel, functor](
                             meh::InputData<InputPartitions> const&... partition_data)
                         {
-                            AnnotateFunction annotation{"focal_operation_partition"};
+                            AnnotateFunction const annotation{
+                                fmt::format("{}: partition", functor_name<Functor>)};
 
                             HPX_UNUSED(input_partitions);
 
@@ -912,7 +914,8 @@ namespace lue {
                                 Array<meh::InputData<InputPartitions>, rank<Kernel>> const&... partition_data,
                             OutputData output_partition_data)
                         {
-                            AnnotateFunction annotation{"focal_operation_partition"};
+                            AnnotateFunction const annotation{
+                                fmt::format("{}: partition", functor_name<Functor>)};
 
                             HPX_UNUSED(input_partitions);
 
@@ -1550,13 +1553,13 @@ namespace lue {
             typename Kernel,
             typename Functor,
             typename... InputPartitions>
-        ArrayPartition<OutputElementT<Functor>, rank<Kernel>> spawn_focal_operation_partition(
+        auto spawn_focal_operation_partition(
             hpx::id_type const locality_id,
             Action const& action,
             Policies const& policies,
             Kernel const& kernel,
             Functor const& functor,
-            InputPartitions&&... input_partitions)
+            InputPartitions&&... input_partitions) -> ArrayPartition<OutputElementT<Functor>, rank<Kernel>>
         {
             static_assert(rank<Kernel> == 2);
 
@@ -1565,8 +1568,6 @@ namespace lue {
                     [locality_id, action, policies, kernel, functor](
                         hpx::tuple<hpx::future<InputPartitions>...>&& input_partitions)
                     {
-                        AnnotateFunction annotation{"focal_operation"};
-
                         auto call_action = [locality_id, action, policies, kernel, functor](
                                                InputPartitions&&... input_partitions) {
                             return action(
@@ -1972,25 +1973,28 @@ namespace lue {
 
 
     template<typename Policies, typename... InputElements, Rank rank, typename Kernel, typename Functor>
-    PartitionedArray<detail::OutputElementT<Functor>, rank> focal_operation(
+    auto focal_operation(
         Policies const& policies,
         Kernel const& kernel,
         Functor functor,
         PartitionedArray<InputElements, rank> const&... arrays)
+        -> PartitionedArray<detail::OutputElementT<Functor>, rank>
     {
         static_assert(sizeof...(arrays) >= 1);
         static_assert(rank == 2);
 
+        AnnotateFunction const annotation{fmt::format("{}: array", functor_name<Functor>)};
+
         return detail::focal_operation_2d(policies, kernel, std::move(functor), arrays...);
     }
 
 
     template<typename Policies, typename InputElement, Rank rank, typename Kernel, typename Functor>
-    PartitionedArray<detail::OutputElementT<Functor>, rank> focal_operation(
+    auto focal_operation(
         Policies const& policies,
         PartitionedArray<InputElement, rank> const& array,
         Kernel const& kernel,
-        Functor functor)
+        Functor functor) -> PartitionedArray<detail::OutputElementT<Functor>, rank>
     {
         return focal_operation(policies, kernel, std::move(functor), array);
     }
@@ -2003,12 +2007,12 @@ namespace lue {
         Rank rank,
         typename Kernel,
         typename Functor>
-    PartitionedArray<detail::OutputElementT<Functor>, rank> focal_operation(
+    auto focal_operation(
         Policies const& policies,
         PartitionedArray<InputElement1, rank> const& array1,
         PartitionedArray<InputElement2, rank> const& array2,
         Kernel const& kernel,
-        Functor functor)
+        Functor functor) -> PartitionedArray<detail::OutputElementT<Functor>, rank>
     {
         return focal_operation(policies, kernel, std::move(functor), array1, array2);
     }
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_sum.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_sum.hpp
index e2dd79d41..62972b2af 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_sum.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/focal_sum.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"focal_sum"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/unary_local_operation.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/unary_local_operation.hpp
index 0fc91bf58..dccc5b0a4 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/unary_local_operation.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/unary_local_operation.hpp
@@ -1,4 +1,5 @@
 #pragma once
+#include "lue/framework/algorithm/functor_traits.hpp"
 #include "lue/framework/algorithm/local_operation_export.hpp"
 #include "lue/framework/algorithm/unary_local_operation.hpp"
 #include "lue/framework/core/annotate.hpp"
@@ -27,7 +28,8 @@ namespace lue {
                     [input_partition, policies, functor](
                         Offset const& offset, InputData const& input_partition_data)
                     {
-                        AnnotateFunction annotation{"unary_local_operation_partition"};
+                        AnnotateFunction const annotation{
+                            fmt::format("{}: partition", functor_name<Functor>)};
 
                         HPX_UNUSED(input_partition);
 
@@ -107,6 +109,8 @@ namespace lue {
         using OutputPartitions = PartitionsT<OutputArray>;
         using OutputPartition = PartitionT<OutputArray>;
 
+        AnnotateFunction const annotation{fmt::format("{}: array", functor_name<Functor>)};
+
         lue_hpx_assert(all_are_valid(input_array.partitions()));
 
         detail::UnaryLocalOperationPartitionAction<Policies, InputPartition, OutputPartition, Functor> action;
@@ -122,11 +126,7 @@ namespace lue {
 
                 [locality_id = localities[partition_idx], action, policies, functor](
                     InputPartition const& input_partition)
-                {
-                    AnnotateFunction annotation{"unary_local_operation"};
-
-                    return action(locality_id, policies, input_partition, functor);
-                },
+                { return action(locality_id, policies, input_partition, functor); },
 
                 input_partitions[partition_idx]);
         }
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/upstream.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/upstream.hpp
index 97712116d..3eb5b2930 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/upstream.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/upstream.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"upstream"};
+
                 using OutputElement = MaterialElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/where.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/where.hpp
index 8769a780c..f2c976513 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/where.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/where.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"where"};
+
                 using OutputElement = InputElement;
 
                 constexpr auto operator()(
@@ -77,7 +79,7 @@ namespace lue {
                                 Data const& true_data,
                                 Data const& false_data)
                             {
-                                AnnotateFunction annotation{"where"};
+                                AnnotateFunction const annotation{"where: partition"};
 
                                 HPX_UNUSED(condition_partition);
                                 HPX_UNUSED(true_partition);
@@ -194,7 +196,7 @@ namespace lue {
                                 ConditionData const& condition_data,
                                 Data const& true_data)
                             {
-                                AnnotateFunction annotation{"where_partition"};
+                                AnnotateFunction const annotation{"where: partition"};
 
                                 HPX_UNUSED(condition_partition);
                                 HPX_UNUSED(true_partition);
@@ -307,7 +309,7 @@ namespace lue {
                                 ConditionData const& condition_data,
                                 Data const& false_data)
                             {
-                                AnnotateFunction annotation{"where_partition"};
+                                AnnotateFunction const annotation{"where: partition"};
 
                                 HPX_UNUSED(condition_partition);
                                 HPX_UNUSED(false_partition);
@@ -413,7 +415,7 @@ namespace lue {
                             [policies, condition_partition, true_scalar, false_scalar](
                                 Offset const& offset, ConditionData const& condition_data)
                             {
-                                AnnotateFunction annotation{"where_partition"};
+                                AnnotateFunction const annotation{"where: partition"};
 
                                 HPX_UNUSED(condition_partition);
 
@@ -515,6 +517,8 @@ namespace lue {
         using Partitions = PartitionsT<Array>;
         using Partition = PartitionT<Array>;
 
+        AnnotateFunction const annotation{"where: array"};
+
         detail::verify_compatible(condition, true_array, false_array);
 
         detail::where::WherePartitionAction<Policies, ConditionPartition, Partition, Partition> action;
@@ -533,10 +537,7 @@ namespace lue {
                 [locality_id = localities[p], action, policies](
                     ConditionPartition const& condition_partition,
                     Partition const& true_partition,
-                    Partition const& false_partition)
-                {
-                    AnnotateFunction annotation{"where"};
-
+                    Partition const& false_partition) {
                     return action(
                         locality_id, policies, condition_partition, true_partition, false_partition);
                 },
@@ -566,6 +567,8 @@ namespace lue {
         using Partitions = PartitionsT<Array>;
         using Partition = PartitionT<Array>;
 
+        AnnotateFunction const annotation{"where: array"};
+
         detail::verify_compatible(condition, true_array);
 
         detail::where::WherePartitionAction<Policies, ConditionPartition, Partition, Element> action;
@@ -583,10 +586,7 @@ namespace lue {
                 [locality_id = localities[p], action, policies](
                     ConditionPartition const& condition_partition,
                     Partition const& true_partition,
-                    hpx::shared_future<Element> const& false_value_f)
-                {
-                    AnnotateFunction annotation{"where"};
-
+                    hpx::shared_future<Element> const& false_value_f) {
                     return action(
                         locality_id, policies, condition_partition, true_partition, false_value_f.get());
                 },
@@ -616,6 +616,8 @@ namespace lue {
         using Partitions = PartitionsT<Array>;
         using Partition = PartitionT<Array>;
 
+        AnnotateFunction const annotation{"where: array"};
+
         detail::verify_compatible(condition, false_array);
 
         detail::where::WherePartitionAction<Policies, ConditionPartition, Element, Partition> action;
@@ -633,10 +635,7 @@ namespace lue {
                 [locality_id = localities[p], action, policies](
                     ConditionPartition const& condition_partition,
                     hpx::shared_future<Element> const& true_value_f,
-                    Partition const& false_partition)
-                {
-                    AnnotateFunction annotation{"where"};
-
+                    Partition const& false_partition) {
                     return action(
                         locality_id, policies, condition_partition, true_value_f.get(), false_partition);
                 },
@@ -665,6 +664,8 @@ namespace lue {
         using Array = PartitionedArray<Element, rank>;
         using Partitions = PartitionsT<Array>;
 
+        AnnotateFunction const annotation{"where: array"};
+
         detail::where::WherePartitionAction<Policies, ConditionPartition, Element, Element> action;
 
         Localities<rank> const& localities{condition.localities()};
@@ -679,10 +680,7 @@ namespace lue {
                 [locality_id = localities[p], action, policies](
                     ConditionPartition const& condition_partition,
                     hpx::shared_future<Element> const& true_value_f,
-                    hpx::shared_future<Element> const& false_value_f)
-                {
-                    AnnotateFunction annotation{"where"};
-
+                    hpx::shared_future<Element> const& false_value_f) {
                     return action(
                         locality_id, policies, condition_partition, true_value_f.get(), false_value_f.get());
                 },
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_diversity.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_diversity.hpp
index d3e4dc187..7fdbf2633 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_diversity.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_diversity.hpp
@@ -18,6 +18,8 @@ namespace lue {
                 static_assert(std::is_integral_v<InputElement>);
                 static_assert(std::is_integral_v<Zone>);
 
+                static constexpr char const* name{"zonal_diversity"};
+
                 using OutputElement = Count;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_majority.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_majority.hpp
index e07367aec..ce3593a8a 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_majority.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_majority.hpp
@@ -17,6 +17,8 @@ namespace lue {
                 static_assert(std::is_integral_v<InputElement>);
                 static_assert(std::is_integral_v<Zone>);
 
+                static constexpr char const* name{"zonal_majority"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_maximum.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_maximum.hpp
index 5e2bded13..885cf03f8 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_maximum.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_maximum.hpp
@@ -18,6 +18,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<Zone>);
 
+                static constexpr char const* name{"zonal_maximum"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_mean.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_mean.hpp
index 3d6998b9a..1717ec4df 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_mean.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_mean.hpp
@@ -18,6 +18,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<Zone>);
 
+                static constexpr char const* name{"zonal_mean"};
+
                 using OutputElement = InputElement;
 
                 using Sum = InputElement;
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_minimum.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_minimum.hpp
index 10c85c990..342f17f57 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_minimum.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_minimum.hpp
@@ -18,6 +18,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<Zone>);
 
+                static constexpr char const* name{"zonal_minimum"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_operation.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_operation.hpp
index c8e07d1d1..325a27e03 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_operation.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_operation.hpp
@@ -43,7 +43,8 @@ namespace lue {
                                 [policies, input_element, zones_partition](
                                     ZonesData const& zones_partition_data)
                                 {
-                                    AnnotateFunction annotation{"zonal_operation_partition"};
+                                    AnnotateFunction const annotation{
+                                        fmt::format("{}: partition", functor_name<Functor>)};
 
                                     HPX_UNUSED(zones_partition);
 
@@ -129,7 +130,8 @@ namespace lue {
                                     InputData const& input_partition_data,
                                     ZonesData const& zones_partition_data)
                                 {
-                                    AnnotateFunction annotation{"zonal_operation_partition"};
+                                    AnnotateFunction const annotation{
+                                        fmt::format("{}: partition", functor_name<Functor>)};
 
                                     HPX_UNUSED(input_partition);
                                     HPX_UNUSED(zones_partition);
@@ -188,10 +190,10 @@ namespace lue {
             @return     Partition with per zone the corresponding statistic
         */
         template<typename Policies, typename ZonesPartition, typename OutputPartition, typename Functor>
-        OutputPartition zonal_operation_partition2(
+        auto zonal_operation_partition2(
             Policies const& policies,
             ZonesPartition const& zones_partition,
-            AggregatorT<Functor> const& aggregator)
+            AggregatorT<Functor> const& aggregator) -> OutputPartition
         {
             using Offset = OffsetT<ZonesPartition>;
             using ZonesData = DataT<ZonesPartition>;
@@ -206,7 +208,8 @@ namespace lue {
                     [policies, zones_partition, aggregator](
                         Offset const& offset, ZonesData const& zones_partition_data)
                     {
-                        AnnotateFunction annotation{"zonal_operation_partition"};
+                        AnnotateFunction const annotation{
+                            fmt::format("{}: partition: reclass", functor_name<Functor>)};
 
                         HPX_UNUSED(zones_partition);
 
@@ -294,11 +297,7 @@ namespace lue {
                     [locality_id = localities[p], action, policies, functor](
                         hpx::shared_future<InputElement> const& input_element,
                         ZonesPartition const& zones_partition)
-                    {
-                        AnnotateFunction annotation{"zonal_operation"};
-
-                        return action(locality_id, policies, input_element.get(), zones_partition, functor);
-                    },
+                    { return action(locality_id, policies, input_element.get(), zones_partition, functor); },
 
                     input_element,
                     zones_partitions[p]);
@@ -315,7 +314,8 @@ namespace lue {
                 .then(
                     [](hpx::future<Aggregators>&& aggregators)
                     {
-                        AnnotateFunction annotation{"zonal_operation"};
+                        AnnotateFunction const annotation{
+                            fmt::format("{}: aggregate", functor_name<Functor>)};
 
                         Aggregator result{};
 
@@ -342,11 +342,7 @@ namespace lue {
 
                     [locality_id = localities[p], action, policies, zones_partition = zones_partitions[p]](
                         hpx::shared_future<Aggregator> const& aggregator)
-                    {
-                        AnnotateFunction annotation{"zonal_operation"};
-
-                        return action(locality_id, policies, zones_partition, aggregator.get());
-                    },
+                    { return action(locality_id, policies, zones_partition, aggregator.get()); },
 
                     aggregator);
             }
@@ -416,7 +412,8 @@ namespace lue {
                     [locality_id = localities[p], action, policies, functor](
                         InputPartition const& input_partition, ZonesPartition const& zones_partition)
                     {
-                        AnnotateFunction annotation{"zonal_operation"};
+                        AnnotateFunction const annotation{
+                            fmt::format("{}: partition: call statistic action", functor_name<Functor>)};
 
                         return action(locality_id, policies, input_partition, zones_partition, functor);
                     },
@@ -436,7 +433,8 @@ namespace lue {
                 .then(
                     [](hpx::future<Aggregators>&& aggregators)
                     {
-                        AnnotateFunction annotation{"zonal_operation"};
+                        AnnotateFunction const annotation{
+                            fmt::format("{}: array: merge", functor_name<Functor>)};
 
                         Aggregator result{};
 
@@ -464,7 +462,8 @@ namespace lue {
                     [locality_id = localities[p], action, policies, zones_partition = zones_partitions[p]](
                         hpx::shared_future<Aggregator> const& aggregator)
                     {
-                        AnnotateFunction annotation{"zonal_operation"};
+                        AnnotateFunction const annotation{
+                            fmt::format("{}: partition: call reclass action", functor_name<Functor>)};
 
                         return action(locality_id, policies, zones_partition, aggregator.get());
                     },
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_sum.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_sum.hpp
index 9562e4d80..b912a54af 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_sum.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/zonal_sum.hpp
@@ -18,6 +18,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<Zone>);
 
+                static constexpr char const* name{"zonal_sum"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/divide.hpp b/source/framework/algorithm/include/lue/framework/algorithm/divide.hpp
index 30cfd5734..1aba26699 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/divide.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/divide.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"divide"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/downstream_distance.hpp b/source/framework/algorithm/include/lue/framework/algorithm/downstream_distance.hpp
index db9e8231c..ad96ac29d 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/downstream_distance.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/downstream_distance.hpp
@@ -22,15 +22,16 @@ namespace lue {
 
 
     template<typename Policies, typename FlowDirectionElement, typename DistanceElement, Rank rank>
-    PartitionedArray<DistanceElement, rank> downstream_distance(
+    auto downstream_distance(
         Policies const& policies,
         PartitionedArray<FlowDirectionElement, rank> const& flow_direction,
-        DistanceElement const cell_size);
+        DistanceElement cell_size) -> PartitionedArray<DistanceElement, rank>;
 
 
     template<typename FlowDirectionElement, typename DistanceElement, Rank rank>
-    PartitionedArray<DistanceElement, rank> downstream_distance(
+    auto downstream_distance(
         PartitionedArray<FlowDirectionElement, rank> const& flow_direction, DistanceElement const cell_size)
+        -> PartitionedArray<DistanceElement, rank>
     {
         using Policies = policy::downstream_distance::DefaultPolicies<FlowDirectionElement, DistanceElement>;
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/equal_to.hpp b/source/framework/algorithm/include/lue/framework/algorithm/equal_to.hpp
index 2c92f4ba5..e53c7a89c 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/equal_to.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/equal_to.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"equal_to"};
+
                 using OutputElement = OutputElement_;
 
                 constexpr auto operator()(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/exp.hpp b/source/framework/algorithm/include/lue/framework/algorithm/exp.hpp
index 1eecb2ea9..4b483a2e7 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/exp.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/exp.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"exp"};
+
                 using OutputElement = InputElement;
 
                 auto operator()(InputElement const& input_element) const noexcept -> OutputElement
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/floor.hpp b/source/framework/algorithm/include/lue/framework/algorithm/floor.hpp
index eda35e766..ba6a2f4d4 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/floor.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/floor.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"floor"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/functor_traits.hpp b/source/framework/algorithm/include/lue/framework/algorithm/functor_traits.hpp
index e4fc45aa5..724ced295 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/functor_traits.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/functor_traits.hpp
@@ -17,11 +17,11 @@ namespace lue {
 
         public:
 
-            static constexpr bool const is_functor{false};
+            static constexpr char const* name{SomeType::name};
     };
 
 
     template<typename Functor>
-    constexpr bool is_functor_v = FunctorTraits<Functor>::is_functor;
+    constexpr char const* functor_name = FunctorTraits<Functor>::name;
 
 }  // namespace lue
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/greater_than.hpp b/source/framework/algorithm/include/lue/framework/algorithm/greater_than.hpp
index 79a352aa6..156fbb682 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/greater_than.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/greater_than.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"greater_than"};
+
                 using OutputElement = OutputElement_;
 
                 constexpr auto operator()(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/greater_than_equal_to.hpp b/source/framework/algorithm/include/lue/framework/algorithm/greater_than_equal_to.hpp
index d28c9ef45..3c1945ebc 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/greater_than_equal_to.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/greater_than_equal_to.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"greater_than_equal_to"};
+
                 using OutputElement = OutputElement_;
 
                 constexpr auto operator()(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/iterate_per_element.hpp b/source/framework/algorithm/include/lue/framework/algorithm/iterate_per_element.hpp
index 5b181ad33..fb87aa0e5 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/iterate_per_element.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/iterate_per_element.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_arithmetic_v<InputElement>);
 
+                static constexpr char const* name{"iterate_per_element"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/less_than.hpp b/source/framework/algorithm/include/lue/framework/algorithm/less_than.hpp
index 278c102b9..a2809c471 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/less_than.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/less_than.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"less_than"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/less_than_equal_to.hpp b/source/framework/algorithm/include/lue/framework/algorithm/less_than_equal_to.hpp
index deb28fbfe..6ae9f26f8 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/less_than_equal_to.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/less_than_equal_to.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"less_than_equal_to"};
+
                 using OutputElement = OutputElement_;
 
                 constexpr auto operator()(
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/log.hpp b/source/framework/algorithm/include/lue/framework/algorithm/log.hpp
index c46a89489..6945322ed 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/log.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/log.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"log"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/log10.hpp b/source/framework/algorithm/include/lue/framework/algorithm/log10.hpp
index 5d8f797a9..99010418b 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/log10.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/log10.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"log10"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/logical_and.hpp b/source/framework/algorithm/include/lue/framework/algorithm/logical_and.hpp
index 96e8046a2..dac5e2e78 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/logical_and.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/logical_and.hpp
@@ -14,6 +14,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<InputElement>);
 
+                static constexpr char const* name{"logical_and"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/logical_exclusive_or.hpp b/source/framework/algorithm/include/lue/framework/algorithm/logical_exclusive_or.hpp
index 053b7cef1..b375319b8 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/logical_exclusive_or.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/logical_exclusive_or.hpp
@@ -14,6 +14,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<InputElement>);
 
+                static constexpr char const* name{"logical_exclusive_or"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/logical_inclusive_or.hpp b/source/framework/algorithm/include/lue/framework/algorithm/logical_inclusive_or.hpp
index ddcf496bc..9d0a8c6e8 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/logical_inclusive_or.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/logical_inclusive_or.hpp
@@ -14,6 +14,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<InputElement>);
 
+                static constexpr char const* name{"logical_inclusive_or"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/logical_not.hpp b/source/framework/algorithm/include/lue/framework/algorithm/logical_not.hpp
index 010942090..abc330dff 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/logical_not.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/logical_not.hpp
@@ -14,6 +14,8 @@ namespace lue {
 
                 static_assert(std::is_integral_v<InputElement>);
 
+                static constexpr char const* name{"logical_not"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/modulus.hpp b/source/framework/algorithm/include/lue/framework/algorithm/modulus.hpp
index bfe281973..ab5e5dc3f 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/modulus.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/modulus.hpp
@@ -14,10 +14,13 @@ namespace lue {
                 static_assert(std::is_integral_v<InputElement>);
                 static_assert(std::is_same_v<InputElement, OutputElement_>);
 
+                static constexpr char const* name{"modulus"};
+
                 using OutputElement = OutputElement_;
 
                 constexpr auto operator()(
-                    InputElement const& input_element1, InputElement const& input_element2) const noexcept -> OutputElement
+                    InputElement const& input_element1,
+                    InputElement const& input_element2) const noexcept -> OutputElement
                 {
                     return input_element1 % input_element2;
                 }
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/multiply.hpp b/source/framework/algorithm/include/lue/framework/algorithm/multiply.hpp
index e65d11575..ff88d24c5 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/multiply.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/multiply.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"multiply"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/negate.hpp b/source/framework/algorithm/include/lue/framework/algorithm/negate.hpp
index a69ff7684..615ef11cb 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/negate.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/negate.hpp
@@ -14,6 +14,8 @@ namespace lue {
 
                 static_assert(std::is_signed_v<InputElement> || std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"negate"};
+
                 using OutputElement = InputElement;
 
                 auto operator()(InputElement const& input_element) const noexcept -> OutputElement
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/normal.hpp b/source/framework/algorithm/include/lue/framework/algorithm/normal.hpp
index ab6c1d40f..89dd56241 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/normal.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/normal.hpp
@@ -65,9 +65,8 @@ namespace lue {
                             // Standard mersenne_twister_engine seeded with the random_device
                             std::mt19937 random_number_engine(random_device());
 
-                            auto distribution = [mean, stddev]() {
-                                return std::normal_distribution<OutputElement>{mean, stddev};
-                            }();
+                            auto distribution = [mean, stddev]()
+                            { return std::normal_distribution<OutputElement>{mean, stddev}; }();
 
                             std::generate(
                                 output_partition_data.begin(),
@@ -266,16 +265,6 @@ namespace lue {
     };
 
 
-    template<typename Element, Rank rank>
-    class FunctorTraits<InstantiateNormal<Element, rank>>
-    {
-
-        public:
-
-            static constexpr bool const is_functor{true};
-    };
-
-
     template<typename Element, typename Policies, typename Shape>
     PartitionedArray<Element, rank<Shape>> normal(
         Policies const& policies,
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/not_equal_to.hpp b/source/framework/algorithm/include/lue/framework/algorithm/not_equal_to.hpp
index f6fcd1d68..f92ffc498 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/not_equal_to.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/not_equal_to.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"not_equal_to"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/pow.hpp b/source/framework/algorithm/include/lue/framework/algorithm/pow.hpp
index e1f819e72..61205ee56 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/pow.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/pow.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"pow"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/round.hpp b/source/framework/algorithm/include/lue/framework/algorithm/round.hpp
index 26d0add0b..4471fb855 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/round.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/round.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"round"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/sin.hpp b/source/framework/algorithm/include/lue/framework/algorithm/sin.hpp
index 9d5b12bf5..5080ddaba 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/sin.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/sin.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"sin"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/sqrt.hpp b/source/framework/algorithm/include/lue/framework/algorithm/sqrt.hpp
index 9e04bf8c8..10e151914 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/sqrt.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/sqrt.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"sqrt"};
+
                 using OutputElement = InputElement;
 
                 auto operator()(InputElement const& input_element) const noexcept -> OutputElement
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/subtract.hpp b/source/framework/algorithm/include/lue/framework/algorithm/subtract.hpp
index 5a966ab1f..8233a401b 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/subtract.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/subtract.hpp
@@ -13,6 +13,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"subtract"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/tan.hpp b/source/framework/algorithm/include/lue/framework/algorithm/tan.hpp
index f74f744ba..b1e5c5ccb 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/tan.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/tan.hpp
@@ -15,6 +15,8 @@ namespace lue {
 
                 static_assert(std::is_floating_point_v<InputElement>);
 
+                static constexpr char const* name{"tan"};
+
                 using OutputElement = InputElement;
 
 
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/uniform.hpp b/source/framework/algorithm/include/lue/framework/algorithm/uniform.hpp
index dbc5bcac7..6b7c3a021 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/uniform.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/uniform.hpp
@@ -314,16 +314,6 @@ namespace lue {
     };
 
 
-    template<typename Element, Rank rank>
-    class FunctorTraits<InstantiateUniform<Element, rank>>
-    {
-
-        public:
-
-            static constexpr bool const is_functor{true};
-    };
-
-
     template<typename Policies, typename Shape>
     requires(!std::is_same_v<policy::OutputElementT<Policies, 0>, std::uint8_t> && !std::is_same_v<policy::OutputElementT<Policies, 0>, std::int8_t>) auto uniform(
         Policies const& policies,
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/valid.hpp b/source/framework/algorithm/include/lue/framework/algorithm/valid.hpp
index 62e537cfd..40ba2395e 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/valid.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/valid.hpp
@@ -12,6 +12,8 @@ namespace lue {
 
             public:
 
+                static constexpr char const* name{"valid"};
+
                 using OutputElement = OutputElement_;
 
 
diff --git a/source/framework/core/include/lue/framework/core/annotate.hpp b/source/framework/core/include/lue/framework/core/annotate.hpp
index e3b889f37..ee523b892 100644
--- a/source/framework/core/include/lue/framework/core/annotate.hpp
+++ b/source/framework/core/include/lue/framework/core/annotate.hpp
@@ -10,7 +10,7 @@ namespace lue {
 #ifdef HPX_HAVE_APEX
     using AnnotateFunction = hpx::scoped_annotation;
 #else
-#define AnnotateFunction [[maybe_unused]] static char const*
+#define AnnotateFunction [[maybe_unused]] static auto
 #endif
 
 }  // namespace lue
diff --git a/source/framework/io/src/raster.cpp b/source/framework/io/src/raster.cpp
index 9bc2154a8..9ffc6f9ee 100644
--- a/source/framework/io/src/raster.cpp
+++ b/source/framework/io/src/raster.cpp
@@ -394,15 +394,6 @@ namespace lue {
     }  // Anonymous namespace
 
 
-    template<typename Element>
-    class FunctorTraits<ReadPartitionsPerLocality<Element>>
-    {
-        public:
-
-            static constexpr bool is_functor{true};
-    };
-
-
     template<typename Element>
     auto read(std::string const& name, Shape<Count, 2> const& partition_shape) -> PartitionedArray<Element, 2>
     {
@@ -478,8 +469,9 @@ namespace lue {
 
     template<typename Element>
     auto write(
-        PartitionedArray<Element, 2> const& array, std::string const& name, std::string const& clone_name)
-        -> hpx::future<void>
+        PartitionedArray<Element, 2> const& array,
+        std::string const& name,
+        std::string const& clone_name) -> hpx::future<void>
     {
         if constexpr (!gdal::supports<Element>())
         {
diff --git a/source/framework/python/CMakeLists.txt b/source/framework/python/CMakeLists.txt
index 985ff238c..57e997d82 100644
--- a/source/framework/python/CMakeLists.txt
+++ b/source/framework/python/CMakeLists.txt
@@ -4,7 +4,7 @@ add_test_conditionally(test)
 
 
 # framework/*
-foreach(module_name __init__.py support.py)
+foreach(module_name __init__.py pcraster_provider.py support.py)
     file(GENERATE
         OUTPUT ${CMAKE_BINARY_DIR}/${LUE_BUILD_PYTHON_PACKAGE_DIR}/$<$<BOOL:${LUE_GENERATOR_IS_MULTI_CONFIG}>:$<CONFIG>>/framework/${module_name}
         INPUT framework/${module_name}
diff --git a/source/framework/python/framework/__init__.py b/source/framework/python/framework/__init__.py
index 4640003f9..70393ee93 100644
--- a/source/framework/python/framework/__init__.py
+++ b/source/framework/python/framework/__init__.py
@@ -1,10 +1,12 @@
 from ..lue_py import framework as lfr
 from ..lue_py.framework import *
 from .support import (
+    DefaultProgressor,
     DynamicModel,
     DynamicModelRunner,
     MonteCarloModel,
     MonteCarloModelRunner,
+    SilentProgressor,
     StaticModel,
     StaticModelRunner,
 )
diff --git a/source/framework/python/framework/pcraster_provider.py b/source/framework/python/framework/pcraster_provider.py
new file mode 100644
index 000000000..c423251bf
--- /dev/null
+++ b/source/framework/python/framework/pcraster_provider.py
@@ -0,0 +1,48 @@
+import os
+
+
+def pcraster_by_lue():
+    import lue.framework as lfr
+    import lue.pcraster as pcr
+    import lue.pcraster.framework as pcrfw
+
+    pcr.runtime_scope = lfr.runtime_scope
+
+    return pcr, pcrfw
+
+
+def pcraster_by_pcraster():
+    import pcraster as pcr
+    import pcraster.framework as pcrfw
+
+    def runtime_scope(user_main):
+        def decorated_function(*args, **kwargs):
+            user_main(*args, **kwargs)
+
+        return decorated_function
+
+    pcr.runtime_scope = runtime_scope
+
+    return pcr, pcrfw
+
+
+def pcraster_by_provider():
+    pcraster_by_provider_name = {
+        "lue": pcraster_by_lue,
+        "pcraster": pcraster_by_pcraster,
+    }
+
+    provider_name = os.getenv("LUE_PCRASTER_PROVIDER_NAME", default="unknown")
+
+    if provider_name not in pcraster_by_provider_name:
+        raise RuntimeError(
+            f'Select a PCRaster provider by setting the LUE_PCRASTER_PROVIDER_NAME environment variable to one of {", ".join(pcraster_by_provider_name)}'
+        )
+
+    pcr, pcrfw = pcraster_by_provider_name[provider_name]()
+    pcr.provider_name = provider_name
+
+    return pcr, pcrfw
+
+
+pcr, pcrfw = pcraster_by_provider()
diff --git a/source/framework/python/framework/support.py b/source/framework/python/framework/support.py
index d9d1079c3..b318ef4cf 100644
--- a/source/framework/python/framework/support.py
+++ b/source/framework/python/framework/support.py
@@ -5,7 +5,7 @@
 import lue.framework as lfr
 
 
-class Progressor(lfr.Progressor):
+class DefaultProgressor(lfr.Progressor):
     def __init__(self):
         lfr.Progressor.__init__(self)
 
@@ -30,6 +30,11 @@ def postprocess(self):
         sys.stdout.flush()
 
 
+class SilentProgressor(lfr.Progressor):
+    def __init__(self):
+        lfr.Progressor.__init__(self)
+
+
 class StaticModel(lfr.Model):
     def __init__(self):
         lfr.Model.__init__(self)
@@ -45,9 +50,9 @@ class StaticModelRunner(object):
     def __init__(self, model):
         self.model = model
 
-    def run(self):
-        # self.model.run_initial()
-        lfr.run_deterministic(self.model, Progressor(), 0)
+    def run(self, *, progressor=DefaultProgressor(), rate_limit=0):
+        assert rate_limit >= 0, rate_limit
+        lfr.run_deterministic(self.model, progressor, 0, rate_limit)
 
 
 class DynamicModel(lfr.Model):
@@ -97,15 +102,11 @@ def __init__(self, model, last_time_step, first_time_step=1):
             1 <= first_time_step <= last_time_step
         ), f"{first_time_step}, {last_time_step}"
 
-    def run(self):
+    def run(self, *, progressor=DefaultProgressor(), rate_limit=0):
+        assert rate_limit >= 0, rate_limit
         assert self.first_time_step == 1, self.first_time_step
-        lfr.run_deterministic(self.model, Progressor(), self.last_time_step)
-
-        # self.model.run_initial()
 
-        # for time_step in range(self.first_time_step, self.last_time_step + 1):
-        #     self.model.current_time_step = time_step
-        #     self.model.run_dynamic()
+        lfr.run_deterministic(self.model, progressor, self.last_time_step, rate_limit)
 
 
 class MonteCarloModel(object):  # (lfr.Model):
@@ -237,7 +238,7 @@ def __init__(self, framework_model, nr_samples, remove_existing_directories):
             self.framework_model.model, remove_existing_directories
         )
 
-    def run(self, rate_limit=0):
+    def run(self, *, progressor=DefaultProgressor(), rate_limit=0):
         assert rate_limit >= 0, rate_limit
         assert (
             self.framework_model.first_time_step == 1
@@ -245,7 +246,7 @@ def run(self, rate_limit=0):
 
         lfr.run_stochastic(
             self.framework_model.model,
-            Progressor(),
+            progressor,
             self.nr_samples,
             self.framework_model.last_time_step,
             rate_limit,
diff --git a/source/framework/python/include/lue/py/framework/core/partitioned_array.hpp b/source/framework/python/include/lue/py/framework/core/partitioned_array.hpp
index caa9e931c..c372d8834 100644
--- a/source/framework/python/include/lue/py/framework/core/partitioned_array.hpp
+++ b/source/framework/python/include/lue/py/framework/core/partitioned_array.hpp
@@ -21,24 +21,27 @@
 
 
 namespace lue::framework {
+    namespace detail {
 
-    template<typename Element, Rank rank>
-    static auto formal_string_representation(PartitionedArray<Element, rank> const& array) -> std::string
-    {
-        return fmt::format(
-            "PartitionedArray<{}, {}>{{shape={}, shape_in_partitions={}}}",
-            as_string<Element>(),
-            rank,
-            as_string(array.shape()),
-            as_string(array.partitions().shape()));
-    }
+        template<typename Element, Rank rank>
+        static auto formal_string_representation(PartitionedArray<Element, rank> const& array) -> std::string
+        {
+            return fmt::format(
+                "PartitionedArray<{}, {}>{{shape={}, shape_in_partitions={}}}",
+                as_string<Element>(),
+                rank,
+                as_string(array.shape()),
+                as_string(array.partitions().shape()));
+        }
 
 
-    template<typename Element, Rank rank>
-    static auto informal_string_representation(PartitionedArray<Element, rank> const& array) -> std::string
-    {
-        return formal_string_representation(array);
-    }
+        template<typename Element, Rank rank>
+        static auto informal_string_representation(PartitionedArray<Element, rank> const& array) -> std::string
+        {
+            return formal_string_representation(array);
+        }
+
+    }  // namespace detail
 
 
     template<typename Element, Rank rank>
@@ -101,9 +104,9 @@ namespace lue::framework {
                             .then([](auto&&) { return hpx::make_ready_future<void>(); });
                     })
 
-                .def("__repr__", [](Array const& array) { return formal_string_representation(array); })
+                .def("__repr__", [](Array const& array) { return detail::formal_string_representation(array); })
 
-                .def("__str__", [](Array const& array) { return informal_string_representation(array); })
+                .def("__str__", [](Array const& array) { return detail::informal_string_representation(array); })
 
                 // bool(a), not a, if a, while a, ...
                 .def(
@@ -113,151 +116,150 @@ namespace lue::framework {
                         // ValueError
                         throw std::invalid_argument("The truth value of an array is ambiguous");
                     })
+                ;
 
-                // a < b
-                .def(
-                    "__lt__",
-                    [](Array const& argument1, Array const& argument2)
-                    { return less_than<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__lt__",
-                    [](Array const& argument1, Element const argument2)
-                    { return less_than<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__lt__",
-                    [](Array const& argument1, Scalar const& argument2)
-                    { return less_than<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-
-                // a <= b
-                .def(
-                    "__le__",
-                    [](Array const& argument1, Array const& argument2)
-                    { return less_than_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__le__",
-                    [](Array const& argument1, Element const argument2)
-                    { return less_than_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__le__",
-                    [](Array const& argument1, Scalar const& argument2)
-                    { return less_than_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-
-                // a == b
-                .def(
-                    "__eq__",
-                    [](Array const& argument1, Array const& argument2)
-                    { return equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__eq__",
-                    [](Array const& argument1, Element const argument2)
-                    { return equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__eq__",
-                    [](Array const& argument1, Scalar const& argument2)
-                    { return equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
+        class_
+            // a < b
+            .def(
+                "__lt__",
+                [](Array const& argument1, Array const& argument2)
+                { return less_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__lt__",
+                [](Array const& argument1, Scalar const& argument2)
+                { return less_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__lt__",
+                [](Array const& argument1, Element const argument2)
+                { return less_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
 
-                // a != b
-                .def(
-                    "__ne__",
-                    [](Array const& argument1, Array const& argument2)
-                    { return not_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__ne__",
-                    [](Array const& argument1, Element const argument2)
-                    { return not_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__ne__",
-                    [](Array const& argument1, Scalar const& argument2)
-                    { return not_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
+            // a <= b
+            .def(
+                "__le__",
+                [](Array const& argument1, Array const& argument2)
+                { return less_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__le__",
+                [](Array const& argument1, Scalar const& argument2)
+                { return less_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__le__",
+                [](Array const& argument1, Element const argument2)
+                { return less_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
 
-                // a > b
-                .def(
-                    "__gt__",
-                    [](Array const& argument1, Array const& argument2)
-                    { return greater_than<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__gt__",
-                    [](Array const& argument1, Element const argument2)
-                    { return greater_than<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__gt__",
-                    [](Array const& argument1, Scalar const& argument2)
-                    { return greater_than<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
+            // a == b
+            .def(
+                "__eq__",
+                [](Array const& argument1, Array const& argument2)
+                { return equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__eq__",
+                [](Array const& argument1, Scalar const& argument2)
+                { return equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__eq__",
+                [](Array const& argument1, Element const argument2)
+                { return equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
 
-                // a >= b
-                .def(
-                    "__ge__",
-                    [](Array const& argument1, Array const& argument2)
-                    { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__ge__",
-                    [](Array const& argument1, Element const argument2)
-                    { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__ge__",
-                    [](Array const& argument1, Scalar const& argument2)
-                    { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
-                    pybind11::is_operator())
+            // a != b
+            .def(
+                "__ne__",
+                [](Array const& argument1, Array const& argument2)
+                { return not_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__ne__",
+                [](Array const& argument1, Scalar const& argument2)
+                { return not_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__ne__",
+                [](Array const& argument1, Element const argument2)
+                { return not_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
 
-            ;
+            // a > b
+            .def(
+                "__gt__",
+                [](Array const& argument1, Array const& argument2)
+                { return greater_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__gt__",
+                [](Array const& argument1, Scalar const& argument2)
+                { return greater_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__gt__",
+                [](Array const& argument1, Element const argument2)
+                { return greater_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
 
+            // a >= b
+            .def(
+                "__ge__",
+                [](Array const& argument1, Array const& argument2)
+                { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__ge__",
+                [](Array const& argument1, Scalar const& argument2)
+                { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__ge__",
+                [](Array const& argument1, Element const argument2)
+                { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
 
-        // a + b, b + a, a += b
-        class_
+            // a + b, b + a, a += b
             .def(
                 "__add__",
                 [](Array const& argument1, Array const& argument2) { return add(argument1, argument2); },
                 pybind11::is_operator())
             .def(
                 "__add__",
-                [](Array const& argument1, Element const argument2) { return add(argument1, argument2); },
+                [](Array const& argument1, Scalar const& argument2) { return add(argument1, argument2); },
                 pybind11::is_operator())
             .def(
                 "__add__",
-                [](Array const& argument1, Scalar const& argument2) { return add(argument1, argument2); },
+                [](Array const& argument1, Element const argument2) { return add(argument1, argument2); },
                 pybind11::is_operator())
+
             .def(
                 "__radd__",
-                [](Array const& argument2, Element const argument1) { return add(argument1, argument2); },
+                [](Array const& argument2, Scalar const& argument1) { return add(argument1, argument2); },
                 pybind11::is_operator())
             .def(
                 "__radd__",
-                [](Array const& argument2, Scalar const& argument1) { return add(argument1, argument2); },
-                pybind11::is_operator());
+                [](Array const& argument2, Element const argument1) { return add(argument1, argument2); },
+                pybind11::is_operator())
 
-        // a - b, b - a, a -= b
-        class_
+            // a - b, b - a, a -= b
             .def(
                 "__sub__",
                 [](Array const& argument1, Array const& argument2) { return subtract(argument1, argument2); },
                 pybind11::is_operator())
             .def(
                 "__sub__",
-                [](Array const& argument1, Element const argument2)
+                [](Array const& argument1, Scalar const& argument2)
                 { return subtract(argument1, argument2); },
                 pybind11::is_operator())
             .def(
                 "__sub__",
-                [](Array const& argument1, Scalar const& argument2)
+                [](Array const& argument1, Element const argument2)
                 { return subtract(argument1, argument2); },
                 pybind11::is_operator())
+
             .def(
                 "__rsub__",
                 [](Array const& argument2, Element const argument1)
@@ -267,50 +269,118 @@ namespace lue::framework {
                 "__rsub__",
                 [](Array const& argument2, Scalar const& argument1)
                 { return subtract(argument1, argument2); },
-                pybind11::is_operator());
+                pybind11::is_operator())
+            ;
 
         if constexpr (std::is_signed_v<Element> || std::is_floating_point_v<Element>)
         {
-            // abs(a)
-            class_.def(
-                "__abs__", [](Array const& self) { return abs(self); }, pybind11::is_operator());
+            class_
+
+                // abs(a)
+                .def(
+                    "__abs__", [](Array const& self) { return abs(self); }, pybind11::is_operator())
 
-            // -a
-            class_.def(
-                "__neg__", [](Array const& self) { return negate(self); }, pybind11::is_operator());
+                // -a
+                .def(
+                    "__neg__", [](Array const& self) { return negate(self); }, pybind11::is_operator())
+                ;
         }
 
         if constexpr (std::is_integral_v<Element>)
         {
-            // a | b
-            class_.def(
-                "__or__",
-                [](Array const& argument1, Array const& argument2)
-                { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
-                pybind11::is_operator());
+            class_
 
-            // a ^ b
-            class_.def(
-                "__xor__",
-                [](Array const& argument1, Array const& argument2)
-                { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
-                pybind11::is_operator());
+                // a | b
+                .def(
+                    "__or__",
+                    [](Array const& argument1, Array const& argument2)
+                    { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__or__",
+                    [](Array const& argument1, Scalar const& argument2)
+                    { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__or__",
+                    [](Array const& argument1, Element const& argument2)
+                    { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
 
-            // a & b
-            class_.def(
-                "__and__",
-                [](Array const& argument1, Array const& argument2)
-                { return logical_and<BooleanElement>(argument1, argument2); },
-                pybind11::is_operator());
+                .def(
+                    "__ror__",
+                    [](Array const& argument1, Scalar const& argument2)
+                    { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__ror__",
+                    [](Array const& argument1, Element const& argument2)
+                    { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+
+                // a ^ b
+                .def(
+                    "__xor__",
+                    [](Array const& argument1, Array const& argument2)
+                    { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__xor__",
+                    [](Array const& argument1, Scalar const& argument2)
+                    { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__xor__",
+                    [](Array const& argument1, Element const& argument2)
+                    { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
 
-            // ~a
-            class_.def(
-                "__invert__",
-                [](Array const& self) { return logical_not<BooleanElement>(self); },
-                pybind11::is_operator());
+                .def(
+                    "__rxor__",
+                    [](Array const& argument1, Scalar const& argument2)
+                    { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__rxor__",
+                    [](Array const& argument1, Element const& argument2)
+                    { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
 
-            // a % b
-            class_
+                // a & b
+                .def(
+                    "__and__",
+                    [](Array const& argument1, Array const& argument2)
+                    { return logical_and<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__and__",
+                    [](Array const& argument1, Scalar const& argument2)
+                    { return logical_and<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__and__",
+                    [](Array const& argument1, Element const& argument2)
+                    { return logical_and<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+
+                .def(
+                    "__rand__",
+                    [](Array const& argument1, Scalar const& argument2)
+                    { return logical_and<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__rand__",
+                    [](Array const& argument1, Element const& argument2)
+                    { return logical_and<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+
+                // ~a
+                .def(
+                    "__invert__",
+                    [](Array const& self) { return logical_not<BooleanElement>(self); },
+                    pybind11::is_operator())
+
+                // a % b
                 .def(
                     "__mod__",
                     [](Array const& argument1, Array const& argument2)
@@ -318,30 +388,32 @@ namespace lue::framework {
                     pybind11::is_operator())
                 .def(
                     "__mod__",
-                    [](Array const& argument1, Element const argument2)
+                    [](Array const& argument1, Scalar const& argument2)
                     { return modulus(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__mod__",
-                    [](Array const& argument1, Scalar const& argument2)
+                    [](Array const& argument1, Element const argument2)
                     { return modulus(argument1, argument2); },
                     pybind11::is_operator())
+
                 .def(
                     "__rmod__",
-                    [](Array const& argument2, Element const argument1)
+                    [](Array const& argument2, Scalar const& argument1)
                     { return modulus(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__rmod__",
-                    [](Array const& argument2, Scalar const& argument1)
+                    [](Array const& argument2, Element const argument1)
                     { return modulus(argument1, argument2); },
-                    pybind11::is_operator());
+                    pybind11::is_operator())
+                ;
         }
 
         if constexpr (std::is_floating_point_v<Element>)
         {
-            // a * b, b * a, a *= b
             class_
+                // a * b, b * a, a *= b
                 .def(
                     "__mul__",
                     [](Array const& argument1, Array const& argument2)
@@ -349,27 +421,27 @@ namespace lue::framework {
                     pybind11::is_operator())
                 .def(
                     "__mul__",
-                    [](Array const& argument1, Element const argument2)
+                    [](Array const& argument1, Scalar const& argument2)
                     { return multiply(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__mul__",
-                    [](Array const& argument1, Scalar const& argument2)
+                    [](Array const& argument1, Element const argument2)
                     { return multiply(argument1, argument2); },
                     pybind11::is_operator())
+
                 .def(
                     "__rmul__",
-                    [](Array const& argument2, Element const argument1)
+                    [](Array const& argument2, Scalar const& argument1)
                     { return multiply(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__rmul__",
-                    [](Array const& argument2, Scalar const& argument1)
+                    [](Array const& argument2, Element const argument1)
                     { return multiply(argument1, argument2); },
-                    pybind11::is_operator());
+                    pybind11::is_operator())
 
-            // a / b, b / a, a /= b
-            class_
+                // a / b, b / a, a /= b
                 .def(
                     "__truediv__",
                     [](Array const& argument1, Array const& argument2)
@@ -377,47 +449,49 @@ namespace lue::framework {
                     pybind11::is_operator())
                 .def(
                     "__truediv__",
-                    [](Array const& argument1, Element const argument2)
+                    [](Array const& argument1, Scalar const& argument2)
                     { return divide(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__truediv__",
-                    [](Array const& argument1, Scalar const& argument2)
+                    [](Array const& argument1, Element const argument2)
                     { return divide(argument1, argument2); },
                     pybind11::is_operator())
+
                 .def(
                     "__rtruediv__",
-                    [](Array const& argument2, Element const argument1)
+                    [](Array const& argument2, Scalar const& argument1)
                     { return divide(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__rtruediv__",
-                    [](Array const& argument2, Scalar const& argument1)
+                    [](Array const& argument2, Element const argument1)
                     { return divide(argument1, argument2); },
-                    pybind11::is_operator());
+                    pybind11::is_operator())
 
-            // a ** b, b ** a, a **= b
-            class_
+                // a ** b, b ** a, a **= b
                 .def(
                     "__pow__",
                     [](Array const& argument1, Array const& argument2) { return pow(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__pow__",
-                    [](Array const& argument1, Element const argument2) { return pow(argument1, argument2); },
+                    [](Array const& argument1, Scalar const& argument2) { return pow(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__pow__",
-                    [](Array const& argument1, Scalar const& argument2) { return pow(argument1, argument2); },
+                    [](Array const& argument1, Element const argument2) { return pow(argument1, argument2); },
                     pybind11::is_operator())
+
                 .def(
                     "__rpow__",
-                    [](Array const& argument2, Element const argument1) { return pow(argument1, argument2); },
+                    [](Array const& argument2, Scalar const& argument1) { return pow(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
                     "__rpow__",
-                    [](Array const& argument2, Scalar const& argument1) { return pow(argument1, argument2); },
-                    pybind11::is_operator());
+                    [](Array const& argument2, Element const argument1) { return pow(argument1, argument2); },
+                    pybind11::is_operator())
+                ;
         }
     }
 
diff --git a/source/framework/python/include/lue/py/framework/core/scalar.hpp b/source/framework/python/include/lue/py/framework/core/scalar.hpp
index e151beb18..0990cf18c 100644
--- a/source/framework/python/include/lue/py/framework/core/scalar.hpp
+++ b/source/framework/python/include/lue/py/framework/core/scalar.hpp
@@ -1,146 +1,53 @@
-// #include "lue/py/framework/algorithm/abs.hpp"
+#include "lue/framework/algorithm/value_policies/abs.hpp"
 #include "lue/framework/algorithm/value_policies/add.hpp"
 #include "lue/framework/algorithm/value_policies/divide.hpp"
-// #include "lue/py/framework/algorithm/equal_to.hpp"
-// #include "lue/py/framework/algorithm/greater_than.hpp"
-// #include "lue/py/framework/algorithm/greater_than_equal_to.hpp"
-// #include "lue/py/framework/algorithm/less_than.hpp"
-// #include "lue/py/framework/algorithm/less_than_equal_to.hpp"
-// #include "lue/py/framework/algorithm/logical_and.hpp"
-// #include "lue/py/framework/algorithm/logical_exclusive_or.hpp"
-// #include "lue/py/framework/algorithm/logical_inclusive_or.hpp"
-// #include "lue/py/framework/algorithm/logical_not.hpp"
+#include "lue/framework/algorithm/value_policies/equal_to.hpp"
+#include "lue/framework/algorithm/value_policies/greater_than.hpp"
+#include "lue/framework/algorithm/value_policies/greater_than_equal_to.hpp"
+#include "lue/framework/algorithm/value_policies/less_than.hpp"
+#include "lue/framework/algorithm/value_policies/less_than_equal_to.hpp"
+#include "lue/framework/algorithm/value_policies/logical_and.hpp"
+#include "lue/framework/algorithm/value_policies/logical_exclusive_or.hpp"
+#include "lue/framework/algorithm/value_policies/logical_inclusive_or.hpp"
+#include "lue/framework/algorithm/value_policies/logical_not.hpp"
+#include "lue/framework/algorithm/value_policies/modulus.hpp"
 #include "lue/framework/algorithm/value_policies/multiply.hpp"
-// #include "lue/py/framework/algorithm/negate.hpp"
-// #include "lue/py/framework/algorithm/not_equal_to.hpp"
+#include "lue/framework/algorithm/value_policies/negate.hpp"
+#include "lue/framework/algorithm/value_policies/not_equal_to.hpp"
 #include "lue/framework/algorithm/value_policies/pow.hpp"
 #include "lue/framework/algorithm/value_policies/subtract.hpp"
 #include "lue/py/framework/stream.hpp"
 #include <pybind11/pybind11.h>
 
 
-using namespace pybind11::literals;
-
-
 namespace lue::framework {
-
-    namespace {
+    namespace detail {
 
         template<typename Element>
-        static std::string formal_string_representation([[maybe_unused]] Scalar<Element> const& scalar)
+        auto formal_string_representation([[maybe_unused]] Scalar<Element> const& scalar) -> std::string
         {
             return fmt::format("Scalar<{}>{{}}", as_string<Element>());
         }
 
 
         template<typename Element>
-        static std::string informal_string_representation(Scalar<Element> const& scalar)
+        auto informal_string_representation(Scalar<Element> const& scalar) -> std::string
         {
             return formal_string_representation(scalar);
         }
 
-
-        template<typename Element>
-        auto create_scalar(pybind11::object const& value) -> Scalar<Element>
-        {
-            return pybind11::cast<Element>(value);
-        }
-
-
-        auto create_scalar_py(pybind11::object const& dtype_args, pybind11::object const& value)
-            -> pybind11::object
-        {
-            pybind11::dtype const dtype{pybind11::dtype::from_args(dtype_args)};
-
-            // Switch on dtype and return a scalar of the right type
-
-            auto const kind = dtype.kind();
-            auto const size = dtype.itemsize();  // bytes
-            pybind11::object result;
-
-            switch (kind)
-            {
-                case 'i':
-                {
-                    // Signed integer
-                    switch (size)
-                    {
-                        case 4:
-                        {
-                            result = pybind11::cast(create_scalar<std::int32_t>(value));
-                            break;
-                        }
-                        case 8:
-                        {
-                            result = pybind11::cast(create_scalar<std::int64_t>(value));
-                            break;
-                        }
-                    }
-
-                    break;
-                }
-                case 'u':
-                {
-                    // Unsigned integer
-                    switch (size)
-                    {
-                        case 1:
-                        {
-                            result = pybind11::cast(create_scalar<std::uint8_t>(value));
-                            break;
-                        }
-                        case 4:
-                        {
-                            result = pybind11::cast(create_scalar<std::uint32_t>(value));
-                            break;
-                        }
-                        case 8:
-                        {
-                            result = pybind11::cast(create_scalar<std::uint64_t>(value));
-                            break;
-                        }
-                    }
-
-                    break;
-                }
-                case 'f':
-                {
-                    // Floating-point
-                    switch (size)
-                    {
-                        case 4:
-                        {
-                            result = pybind11::cast(create_scalar<float>(value));
-                            break;
-                        }
-                        case 8:
-                        {
-                            result = pybind11::cast(create_scalar<double>(value));
-                            break;
-                        }
-                    }
-
-                    break;
-                }
-            }
-
-            if (!result)
-            {
-                throw std::runtime_error(fmt::format("Unsupported dtype (kind={}, itemsize={})", kind, size));
-            }
-
-            return result;
-        }
-
-    }  // Anonymous namespace
+    }  // namespace detail
 
 
     template<typename Element>
     void bind_scalar(pybind11::module& module)
     {
+        // https://docs.python.org/3/reference/datamodel.html
         using namespace lue::value_policies;
+
+        Rank const rank{2};
         using Scalar = lue::Scalar<Element>;
-        using Value = Element;
+        using BooleanElement = std::uint8_t;
 
         auto class_ =
             pybind11::class_<Scalar>(
@@ -165,9 +72,9 @@ namespace lue::framework {
                     [](Scalar const& self) -> hpx::shared_future<Element> { return self.future(); },
                     "future docstring...")
 
-                .def("__repr__", [](Scalar const& scalar) { return formal_string_representation(scalar); })
+                .def("__repr__", [](Scalar const& scalar) { return detail::formal_string_representation(scalar); })
 
-                .def("__str__", [](Scalar const& scalar) { return informal_string_representation(scalar); })
+                .def("__str__", [](Scalar const& scalar) { return detail::informal_string_representation(scalar); })
 
                 // bool(a), not a, if a, while a, ...
                 .def(
@@ -177,141 +84,97 @@ namespace lue::framework {
                         // ValueError
                         throw std::invalid_argument("The truth value of a scalar is ambiguous");
                     })
+                ;
 
-            //                 // a < b
-            //                 .def(
-            //                     "__lt__",
-            //                     [](Scalar const& argument1, Scalar const& argument2)
-            //                     { return lfr::less_than<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__lt__",
-            //                     [](Scalar const& argument1, Element const argument2)
-            //                     { return lfr::less_than<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__lt__",
-            //                     [](Scalar const& argument1, Value const& argument2)
-            //                     { return lfr::less_than<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //
-            //                 // a <= b
-            //                 .def(
-            //                     "__le__",
-            //                     [](Scalar const& argument1, Scalar const& argument2) {
-            //                         return lfr::less_than_equal_to<ElementT<Scalar>, std::uint8_t>(
-            //                             argument1, argument2);
-            //                     },
-            //                     pybind11::is_operator())
-            //                 .def(
-            //                     "__le__",
-            //                     [](Scalar const& argument1, Element const argument2) {
-            //                         return lfr::less_than_equal_to<ElementT<Scalar>, std::uint8_t>(
-            //                             argument1, argument2);
-            //                     },
-            //                     pybind11::is_operator())
-            //                 .def(
-            //                     "__le__",
-            //                     [](Scalar const& argument1, Value const& argument2) {
-            //                         return lfr::less_than_equal_to<ElementT<Scalar>, std::uint8_t>(
-            //                             argument1, argument2);
-            //                     },
-            //                     pybind11::is_operator())
-            //
-            //                 // a == b
-            //                 .def(
-            //                     "__eq__",
-            //                     [](Scalar const& argument1, Scalar const& argument2)
-            //                     { return lfr::equal_to<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__eq__",
-            //                     [](Scalar const& argument1, Element const argument2)
-            //                     { return lfr::equal_to<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__eq__",
-            //                     [](Scalar const& argument1, Value const& argument2)
-            //                     { return lfr::equal_to<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //
-            //                 // a != b
-            //                 .def(
-            //                     "__ne__",
-            //                     [](Scalar const& argument1, Scalar const& argument2)
-            //                     { return lfr::not_equal_to<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__ne__",
-            //                     [](Scalar const& argument1, Element const argument2)
-            //                     { return lfr::not_equal_to<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__ne__",
-            //                     [](Scalar const& argument1, Value const& argument2)
-            //                     { return lfr::not_equal_to<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //
-            //                 // a > b
-            //                 .def(
-            //                     "__gt__",
-            //                     [](Scalar const& argument1, Scalar const& argument2)
-            //                     { return lfr::greater_than<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__gt__",
-            //                     [](Scalar const& argument1, Element const argument2)
-            //                     { return lfr::greater_than<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //                 .def(
-            //                     "__gt__",
-            //                     [](Scalar const& argument1, Value const& argument2)
-            //                     { return lfr::greater_than<ElementT<Scalar>, std::uint8_t>(argument1,
-            //                     argument2); }, pybind11::is_operator())
-            //
-            //                 // a >= b
-            //                 .def(
-            //                     "__ge__",
-            //                     [](Scalar const& argument1, Scalar const& argument2) {
-            //                         return lfr::greater_than_equal_to<ElementT<Scalar>, std::uint8_t>(
-            //                             argument1, argument2);
-            //                     },
-            //                     pybind11::is_operator())
-            //                 .def(
-            //                     "__ge__",
-            //                     [](Scalar const& argument1, Element const argument2) {
-            //                         return lfr::greater_than_equal_to<ElementT<Scalar>, std::uint8_t>(
-            //                             argument1, argument2);
-            //                     },
-            //                     pybind11::is_operator())
-            //                 .def(
-            //                     "__ge__",
-            //                     [](Scalar const& argument1, Value const& argument2) {
-            //                         return lfr::greater_than_equal_to<ElementT<Scalar>, std::uint8_t>(
-            //                             argument1, argument2);
-            //                     },
-            //                     pybind11::is_operator())
-
-            ;
-
-
-        // a + b, b + a, a += b
         class_
+            // a < b
+            .def(
+                "__lt__",
+                [](Scalar const& argument1, Scalar const& argument2)
+                { return less_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__lt__",
+                [](Scalar const& argument1, Element const& argument2)
+                { return less_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+
+            // a <= b
+            .def(
+                "__le__",
+                [](Scalar const& argument1, Scalar const& argument2)
+                { return less_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__le__",
+                [](Scalar const& argument1, Element const& argument2)
+                { return less_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+
+            // a == b
+            .def(
+                "__eq__",
+                [](Scalar const& argument1, Scalar const& argument2)
+                { return equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__eq__",
+                [](Scalar const& argument1, Element const& argument2)
+                { return equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+
+            // a != b
+            .def(
+                "__ne__",
+                [](Scalar const& argument1, Scalar const& argument2)
+                { return not_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__ne__",
+                [](Scalar const& argument1, Element const& argument2)
+                { return not_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+
+            // a > b
+            .def(
+                "__gt__",
+                [](Scalar const& argument1, Scalar const& argument2)
+                { return greater_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__gt__",
+                [](Scalar const& argument1, Element const& argument2)
+                { return greater_than<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+
+            // a >= b
+            .def(
+                "__ge__",
+                [](Scalar const& argument1, Scalar const& argument2)
+                { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+            .def(
+                "__ge__",
+                [](Scalar const& argument1, Element const& argument2)
+                { return greater_than_equal_to<BooleanElement>(argument1, argument2); },
+                pybind11::is_operator())
+
+            // a + b, b + a, a += b
             .def(
                 "__add__",
                 [](Scalar const& argument1, Scalar const& argument2) { return add(argument1, argument2); },
                 pybind11::is_operator())
             .def(
                 "__add__",
-                [](Scalar const& argument1, Value const& argument2) { return add(argument1, argument2); },
+                [](Scalar const& argument1, Element const& argument2) { return add(argument1, argument2); },
                 pybind11::is_operator())
+
             .def(
                 "__radd__",
-                [](Scalar const& argument2, Value const& argument1) { return add(argument1, argument2); },
-                pybind11::is_operator());
+                [](Scalar const& argument2, Element const& argument1) { return add(argument1, argument2); },
+                pybind11::is_operator())
 
-        // a - b, b - a, a -= b
-        class_
+            // a - b, b - a, a -= b
             .def(
                 "__sub__",
                 [](Scalar const& argument1, Scalar const& argument2)
@@ -319,103 +182,126 @@ namespace lue::framework {
                 pybind11::is_operator())
             .def(
                 "__sub__",
-                [](Scalar const& argument1, Element const argument2)
-                { return subtract(argument1, argument2); },
-                pybind11::is_operator())
-            .def(
-                "__sub__",
-                [](Scalar const& argument1, Value const& argument2)
-                { return subtract(argument1, argument2); },
-                pybind11::is_operator())
-            .def(
-                "__rsub__",
-                [](Scalar const& argument2, Element const argument1)
+                [](Scalar const& argument1, Element const& argument2)
                 { return subtract(argument1, argument2); },
                 pybind11::is_operator())
+
             .def(
                 "__rsub__",
-                [](Scalar const& argument2, Value const& argument1)
+                [](Scalar const& argument2, Element const& argument1)
                 { return subtract(argument1, argument2); },
                 pybind11::is_operator());
 
+        if constexpr (std::is_signed_v<Element> || std::is_floating_point_v<Element>)
+        {
+            class_
+                // abs(a)
+                .def(
+                    "__abs__",
+                    [](Scalar const& argument) { return abs(argument); },
+                    pybind11::is_operator())
 
-        //         if constexpr (std::is_signed_v<Element> || std::is_floating_point_v<Element>)
-        //         {
-        //             // abs(a)
-        //             class_.def(
-        //                 "__abs__",
-        //                 [](Scalar const& argument) { return lfr::abs<ElementT<Scalar>>(argument); },
-        //                 pybind11::is_operator());
-        //
-        //             // -a
-        //             class_.def(
-        //                 "__neg__",
-        //                 [](Scalar const& argument) { return lfr::negate<ElementT<Scalar>>(argument); },
-        //                 pybind11::is_operator());
-        //         }
-        //
-        //         if constexpr (std::is_integral_v<Element>)
-        //         {
-        //             // a | b
-        //             class_.def(
-        //                 "__or__",
-        //                 [](Scalar const& argument1, Scalar const& argument2)
-        //                 { return lfr::logical_inclusive_or<ElementT<Scalar>>(argument1, argument2); },
-        //                 pybind11::is_operator());
-        //
-        //             // a ^ b
-        //             class_.def(
-        //                 "__xor__",
-        //                 [](Scalar const& argument1, Scalar const& argument2)
-        //                 { return lfr::logical_exclusive_or<ElementT<Scalar>>(argument1, argument2); },
-        //                 pybind11::is_operator());
-        //
-        //             // a & b
-        //             class_.def(
-        //                 "__and__",
-        //                 [](Scalar const& argument1, Scalar const& argument2)
-        //                 { return lfr::logical_and<ElementT<Scalar>>(argument1, argument2); },
-        //                 pybind11::is_operator());
-        //
-        //             // ~a
-        //             class_.def(
-        //                 "__invert__",
-        //                 [](Scalar const& argument) { return lfr::logical_not<ElementT<Scalar>>(argument);
-        //                 }, pybind11::is_operator());
-        //         }
+                // -a
+                .def(
+                    "__neg__",
+                    [](Scalar const& argument) { return negate(argument); },
+                    pybind11::is_operator())
+                ;
+        }
 
-        if constexpr (std::is_floating_point_v<Element>)
+        if constexpr (std::is_integral_v<Element>)
         {
-            // a * b, b * a, a *= b
             class_
+                // a | b
                 .def(
-                    "__mul__",
+                    "__or__",
                     [](Scalar const& argument1, Scalar const& argument2)
-                    { return multiply(argument1, argument2); },
+                    { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
-                    "__mul__",
-                    [](Scalar const& argument1, Element const argument2)
-                    { return multiply(argument1, argument2); },
+                    "__or__",
+                    [](Scalar const& argument1, Element const& argument2)
+                    { return logical_inclusive_or<BooleanElement>(argument1, argument2); },
                     pybind11::is_operator())
+
+                // a ^ b
+                .def(
+                    "__xor__",
+                    [](Scalar const& argument1, Scalar const& argument2)
+                    { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__xor__",
+                    [](Scalar const& argument1, Element const& argument2)
+                    { return logical_exclusive_or<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+
+                // a & b
+                .def(
+                    "__and__",
+                    [](Scalar const& argument1, Scalar const& argument2)
+                    { return logical_and<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__and__",
+                    [](Scalar const& argument1, Element const& argument2)
+                    { return logical_and<BooleanElement>(argument1, argument2); },
+                    pybind11::is_operator())
+
+                // ~a
+                .def(
+                    "__invert__",
+                    [](Scalar const& argument)
+                    { return logical_not<BooleanElement>(argument); },
+                    pybind11::is_operator())
+
+                // a % b
+                .def(
+                    "__mod__",
+                    [](Scalar const& argument1, Scalar const& argument2)
+                    { return modulus(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__mod__",
+                    [](Scalar const& argument1, Element const& argument2)
+                    { return modulus(argument1, argument2); },
+                    pybind11::is_operator())
+
+                .def(
+                    "__rmod__",
+                    [](Scalar const& argument1, Scalar const& argument2)
+                    { return modulus(argument1, argument2); },
+                    pybind11::is_operator())
+                .def(
+                    "__rmod__",
+                    [](Scalar const& argument1, Element const& argument2)
+                    { return modulus(argument1, argument2); },
+                    pybind11::is_operator())
+                ;
+        }
+
+        if constexpr (std::is_floating_point_v<Element>)
+        {
+            class_
+                // a * b, b * a, a *= b
                 .def(
                     "__mul__",
-                    [](Scalar const& argument1, Value const& argument2)
+                    [](Scalar const& argument1, Scalar const& argument2)
                     { return multiply(argument1, argument2); },
                     pybind11::is_operator())
                 .def(
-                    "__rmul__",
-                    [](Scalar const& argument2, Element const argument1)
+                    "__mul__",
+                    [](Scalar const& argument1, Element const& argument2)
                     { return multiply(argument1, argument2); },
                     pybind11::is_operator())
+
                 .def(
                     "__rmul__",
-                    [](Scalar const& argument2, Value const& argument1)
+                    [](Scalar const& argument2, Element const& argument1)
                     { return multiply(argument1, argument2); },
-                    pybind11::is_operator());
+                    pybind11::is_operator())
 
-            // a / b, b / a, a /= b
-            class_
+                // a / b, b / a, a /= b
                 .def(
                     "__truediv__",
                     [](Scalar const& argument1, Scalar const& argument2)
@@ -423,27 +309,17 @@ namespace lue::framework {
                     pybind11::is_operator())
                 .def(
                     "__truediv__",
-                    [](Scalar const& argument1, Element const argument2)
-                    { return divide(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__truediv__",
-                    [](Scalar const& argument1, Value const& argument2)
+                    [](Scalar const& argument1, Element const& argument2)
                     { return divide(argument1, argument2); },
                     pybind11::is_operator())
+
                 .def(
                     "__rtruediv__",
-                    [](Scalar const& argument2, Element const argument1)
+                    [](Scalar const& argument2, Element const& argument1)
                     { return divide(argument1, argument2); },
                     pybind11::is_operator())
-                .def(
-                    "__rtruediv__",
-                    [](Scalar const& argument2, Value const& argument1)
-                    { return divide(argument1, argument2); },
-                    pybind11::is_operator());
 
-            // a ** b, b ** a, a **= b
-            class_
+                // a ** b, b ** a, a **= b
                 .def(
                     "__pow__",
                     [](Scalar const& argument1, Scalar const& argument2)
@@ -451,35 +327,15 @@ namespace lue::framework {
                     pybind11::is_operator())
                 .def(
                     "__pow__",
-                    [](Scalar const& argument1, Element const argument2)
-                    { return pow(argument1, argument2); },
-                    pybind11::is_operator())
-                .def(
-                    "__pow__",
-                    [](Scalar const& argument1, Value const& argument2) { return pow(argument1, argument2); },
+                    [](Scalar const& argument1, Element const& argument2) { return pow(argument1, argument2); },
                     pybind11::is_operator())
+
                 .def(
                     "__rpow__",
-                    [](Scalar const& argument2, Element const argument1)
-                    { return pow(argument1, argument2); },
+                    [](Scalar const& argument2, Element const& argument1) { return pow(argument1, argument2); },
                     pybind11::is_operator())
-                .def(
-                    "__rpow__",
-                    [](Scalar const& argument2, Value const& argument1) { return pow(argument1, argument2); },
-                    pybind11::is_operator());
+                ;
         }
-
-        module.def(
-            "create_scalar",
-            create_scalar_py,
-            R"(
-    Create new scalar, with a value
-
-    :param numpy.dtype dtype: Description of the type of the scalar value
-    :param value: Value to initialize scalar with
-)",
-            "dtype"_a,
-            "value"_a);
     }
 
 }  // namespace lue::framework
diff --git a/source/framework/python/pcraster/__init__.py b/source/framework/python/pcraster/__init__.py
index 9567c4cb6..ed0f305e0 100644
--- a/source/framework/python/pcraster/__init__.py
+++ b/source/framework/python/pcraster/__init__.py
@@ -27,11 +27,27 @@
 arguments to an operator (e.g.: ``a & b``). The regular Python rules are then in effect.
 """
 
+import enum
+import os
+from dataclasses import dataclass
+
 import numpy as np
 
 import lue.framework as lfr
 
 
+@dataclass(kw_only=True)
+class BoundingBox:
+    """
+    Class for keeping track of the bounding box of a raster
+    """
+
+    north: float = 0.0
+    west: float = 0.0
+    south: float = 0.0
+    east: float = 0.0
+
+
 class Configuration(object):
     """
     Class for storing information that is required for LUE to be able to mimic PCRaster
@@ -40,10 +56,18 @@ class Configuration(object):
     the underlying array and the size of the raster cells. This class fulfills the same role.
     """
 
-    def __init__(self, array_shape=None, partition_shape=None, cell_size=None):
+    def __init__(
+        self,
+        *,
+        bounding_box=None,
+        cell_size=None,
+        array_shape=None,
+        partition_shape=None,
+    ):
+        self.bounding_box = bounding_box
+        self.cell_size = cell_size
         self.array_shape = array_shape
         self.partition_shape = partition_shape
-        self.cell_size = cell_size
 
 
 """
@@ -53,13 +77,61 @@ def __init__(self, array_shape=None, partition_shape=None, cell_size=None):
 configuration = Configuration()
 
 
-def setclone(*args):
-    raise NotImplementedError(
-        "Please assign to this module's configuration variable. E.g.:\n"
-        "lue.pcraster.configuration.array_shape = array_shape\n"
-        "lue.pcraster.configuration.partition_shape = partition_shape\n"
-        "lue.pcraster.configuration.cell_size = 10"
-    )
+def setclone(pathname):
+
+    raster_properties = lfr.probe_gdal(pathname)
+    array_shape = raster_properties["shape"]
+
+    geo_transform = raster_properties["geo_transform"]
+    cell_shape = geo_transform[1], -geo_transform[5]
+    assert cell_shape[0] == cell_shape[1]
+    cell_size = cell_shape[0]
+
+    nr_rows, nr_cols = array_shape
+
+    west = geo_transform[0]
+    north = geo_transform[3]
+    east = west + nr_cols * cell_size
+    south = north - nr_rows * cell_size
+    bounding_box = BoundingBox(north=north, west=west, south=south, east=east)
+
+    configuration.bounding_box = bounding_box
+    configuration.cell_size = cell_size
+    configuration.array_shape = array_shape
+
+    partition_shape = os.getenv("LUE_PARTITION_SHAPE", default=None)
+
+    if partition_shape is not None:
+        partition_shape = tuple(int(extent) for extent in partition_shape.split(","))
+
+    configuration.partition_shape = partition_shape
+
+
+def clone():
+    if not hasattr(Configuration, "cellSize"):
+
+        def cellSize(self):
+            return self.cell_size
+
+        setattr(Configuration, "cellSize", cellSize)
+    if not hasattr(Configuration, "nrRows"):
+
+        def nrRows(self):
+            return self.array_shape[0]
+
+        setattr(Configuration, "nrRows", nrRows)
+    if not hasattr(Configuration, "nrCols"):
+
+        def nrCols(self):
+            return self.array_shape[1]
+
+        setattr(Configuration, "nrCols", nrCols)
+
+    return configuration
+
+
+def setglobaloption(option):
+    print(f"lue.pcraster: discarding global option {option}\n")
 
 
 def numpy_scalar_type(expression):
@@ -67,6 +139,10 @@ def numpy_scalar_type(expression):
         getattr(lfr, "PartitionedArray<uint8, 2>"): np.uint8,
         getattr(lfr, "PartitionedArray<int32, 2>"): np.int32,
         getattr(lfr, "PartitionedArray<float32, 2>"): np.float32,
+        getattr(lfr, "Scalar<uint8>"): np.uint8,
+        getattr(lfr, "Scalar<int32>"): np.int32,
+        getattr(lfr, "Scalar<float32>"): np.float32,
+        bool: np.uint8,
         int: np.int32,
         float: np.float32,
         np.uint8: np.uint8,
@@ -78,42 +154,166 @@ def numpy_scalar_type(expression):
     return element_type_by_type[type(expression)]
 
 
-def is_spatial(argument):
-    return isinstance(
-        argument,
-        (
-            getattr(lfr, "PartitionedArray<uint8, 2>"),  # Boolean, LDD
-            getattr(lfr, "PartitionedArray<int32, 2>"),  # Nominal, ordinal
-            getattr(lfr, "PartitionedArray<float32, 2>"),  # Scalar, directional
-        ),
+VALUESCALE = enum.Enum(
+    "VALUESCALE",
+    ["Boolean", "Directional", "Ldd", "Nominal", "Ordinal", "Scalar"],
+)
+(
+    Boolean,
+    Directional,
+    Ldd,
+    Nominal,
+    Ordinal,
+    Scalar,
+) = VALUESCALE
+
+
+# NOTE to_numpy is a synchronization point
+def pcr2numpy(array, no_data_value):
+    return lfr.to_numpy(array, no_data_value)
+
+
+def numpy2pcr(data_type, array, no_data_value):
+    assert (
+        data_type == Ldd
+        and array.dtype == np.uint8
+        or data_type == Scalar
+        and array.dtype == np.float32
+    ), f"{data_type} vs {array.dtype}"
+
+    # TODO LUE get rid of the wait, if possible
+    # Why is the wait needed? OK, don't change the array while LUE is working with
+    # it, but this should not result in a segfault.
+    array = lfr.from_numpy(
+        array,
+        partition_shape=configuration.partition_shape,
+        no_data_value=no_data_value,
     )
+    lfr.wait(array)
+    return array
+
+
+def lue_is_uint8_raster(argument):
+    """
+    Return whether `argument` is an asynchronous LUE uint8 raster
+    """
+    return isinstance(
+        argument, getattr(lfr, "PartitionedArray<uint8, 2>")
+    )  # Boolean, LDD
+
+
+def lue_is_int32_raster(argument):
+    """
+    Return whether `argument` is an asynchronous LUE int32 raster
+    """
+    return isinstance(
+        argument, getattr(lfr, "PartitionedArray<int32, 2>")
+    )  # Nominal, ordinal
+
+
+def lue_is_float32_raster(argument):
+    """
+    Return whether `argument` is an asynchronous LUE float32 raster
+    """
+    return isinstance(
+        argument, getattr(lfr, "PartitionedArray<float32, 2>")
+    )  # Scalar, directional
+
+
+def lue_is_uint8_scalar(argument):
+    """
+    Return whether `argument` is a LUE asynchronous scalar uint8 value
+    """
+    return isinstance(argument, getattr(lfr, "Scalar<uint8>"))  # Boolean
+
+
+def lue_is_int32_scalar(argument):
+    """
+    Return whether `argument` is a LUE asynchronous scalar int32 value
+    """
+    return isinstance(argument, getattr(lfr, "Scalar<int32>"))  # Nominal, ordinal
+
+
+def lue_is_float32_scalar(argument):
+    """
+    Return whether `argument` is a LUE asynchronous scalar float32 value
+    """
+    return isinstance(argument, getattr(lfr, "Scalar<float32>"))  # Scalar
 
 
 def is_boolean(expression):
-    return isinstance(expression, (getattr(lfr, "PartitionedArray<uint8, 2>")))
+    """
+    Return whether `expression`'s value scale is boolean
+    """
+    return lue_is_uint8_raster(expression) or lue_is_uint8_scalar(expression)
 
 
 def is_ldd(expression):
-    return isinstance(expression, (getattr(lfr, "PartitionedArray<uint8, 2>")))
+    """
+    Return whether `expression`'s value scale is LDD
+    """
+    return lue_is_uint8_raster(expression)  # or lue_is_uint8_scalar(expression)
 
 
 def is_nominal(expression):
-    return isinstance(expression, (getattr(lfr, "PartitionedArray<int32, 2>")))
+    """
+    Return whether `expression`'s value scale is nominal
+    """
+    return lue_is_int32_raster(expression) or lue_is_int32_scalar(expression)
 
 
 def is_ordinal(expression):
-    return isinstance(expression, (getattr(lfr, "PartitionedArray<int32, 2>")))
+    """
+    Return whether `expression`'s value scale is ordinal
+    """
+    return lue_is_int32_raster(expression) or lue_is_int32_scalar(expression)
 
 
 def is_scalar(expression):
-    return isinstance(expression, (getattr(lfr, "PartitionedArray<float32, 2>")))
+    """
+    Return whether `expression`'s value scale is scalar
+    """
+    return lue_is_float32_raster(expression) or lue_is_float32_scalar(expression)
 
 
 def is_directional(expression):
-    return isinstance(expression, (getattr(lfr, "PartitionedArray<float32, 2>")))
+    """
+    Return whether `expression`'s value scale is directional
+    """
+    return lue_is_float32_raster(expression)  # or lue_is_float32_scalar(expression)
+
+
+def is_spatial(argument):
+    """
+    Return whether `argument` is spatial
+
+    A spatial is represented by a LUE asynchronous raster
+    """
+    return (
+        lue_is_uint8_raster(argument)
+        or lue_is_int32_raster(argument)
+        or lue_is_float32_raster(argument)
+    )
 
 
 def is_non_spatial(argument):
+    """
+    Return whether `argument` is non-spatial
+
+    A non-spatial is represented by a LUE asynchronous scalar
+    """
+    return (
+        lue_is_uint8_scalar(argument)
+        or lue_is_int32_scalar(argument)
+        or lue_is_float32_scalar(argument)
+        # or lue_is_value(argument)
+    )
+
+
+def lue_is_value(argument):
+    """
+    Return whether `argument` is a numeric value
+    """
     return isinstance(
         argument,
         (
@@ -131,7 +331,8 @@ def is_non_spatial(argument):
 
 
 def non_spatial_to_spatial(fill_value, template=None):
-    scalar_type = numpy_scalar_type(fill_value)
+    if lue_is_value(fill_value):
+        fill_value = lfr.create_scalar(numpy_scalar_type(fill_value), fill_value)
 
     if template is None:
         array_shape = configuration.array_shape
@@ -140,9 +341,41 @@ def non_spatial_to_spatial(fill_value, template=None):
 
     partition_shape = configuration.partition_shape
 
-    return lfr.create_array(
-        array_shape, scalar_type, fill_value, partition_shape=partition_shape
-    )
+    return lfr.create_array(array_shape, fill_value, partition_shape=partition_shape)
+
+
+def harmonize_types(expression1, expression2):
+
+    # If both arguments are values, fix the element type of the first one
+    # TODO Could determine the largest one and use that
+    if lue_is_value(expression1) and lue_is_value(expression2):
+        expression1 = lfr.create_scalar(numpy_scalar_type(expression1), expression1)
+
+    # If one of the arguments is a value, fix its element type to the element type of the other argument
+    to_scalar_1 = (
+        is_spatial(expression2) or is_non_spatial(expression2)
+    ) and lue_is_value(expression1)
+    to_scalar_2 = (
+        is_spatial(expression1) or is_non_spatial(expression1)
+    ) and lue_is_value(expression2)
+
+    if to_scalar_2:
+        expression2 = lfr.create_scalar(
+            np.dtype(numpy_scalar_type(expression1)), expression2
+        )
+
+    if to_scalar_1:
+        expression1 = lfr.create_scalar(
+            np.dtype(numpy_scalar_type(expression2)), expression1
+        )
+
+    assert not lue_is_value(expression1), expression1
+    assert not lue_is_value(expression2), expression2
+    assert (
+        expression1.dtype == expression2.dtype
+    ), f"{expression1.dtype} != {expression2.dtype}"
+
+    return expression1, expression2
 
 
 def read_if_necessary(*args) -> tuple:
@@ -159,7 +392,7 @@ def translate_window_length(pcraster_window_length):
     - corresponds with an odd number of cells (odd_number * cell_size)
     - is positive
     """
-    assert is_non_spatial(pcraster_window_length), pcraster_window_length
+    assert lue_is_value(pcraster_window_length), pcraster_window_length
 
     window_length = pcraster_window_length / configuration.cell_size
 
@@ -169,6 +402,8 @@ def translate_window_length(pcraster_window_length):
 
     window_length = int(window_length)
 
+    # assert window_length > 1, "Window lengths should be larger than one cell size"
+
     return window_length
 
 
@@ -177,8 +412,10 @@ def readmap(pathname):
 
 
 def report(expression, pathname):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=expression)
+    # TODO Right? Otherwise, don't create spatials here. Use another format to write scalars to.
+    assert is_spatial(expression)
+    # if is_non_spatial(expression):
+    #     expression = non_spatial_to_spatial(fill_value=expression)
 
     lfr.to_gdal(expression, pathname)
 
@@ -188,103 +425,101 @@ def div(expression1, expression2):
 
 
 def abs(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
+    # TODO if is_non_spatial(expression):
+    # TODO     expression = non_spatial_to_spatial(fill_value=np.float32(expression))
 
     return lfr.abs(expression)
 
 
-def accucapacityflux(ldd, material, transportcapacity):
+def accucapacityflux(flow_direction, material, transportcapacity):
     raise NotImplementedError("accucapacityflux")
 
 
-def accucapacitystate(ldd, material, transportcapacity):
+def accucapacitystate(flow_direction, material, transportcapacity):
     raise NotImplementedError("accucapacitystate")
 
 
-def accuflux(ldd, material):
-    raise NotImplementedError("accuflux")
+def accuflux(flow_direction, material):
+    return accuthreshold(flow_direction, material, 0)[1]
 
 
-def accufraction(ldd, material, transportcapacity):
-    ldd, material, transportcapacity = read_if_necessary(
-        ldd, material, transportcapacity
-    )
+def accufraction(flow_direction, material, transportcapacity):
+    flow_direction = ldd(flow_direction)
+    material = scalar(material)
+    transportcapacity = scalar(transportcapacity)
 
-    assert is_spatial(ldd), type(ldd)
+    assert is_spatial(flow_direction), type(flow_direction)
 
     if is_non_spatial(material):
-        material = non_spatial_to_spatial(fill_value=np.float32(material), template=ldd)
+        # TODO Support non-spatial material
+        material = non_spatial_to_spatial(fill_value=material, template=flow_direction)
 
     if is_non_spatial(transportcapacity):
-        transportcapacity = non_spatial_to_spatial(
-            fill_value=np.float32(transportcapacity), template=ldd
-        )
+        # TODO Support non-spatial transportcapacity
+        transportcapacity = non_spatial_to_spatial(fill_value=transportcapacity, template=flow_direction)
 
-    return lfr.accu_fraction(ldd, material, transportcapacity)
+    return lfr.accu_fraction(flow_direction, material, transportcapacity)
 
 
-def accufractionflux(ldd, material, transportcapacity):
-    return accufraction(ldd, material, transportcapacity)[0]
+def accufractionflux(flow_direction, material, transportcapacity):
+    return accufraction(flow_direction, material, transportcapacity)[0]
 
 
-def accufractionstate(ldd, material, transportcapacity):
-    return accufraction(ldd, material, transportcapacity)[1]
+def accufractionstate(flow_direction, material, transportcapacity):
+    return accufraction(flow_direction, material, transportcapacity)[1]
 
 
-def accuthreshold(ldd, material, threshold):
-    if isinstance(ldd, str):
-        ldd = readmap(ldd)
+def accuthreshold(flow_direction, material, threshold):
+    flow_direction = ldd(flow_direction)
+    material = scalar(material)
+    threshold = scalar(threshold)
 
-    assert is_spatial(ldd), type(ldd)
+    assert is_spatial(flow_direction), type(flow_direction)
 
     if is_non_spatial(material):
-        material = non_spatial_to_spatial(fill_value=np.float32(material), template=ldd)
+        # TODO Support non-spatial material
+        material = non_spatial_to_spatial(fill_value=material, template=flow_direction)
 
     if is_non_spatial(threshold):
-        threshold = non_spatial_to_spatial(
-            fill_value=np.float32(threshold), template=ldd
-        )
+        # TODO Support non-spatial threshold
+        threshold = non_spatial_to_spatial(fill_value=threshold, template=flow_direction)
 
-    return lfr.accu_threshold3(ldd, material, threshold)
+    return lfr.accu_threshold3(flow_direction, material, threshold)
 
 
-def accuthresholdflux(ldd, material, threshold):
-    return accuthreshold(ldd, material, threshold)[0]
+def accuthresholdflux(flow_direction, material, threshold):
+    return accuthreshold(flow_direction, material, threshold)[0]
 
 
-def accuthresholdstate(ldd, material, threshold):
-    return accuthreshold(ldd, material, threshold)[1]
+def accuthresholdstate(flow_direction, material, threshold):
+    return accuthreshold(flow_direction, material, threshold)[1]
 
 
-def accutriggerflux(ldd, material, transporttrigger):
+def accutriggerflux(flow_direction, material, transporttrigger):
     raise NotImplementedError("accutriggerflux")
 
 
-def accutriggerstate(ldd, material, transporttrigger):
+def accutriggerstate(flow_direction, material, transporttrigger):
     raise NotImplementedError("accutriggerstate")
 
 
-def accutraveltimeflux(ldd, material, transporttraveltime):
+def accutraveltimeflux(flow_direction, material, transporttraveltime):
     raise NotImplementedError("accutraveltimeflux")
 
 
-def accutraveltimestate(ldd, material, transporttraveltime):
+def accutraveltimestate(flow_direction, material, transporttraveltime):
     raise NotImplementedError("accutraveltimestate")
 
 
-def accutraveltimefractionflux(ldd, material, transporttraveltime):
+def accutraveltimefractionflux(flow_direction, material, transporttraveltime):
     raise NotImplementedError("accutraveltimefractionflux")
 
 
-def accutraveltimefractionstate(ldd, material, transporttraveltime):
+def accutraveltimefractionstate(flow_direction, material, transporttraveltime):
     raise NotImplementedError("accutraveltimefractionstate")
 
 
 def acos(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.acos(expression)
 
 
@@ -357,9 +592,6 @@ def areauniform(areaclass):
 
 
 def asin(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.asin(expression)
 
 
@@ -368,22 +600,23 @@ def aspect(dem):
 
 
 def atan(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.atan(expression)
 
 
 def boolean(expression):
     expression = read_if_necessary(expression)[0]
 
-    if is_spatial(expression):
-        if is_boolean(expression):
-            return expression
-        else:
-            return lfr.where(expression != 0, np.uint8(1), np.uint8(0))
+    if is_boolean(expression):
+        return expression
+    elif is_spatial(expression):
+        # TODO where doesn't handle dypes correctly
+        # https://github.com/computationalgeography/lue/issues/687
+        # return lfr.where(expression != 0, np.uint8(1), np.uint8(0))
+        return lfr.cast(lfr.where(expression != 0, np.uint8(1), np.uint8(0)), np.uint8)
     elif is_non_spatial(expression):
-        return np.uint8(expression)
+        return expression != 0
+    elif lue_is_value(expression):
+        return lfr.create_scalar(np.uint8, expression != 0)
 
     raise RuntimeError("Unsupported argument: {}".format(expression))
 
@@ -392,8 +625,10 @@ def catchment(*args):
     raise NotImplementedError("catchment")
 
 
-def catchmenttotal(*args):
-    raise NotImplementedError("catchmenttotal")
+def catchmenttotal(amount, flow_direction):
+    # TODO This assumes accuthreshold can handle negative values
+    # https://github.com/computationalgeography/lue/issues/673
+    return accuthreshold(flow_direction, amount, 0)[1]
 
 
 def cellarea(*args):
@@ -404,20 +639,27 @@ def celllength(*args):
     return configuration.cell_size
 
 
+def cellvalue(*args):
+    raise NotImplementedError("cellvalue")
+
+
 def clump(expression):
     return lfr.clump(expression, lfr.Connectivity.diagonal)
 
 
 def cos(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.cos(expression)
 
 
 def cover(expression1, expression2, *expressions):
-    if not is_spatial(expression1):
-        expression1 = non_spatial_to_spatial(fill_value=expression1)
+
+    expression1, expression2 = harmonize_types(expression1, expression2)
+
+    # TODO Should where support scalar condition?
+    if is_non_spatial(expression1):
+        expression1 = non_spatial_to_spatial(expression1)
+
+    assert is_spatial(expression1), expression1
 
     result = lfr.where(lfr.valid(expression1), expression1, expression2)
 
@@ -428,11 +670,7 @@ def cover(expression1, expression2, *expressions):
 
 
 def defined(expression):
-    if isinstance(expression, str):
-        expression = readmap(expression)
-
-    if not is_spatial(expression):
-        return non_spatial_to_spatial(fill_value=np.uint8(1))
+    expression = read_if_necessary(expression)[0]
 
     return lfr.valid(expression)
 
@@ -440,29 +678,26 @@ def defined(expression):
 def directional(expression):
     expression = read_if_necessary(expression)[0]
 
-    if is_spatial(expression):
-        if is_directional(expression):
-            return expression
-        else:
-            return lfr.cast(expression, np.float32)
-    elif is_non_spatial(expression):
-        return np.float32(expression)
+    if is_directional(expression):
+        return expression
+    elif is_spatial(expression) or is_non_spatial(expression):
+        return lfr.cast(expression, np.float32)
+    elif lue_is_value(expression):
+        return lfr.create_scalar(np.float32, expression)
 
     raise RuntimeError("Unsupported argument: {}".format(expression))
 
 
-def downstream(ldd, expression):
-    if isinstance(ldd, str):
-        ldd = readmap(ldd)
+def downstream(flow_direction, expression):
+    flow_direction = read_if_necessary(flow_direction)[0]
 
-    return lfr.downstream(ldd, expression)
+    return lfr.downstream(flow_direction, expression)
 
 
-def downstreamdist(ldd):
-    if isinstance(ldd, str):
-        ldd = readmap(ldd)
+def downstreamdist(flow_direction):
+    flow_direction = read_if_necessary(flow_direction)[0]
 
-    return lfr.downstream_distance(ldd, configuration.cell_size, np.float32)
+    return lfr.downstream_distance(flow_direction, configuration.cell_size, np.float32)
 
 
 def dynwavestate(*args):
@@ -490,9 +725,6 @@ def eq(expression1, expression2):
 
 
 def exp(power):
-    if is_non_spatial(power):
-        power = non_spatial_to_spatial(fill_value=np.float32(power))
-
     return lfr.exp(power)
 
 
@@ -521,36 +753,35 @@ def idiv(*args):
 
 
 def ifthen(condition, expression):
-    if is_non_spatial(condition):
-        assert condition < np.iinfo(np.uint8).max, condition
-        condition = non_spatial_to_spatial(fill_value=np.uint8(condition))
-
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=expression)
+    # TODO Should where support scalar condition?
+    if not is_spatial(condition):
+        if lue_is_value(condition):
+            condition = lfr.create_scalar(np.uint8, condition)
+        condition = non_spatial_to_spatial(condition)
+
+    if lue_is_value(expression):
+        expression = lfr.create_scalar(
+            np.dtype(numpy_scalar_type(expression)), expression
+        )
 
     return lfr.where(condition, expression)
 
 
 def ifthenelse(condition, expression1, expression2):
-    if is_non_spatial(condition):
-        assert condition < np.iinfo(np.uint8).max, condition
-        condition = non_spatial_to_spatial(fill_value=np.uint8(condition))
-
-    if is_spatial(expression1):
-        expression_dtype = expression1.dtype
-    elif is_spatial(expression2):
-        expression_dtype = expression2.dtype
-    else:
-        expression_dtype = np.dtype(numpy_scalar_type(expression1))
-
-    if is_non_spatial(expression1):
-        expression1 = non_spatial_to_spatial(
-            fill_value=expression_dtype.type(expression1)
+    # TODO Should where support scalar condition?
+    if not is_spatial(condition):
+        if lue_is_value(condition):
+            condition = lfr.create_scalar(np.uint8, condition)
+        condition = non_spatial_to_spatial(condition)
+
+    if lue_is_value(expression1):
+        expression1 = lfr.create_scalar(
+            np.dtype(numpy_scalar_type(expression1)), expression1
         )
 
-    if is_non_spatial(expression2):
-        expression2 = non_spatial_to_spatial(
-            fill_value=expression_dtype.type(expression2)
+    if lue_is_value(expression2):
+        expression2 = lfr.create_scalar(
+            np.dtype(numpy_scalar_type(expression2)), expression2
         )
 
     return lfr.where(condition, expression1, expression2)
@@ -564,8 +795,19 @@ def inversedistance(*args):
     raise NotImplementedError("inversedistance")
 
 
-def kinematic(ldd, Qold, q, alpha, beta, nrTimeSlices, dT, dX):
-    return lfr.kinematic_wave(ldd, Qold, q, alpha, beta, dT, dX)
+def kinematic(flow_direction, Qold, q, alpha, beta, nrTimeSlices, dT, dX):
+    # TODO kinematic wave should also support Scalar values
+
+    # flow_direction = ldd(flow_direction)
+    # Qold = scalar(Qold)
+    # q = scalar(q)
+    # alpha = scalar(alpha)
+    # beta = scalar(beta)
+    # nrTimeSlices = scalar(nrTimeSlices)
+    # dT = scalar(dT)
+    # dX = scalar(dX)
+
+    return lfr.kinematic_wave(flow_direction, Qold, q, alpha, beta, dT, dX)
 
 
 def kinwavestate(*args):
@@ -579,13 +821,12 @@ def kinwaveflux(*args):
 def ldd(expression):
     expression = read_if_necessary(expression)[0]
 
-    if is_spatial(expression):
-        if is_ldd(expression):
-            return expression
-        else:
-            return lfr.cast(expression, np.uint8)
-    elif is_non_spatial(expression):
-        return np.uint8(expression)
+    if is_ldd(expression):
+        return expression
+    elif is_spatial(expression) or is_non_spatial(expression):
+        return lfr.cast(expression, np.uint8)
+    elif lue_is_value(expression):
+        return lfr.create_scalar(np.uint8, expression)
 
     raise RuntimeError("Unsupported argument: {}".format(expression))
 
@@ -599,7 +840,7 @@ def lddcreatedem(*args):
 
 
 def ldddist(*args):
-    raise NotImplementedError("lddcreatedist")
+    raise NotImplementedError("ldddist")
 
 
 def lddmask(*args):
@@ -615,16 +856,10 @@ def le(expression1, expression2):
 
 
 def ln(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.log(expression)
 
 
 def log10(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.log10(expression)
 
 
@@ -657,7 +892,7 @@ def mapminimum(expression):
 
 
 def mapnormal():
-    return normal(1)
+    raise NotImplementedError("mapnormal")
 
 
 def maptotal(expression):
@@ -665,7 +900,7 @@ def maptotal(expression):
 
 
 def mapuniform():
-    return uniform(1)
+    raise NotImplementedError("mapuniform")
 
 
 def markwhilesumle(*args):
@@ -679,45 +914,35 @@ def markwhilesumge(*args):
 def max(*args):
     if len(args) == 1:
         return args[0]
+    elif len(args) > 2:
+        return max(args[0], max(*args[1:]))
+    else:
+        expression1, expression2 = harmonize_types(*args)
 
-    if len(args) != 2:
-        raise NotImplementedError("max with more than 2 arguments")
-
-    expression1, expression2 = args
-    type1 = numpy_scalar_type(expression1)
-    type2 = numpy_scalar_type(expression2)
-
-    if is_non_spatial(expression1):
-        expression1 = non_spatial_to_spatial(fill_value=type2(expression1))
-
-    if is_non_spatial(expression2):
-        expression2 = non_spatial_to_spatial(fill_value=type1(expression2))
+        # TODO Support non-spatial condition(?)
+        condition = expression1 >= expression2
 
-    condition = lfr.greater_than_equal_to(expression1, expression2)
+        if is_non_spatial(condition):
+            condition = non_spatial_to_spatial(condition)
 
-    return lfr.where(condition, expression1, expression2)
+        return lfr.where(condition, expression1, expression2)
 
 
 def min(*args):
     if len(args) == 1:
         return args[0]
+    elif len(args) > 2:
+        return min(args[0], min(*args[1:]))
+    else:
+        expression1, expression2 = harmonize_types(*args)
 
-    if len(args) != 2:
-        raise NotImplementedError("min with more than 2 arguments")
-
-    expression1, expression2 = args
-    type1 = numpy_scalar_type(expression1)
-    type2 = numpy_scalar_type(expression2)
-
-    if is_non_spatial(expression1):
-        expression1 = non_spatial_to_spatial(fill_value=type2(expression1))
-
-    if is_non_spatial(expression2):
-        expression2 = non_spatial_to_spatial(fill_value=type1(expression2))
+        # TODO Support non-spatial condition(?)
+        condition = expression1 <= expression2
 
-    condition = lfr.less_than_equal_to(expression1, expression2)
+        if is_non_spatial(condition) or lue_is_value(condition):
+            condition = non_spatial_to_spatial(condition)
 
-    return lfr.where(condition, expression1, expression2)
+        return lfr.where(condition, expression1, expression2)
 
 
 def mod(*args):
@@ -735,20 +960,23 @@ def nodirection(*args):
 def nominal(expression):
     expression = read_if_necessary(expression)[0]
 
-    if is_spatial(expression):
-        if is_nominal(expression):
-            return expression
-        else:
-            return lfr.cast(expression, np.int32)
-    elif is_non_spatial(expression):
-        return np.int32(expression)
+    if is_nominal(expression):
+        return expression
+    elif is_spatial(expression) or is_non_spatial(expression):
+        return lfr.cast(expression, np.int32)
+    elif lue_is_value(expression):
+        return lfr.create_scalar(np.int32, expression)
 
     raise RuntimeError("Unsupported argument: {}".format(expression))
 
 
 def normal(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.uint8(expression))
+    # TODO Support scalar input(?)
+    if is_non_spatial(expression) or lue_is_value(expression):
+        if lue_is_value(expression):
+            expression = lfr.create_scalar(np.uint8, expression)
+
+        expression = non_spatial_to_spatial(fill_value=expression)
 
     return lfr.where(expression, lfr.normal(expression, np.float32, 0, 1))
 
@@ -768,13 +996,12 @@ def order(*args):
 def ordinal(expression):
     expression = read_if_necessary(expression)[0]
 
-    if is_spatial(expression):
-        if is_ordinal(expression):
-            return expression
-        else:
-            return lfr.cast(expression, np.int32)
-    elif is_non_spatial(expression):
-        return np.int32(expression)
+    if is_ordinal(expression):
+        return expression
+    elif is_spatial(expression) or is_non_spatial(expression):
+        return lfr.cast(expression, np.int32)
+    elif lue_is_value(expression):
+        return lfr.create_scalar(np.int32, expression)
 
     raise RuntimeError("Unsupported argument: {}".format(expression))
 
@@ -800,36 +1027,26 @@ def profcurv(*args):
 
 
 def rounddown(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.floor(expression)
 
 
 def roundoff(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.round(expression)
 
 
 def roundup(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.ceil(expression)
 
 
 def scalar(expression):
     expression = read_if_necessary(expression)[0]
 
-    if is_spatial(expression):
-        if is_scalar(expression):
-            return expression
-        else:
-            return lfr.cast(expression, np.float32)
-    elif is_non_spatial(expression):
-        return np.float32(expression)
+    if is_scalar(expression):
+        return expression
+    elif is_spatial(expression) or is_non_spatial(expression):
+        return lfr.cast(expression, np.float32)
+    elif lue_is_value(expression):
+        return lfr.create_scalar(np.float32, expression)
 
     raise RuntimeError("Unsupported argument: {}".format(expression))
 
@@ -843,9 +1060,6 @@ def shift0(*args):
 
 
 def sin(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.sin(expression)
 
 
@@ -889,16 +1103,10 @@ def spreadzone(*args):
 
 
 def sqr(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return expression * expression
 
 
 def sqrt(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.log(expression)
 
 
@@ -915,9 +1123,6 @@ def succ(*args):
 
 
 def tan(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.float32(expression))
-
     return lfr.tan(expression)
 
 
@@ -937,8 +1142,18 @@ def timeinputordinal(*args):
     raise NotImplementedError("timeinputordinal")
 
 
-def timeinputscalar(*args):
-    raise NotImplementedError("timeinputscalar")
+def timeinputscalar(timeseries_pathname: str, id_expression, timestep: int):
+    lines = open(timeseries_pathname).readlines()
+    nr_rows_to_skip = int(lines[1]) + 2
+    lines = lines[nr_rows_to_skip:]
+    line = lines[timestep - 1]
+    values = line.split()[1:]
+    lookup_table = {id_ + 1: float(values[id_]) for id_ in range(len(values))}
+
+    if isinstance(id_expression, int):
+        return lookup_table[id_expression]
+    else:
+        return lfr.reclassify(id_expression, lookup_table, dtype=np.float32)
 
 
 def timeinputdirectional(*args):
@@ -974,26 +1189,35 @@ def transient(*args):
 
 
 def uniform(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.uint8(expression))
+    # TODO Support scalar input(?)
+    if is_non_spatial(expression) or lue_is_value(expression):
+        if lue_is_value(expression):
+            expression = lfr.create_scalar(np.uint8, expression)
+
+        expression = non_spatial_to_spatial(fill_value=expression)
 
     return lfr.where(expression, lfr.uniform(expression, np.float32, 0, 1))
 
 
 def uniqueid(expression):
-    if is_non_spatial(expression):
-        expression = non_spatial_to_spatial(fill_value=np.uint8(expression))
+    # TODO Support scalar input(?)
+    if is_non_spatial(expression) or lue_is_value(expression):
+        if lue_is_value(expression):
+            expression = lfr.create_scalar(np.uint8, expression)
+
+        expression = non_spatial_to_spatial(fill_value=expression)
 
     return lfr.where(expression, lfr.unique_id(expression, dtype=np.int32), 0)
 
 
-def upstream(ldd, material):
-    ldd, material = read_if_necessary(ldd, material)
+def upstream(flow_direction, material):
+    flow_direction, material = read_if_necessary(flow_direction, material)
 
-    if is_non_spatial(material):
+    if is_non_spatial(material) or lue_is_value(material):
+        # TODO Support non-spatial material
         material = non_spatial_to_spatial(fill_value=np.float32(material))
 
-    return lfr.upstream(ldd, material)
+    return lfr.upstream(flow_direction, material)
 
 
 def view(*args):
@@ -1063,13 +1287,38 @@ def windowtotal(expression, window_length):
     return lfr.focal_sum(expression, window)
 
 
-def xcoordinate(*args):
-    raise NotImplementedError("xcoordinate")
+def xcoordinate(expression):
+    # TODO Support non-spatial input?
+    if is_non_spatial(expression) or lue_is_value(expression):
+        expression = non_spatial_to_spatial(expression)
+
+    return (
+        configuration.bounding_box.west
+        + 0.5 * configuration.cell_size
+        + (
+            lfr.cast(lfr.cell_index(defined(expression), 1), np.float32)
+            * configuration.cell_size
+        )
+    )
 
 
 def pcrxor(expression1, expression2):
     return expression1 ^ expression2
 
 
-def ycoordinate(*args):
-    raise NotImplementedError("ycoordinate")
+def ycoordinate(expression):
+    # TODO Support non-spatial input?
+    if is_non_spatial(expression) or lue_is_value(expression):
+        expression = non_spatial_to_spatial(expression)
+
+    return (
+        configuration.bounding_box.south
+        + 0.5 * configuration.cell_size
+        + (
+            (
+                -lfr.cast(lfr.cell_index(defined(expression), 0), np.float32)
+                + (configuration.array_shape[0] - 1)
+            )
+            * configuration.cell_size
+        )
+    )
diff --git a/source/framework/python/pcraster/framework/__init__.py b/source/framework/python/pcraster/framework/__init__.py
index 7db18ebfc..ea16593c9 100644
--- a/source/framework/python/pcraster/framework/__init__.py
+++ b/source/framework/python/pcraster/framework/__init__.py
@@ -99,10 +99,17 @@ class DynamicFramework(lfr.DynamicModelRunner):
 
     def __init__(self, model, last_time_step=0, first_time_step=1):
         super().__init__(model, last_time_step, first_time_step)
+        self.quiet = False
 
     def setQuiet(self, quiet=True):
-        # Skip for now
-        pass
+        self.quiet = quiet
+
+    def run(self, *, rate_limit=0):
+        if self.quiet:
+            progressor = lfr.SilentProgressor()
+        else:
+            progressor = lfr.DefaultProgressor()
+        super().run(progressor=progressor, rate_limit=rate_limit)
 
 
 class MonteCarloModel(lfr.MonteCarloModel):
diff --git a/source/framework/python/src/core/create_partitioned_array.cpp b/source/framework/python/src/core/create_partitioned_array.cpp
index ccb6aaf27..f5bc169e4 100644
--- a/source/framework/python/src/core/create_partitioned_array.cpp
+++ b/source/framework/python/src/core/create_partitioned_array.cpp
@@ -23,7 +23,11 @@ namespace lue::framework {
                 Policies{},
                 array_shape,
                 partition_shape,
-                Functor{hpx::make_ready_future<Element>(pybind11::cast<Element>(fill_value))});
+                Functor{
+                    pybind11::isinstance<Scalar<Element>>(fill_value) ?
+                    pybind11::cast<Scalar<Element>>(fill_value).future() :
+                    hpx::make_ready_future<Element>(pybind11::cast<Element>(fill_value))
+                });
         }
 
 
@@ -193,6 +197,192 @@ namespace lue::framework {
             return result;
         }
 
+
+        template<typename Element, Rank rank>
+        auto create_array2(
+            ShapeT<PartitionedArray<Element, rank>> const& array_shape,
+            ShapeT<PartitionedArray<Element, rank>> const& partition_shape,
+            pybind11::object const& fill_value) -> PartitionedArray<Element, rank>
+        {
+            using Policies = lue::policy::create_partitioned_array::DefaultValuePolicies<Element>;
+            using Functor = lue::InstantiateFilled<Element, rank>;
+
+            return lue::create_partitioned_array(
+                Policies{},
+                array_shape,
+                partition_shape,
+                Functor{pybind11::cast<Scalar<Element>>(fill_value).future()});
+        }
+
+
+        template<Rank rank>
+        auto create_array2(
+            StaticShape<rank> const& array_shape,
+            StaticShape<rank> const& partition_shape,
+            pybind11::dtype const& dtype,
+            pybind11::object const& fill_value) -> pybind11::object
+        {
+            // Switch on dtype and call a function that returns an array of the
+            // right value type
+            auto const kind = dtype.kind();
+            auto const size = dtype.itemsize();  // bytes
+            pybind11::object result;
+
+            switch (kind)
+            {
+                case 'i':
+                {
+                    // Signed integer
+                    switch (size)
+                    {
+                        case 4:
+                        {
+                            result = pybind11::cast(
+                                create_array2<std::int32_t, rank>(array_shape, partition_shape, fill_value));
+                            break;
+                        }
+                        case 8:
+                        {
+                            result = pybind11::cast(
+                                create_array2<std::int64_t, rank>(array_shape, partition_shape, fill_value));
+                            break;
+                        }
+                    }
+
+                    break;
+                }
+                case 'u':
+                {
+                    // Unsigned integer
+                    switch (size)
+                    {
+                        case 1:
+                        {
+                            result = pybind11::cast(
+                                create_array2<std::uint8_t, rank>(array_shape, partition_shape, fill_value));
+                            break;
+                        }
+                        case 4:
+                        {
+                            result = pybind11::cast(
+                                create_array2<std::uint32_t, rank>(array_shape, partition_shape, fill_value));
+                            break;
+                        }
+                        case 8:
+                        {
+                            result = pybind11::cast(
+                                create_array2<std::uint64_t, rank>(array_shape, partition_shape, fill_value));
+                            break;
+                        }
+                    }
+
+                    break;
+                }
+                case 'f':
+                {
+                    // Floating-point
+                    switch (size)
+                    {
+                        case 4:
+                        {
+                            result = pybind11::cast(
+                                create_array2<float, rank>(array_shape, partition_shape, fill_value));
+                            break;
+                        }
+                        case 8:
+                        {
+                            result = pybind11::cast(
+                                create_array2<double, rank>(array_shape, partition_shape, fill_value));
+                            break;
+                        }
+                    }
+
+                    break;
+                }
+            }
+
+            if (!result)
+            {
+                throw std::runtime_error(fmt::format("Unsupported dtype (kind={}, itemsize={})", kind, size));
+            }
+
+            return result;
+        }
+
+
+        auto create_array2_py(
+            pybind11::tuple const& array_shape,
+            pybind11::object const& fill_value,
+            std::optional<pybind11::tuple> const& partition_shape) -> pybind11::object
+        {
+            DynamicShape const dynamic_array_shape{tuple_to_shape(array_shape)};
+
+            // The fill_value should be a Scalar<Element> instance. Ask for its dtype information.
+            pybind11::dtype const dtype{fill_value.attr("dtype")};
+
+            Rank const rank{dynamic_array_shape.size()};
+
+            if (partition_shape)
+            {
+                DynamicShape const dynamic_partition_shape{tuple_to_shape(*partition_shape)};
+
+                if (dynamic_array_shape.size() != dynamic_partition_shape.size())
+                {
+                    throw std::runtime_error(fmt::format(
+                        "Rank of array shape and partition shape must be equal ({} != {})",
+                        dynamic_array_shape.size(),
+                        dynamic_partition_shape.size()));
+                }
+            }
+
+            pybind11::object result;
+
+            if (rank == 1)
+            {
+                StaticShape<1> const static_array_shape{
+                    dynamic_shape_to_static_shape<1>(dynamic_array_shape)};
+                StaticShape<1> static_partition_shape{};
+
+                if (partition_shape)
+                {
+                    static_partition_shape =
+                        dynamic_shape_to_static_shape<1>(tuple_to_shape(*partition_shape));
+                }
+                else
+                {
+                    static_partition_shape = default_partition_shape(static_array_shape);
+                }
+
+                result = create_array2<1>(static_array_shape, static_partition_shape, dtype, fill_value);
+            }
+            else if (rank == 2)
+            {
+                StaticShape<2> const static_array_shape{
+                    dynamic_shape_to_static_shape<2>(dynamic_array_shape)};
+                StaticShape<2> static_partition_shape{};
+
+                if (partition_shape)
+                {
+                    static_partition_shape =
+                        dynamic_shape_to_static_shape<2>(tuple_to_shape(*partition_shape));
+                }
+                else
+                {
+                    static_partition_shape = default_partition_shape(static_array_shape);
+                }
+
+                result = create_array2<2>(static_array_shape, static_partition_shape, dtype, fill_value);
+            }
+            else
+            {
+                throw std::runtime_error("Currently only arrays of rank 1 or 2 are supported");
+            }
+
+            lue_hpx_assert(result);
+
+            return result;
+        }
+
     }  // Anonymous namespace
 
 
@@ -221,6 +411,28 @@ namespace lue::framework {
             pybind11::kw_only(),
             "partition_shape"_a = std::optional<pybind11::tuple>{},
             pybind11::return_value_policy::move);
+
+        module.def(
+            "create_array",
+            create_array2_py,
+            R"(
+    Create new array, filled with a value
+
+    :param tuple array_shape: Shape of the array
+    :param fill_value: Value to fill array with
+    :param tuple partition_shape: Shape of the array partitions. When not
+        passed in, a default shape will be used which might not result in the
+        best performance and scalability.
+    :rtype: PartitionedArray specialization
+
+    The type of the array returned depends on the rank of the array and
+    the type of the array elements.
+)",
+            "array_shape"_a,
+            "fill_value"_a,
+            pybind11::kw_only(),
+            "partition_shape"_a = std::optional<pybind11::tuple>{},
+            pybind11::return_value_policy::move);
     }
 
 }  // namespace lue::framework
diff --git a/source/framework/python/src/core/scalar.cpp b/source/framework/python/src/core/scalar.cpp
index d7dc19caa..590772dda 100644
--- a/source/framework/python/src/core/scalar.cpp
+++ b/source/framework/python/src/core/scalar.cpp
@@ -1,10 +1,120 @@
 #include "lue/py/framework/core/scalar.hpp"
 
 
+using namespace pybind11::literals;
+
+
 namespace lue::framework {
+    namespace {
+
+        template<typename Element>
+        auto create_scalar(pybind11::object const& value) -> Scalar<Element>
+        {
+            return pybind11::cast<Element>(value);
+        }
+
+
+        auto create_scalar_py(pybind11::object const& dtype_args, pybind11::object const& value)
+            -> pybind11::object
+        {
+            pybind11::dtype const dtype{pybind11::dtype::from_args(dtype_args)};
+
+            // Switch on dtype and return a scalar of the right type
+
+            auto const kind = dtype.kind();
+            auto const size = dtype.itemsize();  // bytes
+            pybind11::object result;
+
+            switch (kind)
+            {
+                case 'i':
+                {
+                    // Signed integer
+                    switch (size)
+                    {
+                        case 4:
+                        {
+                            result = pybind11::cast(create_scalar<std::int32_t>(value));
+                            break;
+                        }
+                        case 8:
+                        {
+                            result = pybind11::cast(create_scalar<std::int64_t>(value));
+                            break;
+                        }
+                    }
+
+                    break;
+                }
+                case 'u':
+                {
+                    // Unsigned integer
+                    switch (size)
+                    {
+                        case 1:
+                        {
+                            result = pybind11::cast(create_scalar<std::uint8_t>(value));
+                            break;
+                        }
+                        case 4:
+                        {
+                            result = pybind11::cast(create_scalar<std::uint32_t>(value));
+                            break;
+                        }
+                        case 8:
+                        {
+                            result = pybind11::cast(create_scalar<std::uint64_t>(value));
+                            break;
+                        }
+                    }
+
+                    break;
+                }
+                case 'f':
+                {
+                    // Floating-point
+                    switch (size)
+                    {
+                        case 4:
+                        {
+                            result = pybind11::cast(create_scalar<float>(value));
+                            break;
+                        }
+                        case 8:
+                        {
+                            result = pybind11::cast(create_scalar<double>(value));
+                            break;
+                        }
+                    }
+
+                    break;
+                }
+            }
+
+            if (!result)
+            {
+                throw std::runtime_error(fmt::format("Unsupported dtype (kind={}, itemsize={})", kind, size));
+            }
+
+            return result;
+        }
+
+    }  // Anonymous namespace
 
     void bind_scalar(pybind11::module& module)
     {
+        module.def(
+            "create_scalar",
+            create_scalar_py,
+            R"(
+    Create new scalar, with a value
+
+    :param numpy.dtype dtype: Description of the type of the scalar value
+    :param value: Value to initialize scalar with
+)",
+            "dtype"_a,
+            "value"_a);
+
         bind_scalar<std::uint8_t>(module);
         bind_scalar<std::uint32_t>(module);
         bind_scalar<std::uint64_t>(module);
diff --git a/source/framework/python/src/gdal/from_gdal.cpp b/source/framework/python/src/gdal/from_gdal.cpp
index 398c75e15..5d3e635df 100644
--- a/source/framework/python/src/gdal/from_gdal.cpp
+++ b/source/framework/python/src/gdal/from_gdal.cpp
@@ -13,14 +13,14 @@ namespace lue::framework {
     namespace {
 
         template<typename Element>
-        pybind11::object from_gdal(std::string const& name, StaticShape<2> const& partition_shape)
+        auto from_gdal(std::string const& name, StaticShape<2> const& partition_shape) -> pybind11::object
         {
             return pybind11::cast(read<Element>(name, partition_shape));
         }
 
 
-        pybind11::object from_gdal_py(
-            std::string const& name, std::optional<pybind11::tuple> const& partition_shape)
+        auto from_gdal_py(std::string const& name, std::optional<pybind11::tuple> const& partition_shape)
+            -> pybind11::object
         {
             gdal::DatasetPtr dataset{gdal::open_dataset(name, GDALAccess::GA_ReadOnly)};
             auto const raster_shape{gdal::shape(*dataset)};
@@ -95,6 +95,21 @@ namespace lue::framework {
             return result;
         }
 
+
+        auto probe_gdal(std::string const& name) -> pybind11::dict
+        {
+            // Open GDAL dataset and return its properties
+            gdal::DatasetPtr dataset{gdal::open_dataset(name, GDALAccess::GA_ReadOnly)};
+            auto const raster_shape{gdal::shape(*dataset)};
+            auto const geo_transform{gdal::geo_transform(*dataset)};
+
+            pybind11::dict properties{};
+            properties("shape"_a = raster_shape);
+            properties("geo_transform"_a = geo_transform);
+
+            return properties;
+        }
+
     }  // Anonymous namespace
 
 
@@ -116,6 +131,17 @@ namespace lue::framework {
             pybind11::kw_only(),
             "partition_shape"_a = std::optional<pybind11::tuple>{},
             pybind11::return_value_policy::move);
+
+        module.def(
+            "probe_gdal",
+            probe_gdal,
+            R"(
+    Open raster using GDAL and return its properties
+
+    :param str name: Name of dataset to probe
+    :rtype: Dictionary with information found
+)",
+            "name"_a);
     }
 
 }  // namespace lue::framework


Test patch:
diff --git a/source/framework/algorithm/test/create_partitioned_array_test.cpp b/source/framework/algorithm/test/create_partitioned_array_test.cpp
index 532ba4420..ef6312b68 100644
--- a/source/framework/algorithm/test/create_partitioned_array_test.cpp
+++ b/source/framework/algorithm/test/create_partitioned_array_test.cpp
@@ -55,9 +55,8 @@ namespace {
 
                 return hpx::async(
 
-                    [locality_id, offset, partition_shape, fill_value]() {
-                        return Partition{locality_id, offset, partition_shape, fill_value};
-                    }
+                    [locality_id, offset, partition_shape, fill_value]()
+                    { return Partition{locality_id, offset, partition_shape, fill_value}; }
 
                 );
             }
@@ -87,9 +86,8 @@ namespace {
 
         partitions[0] = hpx::async(
 
-            [locality_id, offset = offsets[0], partition_shape = partition_shapes[0], fill_value]() {
-                return Partition{locality_id, offset, partition_shape, fill_value};
-            }
+            [locality_id, offset = offsets[0], partition_shape = partition_shapes[0], fill_value]()
+            { return Partition{locality_id, offset, partition_shape, fill_value}; }
 
         );
 
@@ -98,9 +96,8 @@ namespace {
             partitions[idx] = partitions[idx - 1].then(
 
                 [locality_id, offset = offsets[idx], partition_shape = partition_shapes[idx], fill_value](
-                    auto const& /* previous_partition */) {
-                    return Partition{locality_id, offset, partition_shape, fill_value};
-                }
+                    auto const& /* previous_partition */)
+                { return Partition{locality_id, offset, partition_shape, fill_value}; }
 
             );
         }
@@ -167,34 +164,6 @@ namespace {
 }  // Anonymous namespace
 
 
-namespace lue {
-
-    template<typename Element, Rank rank>
-    class FunctorTraits<NumberPartitionsIndividually<Element, rank>>
-    {
-
-        public:
-
-            static constexpr bool const is_functor{true};
-    };
-
-}  // namespace lue
-
-
-namespace lue {
-
-    template<typename Element, Rank rank>
-    class FunctorTraits<NumberPartitionsPerLocality<Element, rank>>
-    {
-
-        public:
-
-            static constexpr bool const is_functor{true};
-    };
-
-}  // namespace lue
-
-
 BOOST_AUTO_TEST_CASE(instantiate_partitions_individually)
 {
     // Create a partitioned array in which the elements of the partitions
diff --git a/source/framework/algorithm/test/focal_operation_test.cpp b/source/framework/algorithm/test/focal_operation_test.cpp
index abc041b50..bd23ecc79 100644
--- a/source/framework/algorithm/test/focal_operation_test.cpp
+++ b/source/framework/algorithm/test/focal_operation_test.cpp
@@ -15,6 +15,8 @@ class Functor
 
         using OutputElement = OutputElement_;
 
+        static constexpr char const* name{"functor"};
+
         Functor() = default;
 
         template<typename Kernel, typename OutputPolicies, typename InputPolicies, typename Subspan1>
diff --git a/source/framework/python/test/core/scalar_test.py b/source/framework/python/test/core/scalar_test.py
index 80f954e2e..bb28e1a43 100644
--- a/source/framework/python/test/core/scalar_test.py
+++ b/source/framework/python/test/core/scalar_test.py
@@ -99,3 +99,81 @@ def test_operator_subtract_overloads(self):
             _ = value - scalar
             scalar -= value
             value -= scalar
+
+    # TODO @lue_test.framework_test_case
+    # TODO def test_operator_gt_overloads(self):
+    # TODO     for type_ in self.numeric_types:
+    # TODO         array = self.array[type_]
+    # TODO         scalar = self.scalar[type_]
+    # TODO         value = self.value[type_]
+
+    # TODO         _ = array > array
+    # TODO         _ = array > scalar
+    # TODO         _ = array > value
+    # TODO         _ = scalar > array
+    # TODO         # _ = value > array
+
+    # TODO @lue_test.framework_test_case
+    # TODO def test_operator_ge_overloads(self):
+    # TODO     for type_ in self.numeric_types:
+    # TODO         array = self.array[type_]
+    # TODO         scalar = self.scalar[type_]
+    # TODO         value = self.value[type_]
+
+    # TODO         _ = array >= array
+    # TODO         _ = array >= scalar
+    # TODO         _ = array >= value
+    # TODO         _ = scalar >= array
+    # TODO         # _ = value >= array
+
+    # TODO @lue_test.framework_test_case
+    # TODO def test_operator_eq_overloads(self):
+    # TODO     for type_ in self.numeric_types:
+    # TODO         array = self.array[type_]
+    # TODO         scalar = self.scalar[type_]
+    # TODO         value = self.value[type_]
+
+    # TODO         _ = array == array
+    # TODO         _ = array == scalar
+    # TODO         _ = array == value
+    # TODO         _ = scalar == array
+    # TODO         # _ = value == array
+
+    # TODO @lue_test.framework_test_case
+    # TODO def test_operator_ne_overloads(self):
+    # TODO     for type_ in self.numeric_types:
+    # TODO         array = self.array[type_]
+    # TODO         scalar = self.scalar[type_]
+    # TODO         value = self.value[type_]
+
+    # TODO         _ = array != array
+    # TODO         _ = array != scalar
+    # TODO         _ = array != value
+    # TODO         _ = scalar != array
+    # TODO         # _ = value != array
+
+    # TODO @lue_test.framework_test_case
+    # TODO def test_operator_lt_overloads(self):
+    # TODO     for type_ in self.numeric_types:
+    # TODO         array = self.array[type_]
+    # TODO         scalar = self.scalar[type_]
+    # TODO         value = self.value[type_]
+
+    # TODO         _ = array < array
+    # TODO         _ = array < scalar
+    # TODO         _ = array < value
+    # TODO         _ = scalar < array
+    # TODO         # _ = value < array
+
+    # TODO @lue_test.framework_test_case
+    # TODO def test_operator_le_overloads(self):
+    # TODO     for type_ in self.numeric_types:
+    # TODO         array = self.array[type_]
+    # TODO         scalar = self.scalar[type_]
+    # TODO         value = self.value[type_]
+
+    # TODO         _ = array <= array
+    # TODO         _ = array <= scalar
+    # TODO         _ = array <= value
+    # TODO         _ = scalar <= array
+    # TODO         # _ = value <= array
diff --git a/source/framework/python/test/pcraster/conversion_operation_test.py b/source/framework/python/test/pcraster/conversion_operation_test.py
index 4ae307e9f..0d8105e7d 100644
--- a/source/framework/python/test/pcraster/conversion_operation_test.py
+++ b/source/framework/python/test/pcraster/conversion_operation_test.py
@@ -35,14 +35,14 @@ def test_boolean(self):
     def test_nominal(self):
         for expression_type in [
             int,
-            float,
+            # float,
             np.int32,
             np.int64,
             np.uint8,
             np.uint32,
             np.uint64,
-            np.float32,
-            np.float64,
+            # np.float32,
+            # np.float64,
         ]:
             non_spatial = expression_type(5)
 
diff --git a/source/framework/python/test/pcraster/global_operation_test.py b/source/framework/python/test/pcraster/global_operation_test.py
index 4b30d67d7..b9b8b6a92 100644
--- a/source/framework/python/test/pcraster/global_operation_test.py
+++ b/source/framework/python/test/pcraster/global_operation_test.py
@@ -23,9 +23,10 @@ def test_mapminimum(self):
 
             _ = lpr.mapminimum(expression)
 
-    @lue_test.framework_test_case
-    def test_mapnormal(self):
-        _ = lpr.mapnormal()
+    # TODO Return a non-spatial
+    # @lue_test.framework_test_case
+    # def test_mapnormal(self):
+    #     _ = lpr.mapnormal()
 
     @lue_test.framework_test_case
     def test_maptotal(self):
@@ -34,6 +35,7 @@ def test_maptotal(self):
 
             _ = lpr.maptotal(expression)
 
-    @lue_test.framework_test_case
-    def test_mapuniform(self):
-        _ = lpr.mapuniform()
+    # TODO Return a non-spatial
+    # @lue_test.framework_test_case
+    # def test_mapuniform(self):
+    #     _ = lpr.mapuniform()
diff --git a/source/framework/python/test/pcraster/local_operation_test.py b/source/framework/python/test/pcraster/local_operation_test.py
index f1956b62e..7b4e1a424 100644
--- a/source/framework/python/test/pcraster/local_operation_test.py
+++ b/source/framework/python/test/pcraster/local_operation_test.py
@@ -6,13 +6,16 @@
 from .operation_test import OperationTest, setUpModule, tearDownModule
 
 
+# TODO Add tests for overload accepting value instead of non_spatial
+
+
 class LocalOperationTest(OperationTest):
     @lue_test.framework_test_case
     def test_add(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial + spatial
@@ -25,7 +28,7 @@ def test_max(self):
         for expression_type in [np.uint8, np.int32, np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.max(spatial, spatial)
@@ -38,7 +41,7 @@ def test_min(self):
         for expression_type in [np.uint8, np.int32, np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.min(spatial, spatial)
@@ -51,26 +54,26 @@ def test_subtract(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial - spatial
             _ = spatial - non_spatial
             _ = non_spatial - spatial
-            # _ = non_spatial - non_spatial
+            _ = non_spatial - non_spatial
 
     @lue_test.framework_test_case
     def test_divide(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial / spatial
             _ = spatial / non_spatial
             _ = non_spatial / spatial
-            # _ = non_spatial / non_spatial
+            _ = non_spatial / non_spatial
 
             _ = lpr.div(spatial, spatial)
             _ = lpr.div(spatial, non_spatial)
@@ -82,33 +85,33 @@ def test_multiply(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial * spatial
             _ = spatial * non_spatial
             _ = non_spatial * spatial
-            # _ = non_spatial * non_spatial
+            _ = non_spatial * non_spatial
 
     @lue_test.framework_test_case
     def test_power(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial**spatial
             _ = spatial**non_spatial
             _ = non_spatial**spatial
-            # _ = non_spatial**non_spatial
+            _ = non_spatial**non_spatial
 
     @lue_test.framework_test_case
     def test_abs(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.abs(spatial)
@@ -119,7 +122,7 @@ def test_acos(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.acos(spatial)
@@ -130,25 +133,25 @@ def test_and(self):
         for expression_type in [np.uint8]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.pcrand(spatial, spatial)
-            _ = lpr.pcrand(spatial, spatial)
-            _ = lpr.pcrand(non_spatial, non_spatial)
+            _ = lpr.pcrand(spatial, non_spatial)
+            _ = lpr.pcrand(non_spatial, spatial)
             _ = lpr.pcrand(non_spatial, non_spatial)
 
             _ = spatial & spatial
-            _ = spatial & spatial
+            _ = spatial & non_spatial
+            _ = non_spatial & spatial
             _ = non_spatial & non_spatial
-            # _ = non_spatial & non_spatial
 
     @lue_test.framework_test_case
     def test_asin(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.asin(spatial)
@@ -168,7 +171,7 @@ def test_cos(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.cos(spatial)
@@ -184,15 +187,15 @@ def test_cover(self):
 
             _ = lpr.cover(spatial, spatial)
             _ = lpr.cover(spatial, non_spatial)
-            _ = lpr.cover(non_spatial, non_spatial)
             _ = lpr.cover(non_spatial, spatial)
+            _ = lpr.cover(non_spatial, non_spatial)
 
     @lue_test.framework_test_case
     def test_defined(self):
         for expression_type in [np.uint8, np.int32, np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.defined(spatial)
@@ -203,13 +206,13 @@ def test_eq(self):
         for expression_type in [np.uint8, np.int32, np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial == spatial
             _ = spatial == non_spatial
             _ = non_spatial == spatial
-            # _ = non_spatial == non_spatial
+            _ = non_spatial == non_spatial
 
             _ = lpr.eq(spatial, spatial)
             _ = lpr.eq(spatial, non_spatial)
@@ -221,7 +224,7 @@ def test_exp(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.exp(spatial)
@@ -232,13 +235,13 @@ def test_ge(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial >= spatial
             _ = spatial >= non_spatial
             _ = non_spatial >= spatial
-            # _ = non_spatial >= non_spatial
+            _ = non_spatial >= non_spatial
 
             _ = lpr.ge(spatial, spatial)
             _ = lpr.ge(spatial, non_spatial)
@@ -250,13 +253,13 @@ def test_gt(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial > spatial
             _ = spatial > non_spatial
             _ = non_spatial > spatial
-            # _ = non_spatial > non_spatial
+            _ = non_spatial > non_spatial
 
             _ = lpr.gt(spatial, spatial)
             _ = lpr.gt(spatial, non_spatial)
@@ -267,13 +270,13 @@ def test_gt(self):
     def test_ifthen(self):
         spatial_condition, non_spatial_condition = (
             self.spatial[np.uint8],
-            self.non_spatial[np.uint8].item(),
+            self.non_spatial[np.uint8],
         )
 
         for expression_type in [np.uint8, np.int32, np.float32]:
             spatial_expression, non_spatial_expression = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.ifthen(spatial_condition, spatial_expression)
@@ -285,13 +288,13 @@ def test_ifthen(self):
     def test_ifthenelse(self):
         spatial_condition, non_spatial_condition = (
             self.spatial[np.uint8],
-            self.non_spatial[np.uint8].item(),
+            self.non_spatial[np.uint8],
         )
 
         for expression_type in [np.uint8, np.int32, np.float32]:
             spatial_expression, non_spatial_expression = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.ifthenelse(
@@ -324,13 +327,13 @@ def test_le(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial <= spatial
             _ = spatial <= non_spatial
             _ = non_spatial <= spatial
-            # _ = non_spatial <= non_spatial
+            _ = non_spatial <= non_spatial
 
             _ = lpr.le(spatial, spatial)
             _ = lpr.le(spatial, non_spatial)
@@ -353,13 +356,13 @@ def test_lt(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial < spatial
             _ = spatial < non_spatial
             _ = non_spatial < spatial
-            # _ = non_spatial < non_spatial
+            _ = non_spatial < non_spatial
 
             _ = lpr.lt(spatial, spatial)
             _ = lpr.lt(spatial, non_spatial)
@@ -371,13 +374,13 @@ def test_ne(self):
         for expression_type in [np.uint8, np.int32, np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = spatial != spatial
             _ = spatial != non_spatial
             _ = non_spatial != spatial
-            # _ = non_spatial != non_spatial
+            _ = non_spatial != non_spatial
 
             _ = lpr.ne(spatial, spatial)
             _ = lpr.ne(spatial, non_spatial)
@@ -389,7 +392,7 @@ def test_normal(self):
         for expression_type in [np.uint8]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.normal(spatial)
@@ -400,39 +403,39 @@ def test_not(self):
         for expression_type in [np.uint8]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.pcrnot(spatial)
             _ = lpr.pcrnot(non_spatial)
 
             _ = ~spatial
-            # _ = ~non_spatial
+            _ = ~non_spatial
 
     @lue_test.framework_test_case
     def test_or(self):
         for expression_type in [np.uint8]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.pcror(spatial, spatial)
-            _ = lpr.pcror(spatial, spatial)
-            _ = lpr.pcror(non_spatial, non_spatial)
+            _ = lpr.pcror(spatial, non_spatial)
+            _ = lpr.pcror(non_spatial, spatial)
             _ = lpr.pcror(non_spatial, non_spatial)
 
             _ = spatial | spatial
-            _ = spatial | spatial
-            # _ = non_spatial | non_spatial
-            # _ = non_spatial | non_spatial
+            _ = spatial | non_spatial
+            _ = non_spatial | spatial
+            _ = non_spatial | non_spatial
 
     @lue_test.framework_test_case
     def test_rounddown(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.rounddown(spatial)
@@ -443,7 +446,7 @@ def test_roundoff(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.roundoff(spatial)
@@ -454,7 +457,7 @@ def test_roundup(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.roundup(spatial)
@@ -465,7 +468,7 @@ def test_sin(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.sin(spatial)
@@ -487,7 +490,7 @@ def test_sqr(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.sqr(spatial)
@@ -498,7 +501,7 @@ def test_sqrt(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.sqrt(spatial)
@@ -509,7 +512,7 @@ def test_tan(self):
         for expression_type in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.tan(spatial)
@@ -520,7 +523,7 @@ def test_uniform(self):
         for expression_type in [np.uint8]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.uniform(spatial)
@@ -531,7 +534,7 @@ def test_uniqueid(self):
         for expression_type in [np.uint8]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.uniqueid(spatial)
@@ -542,15 +545,15 @@ def test_xor(self):
         for expression_type in [np.uint8]:
             spatial, non_spatial = (
                 self.spatial[expression_type],
-                self.non_spatial[expression_type].item(),
+                self.non_spatial[expression_type],
             )
 
             _ = lpr.pcrxor(spatial, spatial)
-            _ = lpr.pcrxor(spatial, spatial)
-            _ = lpr.pcrxor(non_spatial, non_spatial)
+            _ = lpr.pcrxor(spatial, non_spatial)
+            _ = lpr.pcrxor(non_spatial, spatial)
             _ = lpr.pcrxor(non_spatial, non_spatial)
 
             _ = spatial | spatial
-            _ = spatial | spatial
-            # _ = non_spatial | non_spatial
-            # _ = non_spatial | non_spatial
+            _ = spatial | non_spatial
+            _ = non_spatial | spatial
+            _ = non_spatial | non_spatial
diff --git a/source/framework/python/test/pcraster/operation_test.py b/source/framework/python/test/pcraster/operation_test.py
index 2d4b248ee..b0cab9fd1 100644
--- a/source/framework/python/test/pcraster/operation_test.py
+++ b/source/framework/python/test/pcraster/operation_test.py
@@ -27,60 +27,81 @@ def setUpClass(cls):
         lpr.configuration.partition_shape = cls.partition_shape
         lpr.configuration.cell_size = 10
 
+        cls.numeric_types = lue_test.numeric_types
+        cls.dtype_by_type = lue_test.dtype_by_type
+        cls.value = lue_test.value_by_type
         cls.non_spatial = {
-            np.uint8: np.uint8(1),
-            np.uint32: np.uint32(32),
-            np.uint64: np.uint32(64),
-            np.int32: np.int32(-32),
-            np.int64: np.int32(-64),
-            np.float32: np.float32(3.2),
-            np.float64: np.float64(6.4),
+            type_: lfr.create_scalar(cls.dtype_by_type[type_], cls.value[type_])
+            for type_ in cls.numeric_types
         }
         cls.spatial = {
-            np.uint8: lfr.create_array(
+            type_: lfr.create_array(
                 cls.array_shape,
-                np.uint8,
-                fill_value=cls.non_spatial[np.uint8],
+                type_,
+                fill_value=value,
                 partition_shape=cls.partition_shape,
-            ),
-            np.uint32: lfr.create_array(
-                cls.array_shape,
-                np.uint32,
-                fill_value=cls.non_spatial[np.uint32],
-                partition_shape=cls.partition_shape,
-            ),
-            np.uint64: lfr.create_array(
-                cls.array_shape,
-                np.uint64,
-                fill_value=cls.non_spatial[np.uint64],
-                partition_shape=cls.partition_shape,
-            ),
-            np.int32: lfr.create_array(
-                cls.array_shape,
-                np.int32,
-                fill_value=cls.non_spatial[np.int32],
-                partition_shape=cls.partition_shape,
-            ),
-            np.int64: lfr.create_array(
-                cls.array_shape,
-                np.int64,
-                fill_value=cls.non_spatial[np.int64],
-                partition_shape=cls.partition_shape,
-            ),
-            np.float32: lfr.create_array(
-                cls.array_shape,
-                np.float32,
-                fill_value=cls.non_spatial[np.float32],
-                partition_shape=cls.partition_shape,
-            ),
-            np.float64: lfr.create_array(
-                cls.array_shape,
-                np.float64,
-                fill_value=cls.non_spatial[np.float64],
-                partition_shape=cls.partition_shape,
-            ),
+            )
+            for type_, value in cls.value.items()
         }
+
+
+
+        # cls.non_spatial = {
+        #     np.uint8: np.uint8(1),
+        #     np.uint32: np.uint32(32),
+        #     np.uint64: np.uint64(64),
+        #     np.int32: np.int32(-32),
+        #     np.int64: np.int64(-64),
+        #     np.float32: np.float32(3.2),
+        #     np.float64: np.float64(6.4),
+        # }
+        # cls.spatial = {
+        #     np.uint8: lfr.create_array(
+        #         cls.array_shape,
+        #         np.uint8,
+        #         fill_value=cls.non_spatial[np.uint8],
+        #         partition_shape=cls.partition_shape,
+        #     ),
+        #     np.uint32: lfr.create_array(
+        #         cls.array_shape,
+        #         np.uint32,
+        #         fill_value=cls.non_spatial[np.uint32],
+        #         partition_shape=cls.partition_shape,
+        #     ),
+        #     np.uint64: lfr.create_array(
+        #         cls.array_shape,
+        #         np.uint64,
+        #         fill_value=cls.non_spatial[np.uint64],
+        #         partition_shape=cls.partition_shape,
+        #     ),
+        #     np.int32: lfr.create_array(
+        #         cls.array_shape,
+        #         np.int32,
+        #         fill_value=cls.non_spatial[np.int32],
+        #         partition_shape=cls.partition_shape,
+        #     ),
+        #     np.int64: lfr.create_array(
+        #         cls.array_shape,
+        #         np.int64,
+        #         fill_value=cls.non_spatial[np.int64],
+        #         partition_shape=cls.partition_shape,
+        #     ),
+        #     np.float32: lfr.create_array(
+        #         cls.array_shape,
+        #         np.float32,
+        #         fill_value=cls.non_spatial[np.float32],
+        #         partition_shape=cls.partition_shape,
+        #     ),
+        #     np.float64: lfr.create_array(
+        #         cls.array_shape,
+        #         np.float64,
+        #         fill_value=cls.non_spatial[np.float64],
+        #         partition_shape=cls.partition_shape,
+        #     ),
+        # }
+
         direction = 3
+
         cls.ldd = lfr.create_array(
             cls.array_shape, np.uint8, direction, partition_shape=cls.partition_shape
         )
diff --git a/source/framework/python/test/pcraster/routing_operation_test.py b/source/framework/python/test/pcraster/routing_operation_test.py
index fd56fcda6..2add4d547 100644
--- a/source/framework/python/test/pcraster/routing_operation_test.py
+++ b/source/framework/python/test/pcraster/routing_operation_test.py
@@ -15,7 +15,7 @@ def test_accufraction(self):
         for type_ in [np.float32]:
             spatial_material, non_spatial_material = (
                 self.spatial[type_],
-                self.non_spatial[type_].item(),
+                self.non_spatial[type_],
             )
 
             # Fraction must be with [0, 1]
@@ -41,7 +41,7 @@ def test_accuthreshold(self):
         for type_ in [np.float32]:
             spatial_material, non_spatial_material = (
                 self.spatial[type_],
-                self.non_spatial[type_].item(),
+                self.non_spatial[type_],
             )
 
             non_spatial_threshold = type_(0.5) * non_spatial_material
@@ -79,7 +79,7 @@ def test_kinematic(self):
         for type_ in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[type_],
-                self.non_spatial[type_].item(),
+                self.value[type_],  # TODO non_spatial (Scalar) is not supported yet
             )
 
             _ = lpr.kinematic(
@@ -98,7 +98,7 @@ def test_lddcreate(self):
         for type_ in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[type_],
-                self.non_spatial[type_].item(),
+                self.non_spatial[type_],
             )
 
             _ = lpr.lddcreate(
@@ -116,7 +116,7 @@ def test_upstream(self):
         for type_ in [np.float32]:
             spatial, non_spatial = (
                 self.spatial[type_],
-                self.non_spatial[type_].item(),
+                self.value[type_],  # TODO non_spatial (Scalar) is not supported yet
             )
 
             _ = lpr.upstream(ldd, non_spatial)