From ff1ff37627d1ae0776e45a1c6588c81e4e108e6c Mon Sep 17 00:00:00 2001 From: Kor de Jong Date: Thu, 23 Jan 2025 17:07:10 +0100 Subject: [PATCH] Add support for reading hyperslab using gdal --- .../framework/core/domain_decomposition.hpp | 2 +- .../include/lue/framework/core/hyperslab.hpp | 113 +++++++++++++ source/framework/core/test/CMakeLists.txt | 1 + source/framework/core/test/hyperslab_test.cpp | 132 +++++++++++++++ .../io/include/lue/framework/io/gdal.hpp | 7 + source/framework/io/source/gdal.cpp.in | 116 +++++++++++-- source/framework/io/test/gdal_test.cpp | 156 +++++++++++++++++- source/framework/python/CMakeLists.txt | 1 + .../python/source/core/hyperslab.cpp | 48 ++++++ .../python/source/gdal/from_gdal.cpp | 143 +++++++++++++++- source/framework/python/source/submodule.cpp | 2 + source/framework/python/test/gdal_test.py | 50 +++++- 12 files changed, 739 insertions(+), 32 deletions(-) create mode 100644 source/framework/core/include/lue/framework/core/hyperslab.hpp create mode 100644 source/framework/core/test/hyperslab_test.cpp create mode 100644 source/framework/python/source/core/hyperslab.cpp diff --git a/source/framework/core/include/lue/framework/core/domain_decomposition.hpp b/source/framework/core/include/lue/framework/core/domain_decomposition.hpp index 084d970b2..9cba28229 100644 --- a/source/framework/core/include/lue/framework/core/domain_decomposition.hpp +++ b/source/framework/core/include/lue/framework/core/domain_decomposition.hpp @@ -361,7 +361,7 @@ namespace lue { template - Shape default_partition_shape(Shape const& array_shape) + auto default_partition_shape(Shape const& array_shape) -> Shape { Count const nr_worker_threads{static_cast(hpx::get_num_worker_threads())}; diff --git a/source/framework/core/include/lue/framework/core/hyperslab.hpp b/source/framework/core/include/lue/framework/core/hyperslab.hpp new file mode 100644 index 000000000..f355a6dd8 --- /dev/null +++ b/source/framework/core/include/lue/framework/core/hyperslab.hpp @@ -0,0 +1,113 @@ +#pragma once +#include "lue/framework/core/offset.hpp" +#include +#include +#include + + +namespace lue { + + template + class Hyperslab + { + + private: + + using Array = std::array; + + public: + + using Offsets = Array; + using Counts = Array; + using Strides = Array; + using Center = Array; + using Shape = Array; + + + Hyperslab(Center const& center, Shape const& shape): + + _counts{shape} + + { + Shape half_shape{shape}; + + std::transform( + half_shape.begin(), + half_shape.end(), + half_shape.begin(), + [](auto const count) { return 0.5 * (count % 2 == 1 ? count - 1 : count); }); + + // Equality in the sense that shape can be centered and not extend into areas with smaller + // than zero idxs + if (!std::equal( + half_shape.begin(), + half_shape.end(), + center.begin(), + [](auto const count, auto const idx) + { + // Subtracting this many cells from the center cell results in a negative idx... + return idx >= count; + })) + { + throw std::runtime_error("Hyperslab shape extents beyond the limits of the array"); + } + + std::transform( + center.begin(), + center.end(), + half_shape.begin(), + _offsets.begin(), + [](auto const idx, auto const count) + { + assert(idx >= count); + return idx - count; + }); + _strides.fill(1); + } + + auto offsets() const -> Offsets const& + { + return _offsets; + } + + auto counts() const -> Counts const& + { + return _counts; + } + + auto strides() const -> Strides const& + { + return _strides; + } + + [[nodiscard]] auto is_strided() const -> bool + { + return std::ranges::any_of(_strides, [](auto const stride) { return stride > 1; }); + } + + // Returns whether the indices represented by the hyperslab all fall within the bounds of the + // array whose shape is passed in. By definition (see constructor), the minimum indices of the + // hyperslab are fine, so only the maximum indices are checked here. + auto is_contained_by(Shape const& array_shape) const -> bool + { + for (std::size_t idx = 0; idx < rank; ++idx) + { + if (!(_offsets[idx] + _counts[idx] - 1 < array_shape[idx])) + { + return false; + } + } + + return true; + } + + private: + + Offsets _offsets{}; + + Counts _counts; + + Strides _strides{}; + }; + +} // namespace lue diff --git a/source/framework/core/test/CMakeLists.txt b/source/framework/core/test/CMakeLists.txt index 3c6d75d62..d60e9c2a1 100644 --- a/source/framework/core/test/CMakeLists.txt +++ b/source/framework/core/test/CMakeLists.txt @@ -4,6 +4,7 @@ set(names configuration_entry domain_decomposition hilbert_curve + hyperslab index_util linear_curve math diff --git a/source/framework/core/test/hyperslab_test.cpp b/source/framework/core/test/hyperslab_test.cpp new file mode 100644 index 000000000..ee0535a3a --- /dev/null +++ b/source/framework/core/test/hyperslab_test.cpp @@ -0,0 +1,132 @@ +#define BOOST_TEST_MODULE lue framework core hyperslab +#include "lue/framework/core/hyperslab.hpp" +#include "lue/stream.hpp" +#include + + +BOOST_AUTO_TEST_CASE(construct_2d_odd_shape) +{ + lue::Rank const rank{2}; + using Hyperslab = lue::Hyperslab; + + Hyperslab::Center center{20, 30}; + Hyperslab::Shape shape{7, 5}; + Hyperslab hyperslab{center, shape}; + + Hyperslab::Offsets const offsets{17, 28}; + Hyperslab::Counts const counts{shape}; + Hyperslab::Strides const strides{{1, 1}}; + + BOOST_CHECK_EQUAL(hyperslab.offsets(), offsets); + BOOST_CHECK_EQUAL(hyperslab.counts(), counts); + BOOST_CHECK_EQUAL(hyperslab.strides(), strides); +} + + +BOOST_AUTO_TEST_CASE(construct_2d_even_shape) +{ + lue::Rank const rank{2}; + using Hyperslab = lue::Hyperslab; + + Hyperslab::Center center{20, 30}; + Hyperslab::Shape shape{6, 4}; + Hyperslab hyperslab{center, shape}; + + Hyperslab::Offsets const offsets{17, 28}; + Hyperslab::Counts const counts{shape}; + Hyperslab::Strides const strides{{1, 1}}; + + BOOST_CHECK_EQUAL(hyperslab.offsets(), offsets); + BOOST_CHECK_EQUAL(hyperslab.counts(), counts); + BOOST_CHECK_EQUAL(hyperslab.strides(), strides); +} + + +BOOST_AUTO_TEST_CASE(valid_hyperslab) +{ + lue::Rank const rank{2}; + using Hyperslab = lue::Hyperslab; + + { + // Empty hyperslab + Hyperslab::Center const center{30, 20}; + Hyperslab::Shape const shape{0, 0}; + BOOST_CHECK_NO_THROW(Hyperslab(center, shape)); + } + + { + // Empty hyperslab + Hyperslab::Center const center{30, 20}; + Hyperslab::Shape const shape{60, 0}; + BOOST_CHECK_NO_THROW(Hyperslab(center, shape)); + } + + { + // Empty hyperslab + Hyperslab::Center const center{30, 20}; + Hyperslab::Shape const shape{0, 40}; + BOOST_CHECK_NO_THROW(Hyperslab(center, shape)); + } + + { + // Hyperslab corresponds with whole array + Hyperslab::Center const center{30, 20}; + Hyperslab::Shape const shape{60, 40}; + BOOST_CHECK_NO_THROW(Hyperslab(center, shape)); + } +} + + +BOOST_AUTO_TEST_CASE(invalid_hyperslab) +{ + lue::Rank const rank{2}; + using Hyperslab = lue::Hyperslab; + + { + Hyperslab::Center const center{30, 20}; + Hyperslab::Shape const shape{62, 42}; + + BOOST_CHECK_EXCEPTION( + Hyperslab(center, shape), + std::runtime_error, + [](auto const& exception) + { return std::string(exception.what()).find("extents beyond the limits") != std::string::npos; }); + } + + { + Hyperslab::Center const center{30, 20}; + Hyperslab::Shape const shape{62, 40}; + + BOOST_CHECK_EXCEPTION( + Hyperslab(center, shape), + std::runtime_error, + [](auto const& exception) + { return std::string(exception.what()).find("extents beyond the limits") != std::string::npos; }); + } + + { + Hyperslab::Center const center{30, 20}; + Hyperslab::Shape const shape{60, 42}; + + BOOST_CHECK_EXCEPTION( + Hyperslab(center, shape), + std::runtime_error, + [](auto const& exception) + { return std::string(exception.what()).find("extents beyond the limits") != std::string::npos; }); + } +} + + +BOOST_AUTO_TEST_CASE(is_contained_by) +{ + lue::Rank const rank{2}; + using Hyperslab = lue::Hyperslab; + + Hyperslab::Center center{20, 30}; + Hyperslab::Shape shape{7, 5}; + Hyperslab hyperslab{center, shape}; + + BOOST_CHECK(hyperslab.is_contained_by(Hyperslab::Shape{24, 33})); + BOOST_CHECK(!hyperslab.is_contained_by(Hyperslab::Shape{23, 33})); + BOOST_CHECK(!hyperslab.is_contained_by(Hyperslab::Shape{24, 32})); +} diff --git a/source/framework/io/include/lue/framework/io/gdal.hpp b/source/framework/io/include/lue/framework/io/gdal.hpp index 0df603216..034b9e7d2 100644 --- a/source/framework/io/include/lue/framework/io/gdal.hpp +++ b/source/framework/io/include/lue/framework/io/gdal.hpp @@ -1,4 +1,5 @@ #pragma once +#include "lue/framework/core/hyperslab.hpp" #include "lue/framework/partitioned_array.hpp" @@ -8,6 +9,12 @@ namespace lue { auto from_gdal(std::string const& name, Shape const& partition_shape) -> PartitionedArray; + template + auto from_gdal( + std::string const& name, + Hyperslab<2> const& hyperslab, + Shape const& partition_shape) -> PartitionedArray; + template auto to_gdal( PartitionedArray const& array, diff --git a/source/framework/io/source/gdal.cpp.in b/source/framework/io/source/gdal.cpp.in index 5d7a053f2..99ab590e3 100644 --- a/source/framework/io/source/gdal.cpp.in +++ b/source/framework/io/source/gdal.cpp.in @@ -9,6 +9,7 @@ #include #include #include +#include namespace lue { @@ -78,6 +79,22 @@ namespace lue { namespace { + void verify_compatibility( + Shape const& array_shape, + Hyperslab<2> const& hyperslab) + { + if(hyperslab.is_strided()) + { + throw std::runtime_error("Handling strided hyperslab is not supported yet"); + } + + if(!hyperslab.is_contained_by(array_shape)) + { + throw std::runtime_error("Hyperslab extents beyond array"); + } + } + + class Band { @@ -310,9 +327,10 @@ namespace lue { static constexpr bool instantiate_per_locality{true}; - ReadPartitionsPerLocality(std::string const& name): + ReadPartitionsPerLocality(std::string name, Hyperslab::Offsets const& offsets={0, 0}): - _name{name} + _name{std::move(name)}, + _offsets{offsets} { } @@ -329,6 +347,15 @@ namespace lue { using Action = ReadPartitionsPerLocalityAction; Action action{}; + for(auto& offset: offsets) + { + std::transform(_offsets.begin(), _offsets.end(), offset.begin(), offset.begin(), + [](auto offset1, auto const offset2) + { + return offset1 + offset2; + }); + } + return hpx::async( [locality_id, @@ -345,6 +372,9 @@ namespace lue { private: std::string _name; + + // Offset of all partitions within the raster + Hyperslab::Offsets _offsets; }; @@ -389,6 +419,25 @@ namespace lue { { }; + + template + auto meta(std::string const& name) -> std::tuple, Element, bool> + { + gdal::Raster raster{gdal::open_dataset(name, ::GA_ReadOnly)}; + gdal::Raster::Band raster_band{raster.band(1)}; + + Shape array_shape{}; + + auto const raster_shape = raster.shape(); + array_shape[0] = static_cast(raster_shape[0]); + array_shape[1] = static_cast(raster_shape[1]); + + auto [no_data_value, no_data_value_is_valid] = raster_band.no_data_value(); + + return {array_shape, no_data_value, no_data_value_is_valid}; + } + + } // Anonymous namespace @@ -416,19 +465,55 @@ namespace lue { // - Determine the no-data value to use // Only then can we call the logic that creates a partitioned array asynchronously. - Shape array_shape; - bool no_data_value_is_valid{}; - Element no_data_value; + auto const [array_shape, no_data_value, no_data_value_is_valid] = meta(name); + + using Functor = ReadPartitionsPerLocality; + + if (!no_data_value_is_valid) + { + // No no-data value is set in the band. Don't bother with no-data. + + using Policies = policy::from_gdal::DefaultPolicies; + Policies policies{}; + return create_partitioned_array(policies, array_shape, partition_shape, Functor{name}); + } + // A no-data value is set in the band. Let's be nice and use it. Upon reading, + // it will be converted to our own conventions for handling no-data. + + if (std::is_floating_point_v && std::isnan(no_data_value)) { - gdal::Raster raster{gdal::open_dataset(name, ::GA_ReadOnly)}; - gdal::Raster::Band raster_band{raster.band(1)}; + // No-data value is a NaN. This is what LUE uses as no-data for floating point + // by default. + using Policies = policy::from_gdal::ValuePolicies; + Policies policies{}; + return create_partitioned_array(policies, array_shape, partition_shape, Functor{name}); + } - auto const raster_shape = raster.shape(); - array_shape[0] = static_cast(raster_shape[0]); - array_shape[1] = static_cast(raster_shape[1]); - std::tie(no_data_value, no_data_value_is_valid) = raster_band.no_data_value(); + { + // No-data value is not NaN. Use a special no-data policy for this special value. + using Policies = policy::from_gdal::SpecialValuePolicy; + Policies policies{no_data_value}; + return create_partitioned_array(policies, array_shape, partition_shape, Functor{name}); } + } + } + + + template + auto from_gdal(std::string const& name, Hyperslab<2> const& hyperslab, Shape const& partition_shape) + -> PartitionedArray + { + if constexpr (!gdal::supports()) + { + gdal::verify_support(); + return PartitionedArray{}; + } + else + { + auto const [array_shape, no_data_value, no_data_value_is_valid] = meta(name); + + verify_compatibility(array_shape, hyperslab); using Functor = ReadPartitionsPerLocality; @@ -438,7 +523,7 @@ namespace lue { using Policies = policy::from_gdal::DefaultPolicies; Policies policies{}; - return create_partitioned_array(policies, array_shape, partition_shape, Functor{name}); + return create_partitioned_array(policies, hyperslab.counts(), partition_shape, Functor{name, hyperslab.offsets()}); } // A no-data value is set in the band. Let's be nice and use it. Upon reading, @@ -450,14 +535,14 @@ namespace lue { // by default. using Policies = policy::from_gdal::ValuePolicies; Policies policies{}; - return create_partitioned_array(policies, array_shape, partition_shape, Functor{name}); + return create_partitioned_array(policies, hyperslab.counts(), partition_shape, Functor{name, hyperslab.offsets()}); } { // No-data value is not NaN. Use a special no-data policy for this special value. using Policies = policy::from_gdal::SpecialValuePolicy; Policies policies{no_data_value}; - return create_partitioned_array(policies, array_shape, partition_shape, Functor{name}); + return create_partitioned_array(policies, hyperslab.counts(), partition_shape, Functor{name, hyperslab.offsets()}); } } } @@ -604,6 +689,9 @@ namespace lue { template lue::PartitionedArray<{{ Element }}, {{ rank }}> lue::from_gdal<{{ Element }}>( std::string const&, lue::Shape const&); + template lue::PartitionedArray<{{ Element }}, {{ rank }}> lue::from_gdal<{{ Element }}>( + std::string const&, lue::Hyperslab<{{ rank }}> const&, lue::Shape const&); + template hpx::future lue::to_gdal<{{ Element }}>( lue::PartitionedArray<{{ Element }}, {{ rank }}> const&, std::string const&, std::string const&); diff --git a/source/framework/io/test/gdal_test.cpp b/source/framework/io/test/gdal_test.cpp index fe714d1e6..1a8a34949 100644 --- a/source/framework/io/test/gdal_test.cpp +++ b/source/framework/io/test/gdal_test.cpp @@ -1,5 +1,8 @@ #define BOOST_TEST_MODULE lue framework io gdal #include "lue/framework/algorithm/create_partitioned_array.hpp" +#include "lue/framework/algorithm/range.hpp" +#include "lue/framework/algorithm/value_policies/add.hpp" +#include "lue/framework/algorithm/value_policies/cell_index.hpp" #include "lue/framework/io/gdal.hpp" #include "lue/framework/test/hpx_unit_test.hpp" #include "lue/framework.hpp" @@ -7,7 +10,7 @@ #include -BOOST_AUTO_TEST_CASE(round_trip_1) +BOOST_AUTO_TEST_CASE(array_all_valid) { // Signed int, all valid using Element = lue::SignedIntegralElement<0>; @@ -17,7 +20,7 @@ BOOST_AUTO_TEST_CASE(round_trip_1) Element const fill_value{5}; Array array_written{lue::create_partitioned_array(array_shape, partition_shape, fill_value)}; - std::string const name{"lue_framework_io_raster_round_trip_1.tif"}; + std::string const name{"lue_framework_io_gdal_array_int_all_valid.tif"}; lue::to_gdal(array_written, name).wait(); @@ -29,7 +32,7 @@ BOOST_AUTO_TEST_CASE(round_trip_1) } -BOOST_AUTO_TEST_CASE(round_trip_2) +BOOST_AUTO_TEST_CASE(array_none_valid) { // Signed int, none valid using Element = lue::SignedIntegralElement<0>; @@ -39,7 +42,7 @@ BOOST_AUTO_TEST_CASE(round_trip_2) Element const fill_value{lue::policy::no_data_value}; Array array_written{lue::create_partitioned_array(array_shape, partition_shape, fill_value)}; - std::string const name{"lue_framework_io_raster_round_trip_2.tif"}; + std::string const name{"lue_framework_io_gdal_array_int_none_valid.tif"}; lue::to_gdal(array_written, name).wait(); @@ -51,7 +54,7 @@ BOOST_AUTO_TEST_CASE(round_trip_2) } -BOOST_AUTO_TEST_CASE(round_trip_3) +BOOST_AUTO_TEST_CASE(array_float_all_valid) { // Float, all valid using Element = lue::FloatingPointElement<0>; @@ -61,7 +64,7 @@ BOOST_AUTO_TEST_CASE(round_trip_3) Element const fill_value{5}; Array array_written{lue::create_partitioned_array(array_shape, partition_shape, fill_value)}; - std::string const name{"lue_framework_io_raster_round_trip_3.tif"}; + std::string const name{"lue_framework_io_gdal_array_float_all_valid.tif"}; lue::to_gdal(array_written, name).wait(); @@ -73,7 +76,7 @@ BOOST_AUTO_TEST_CASE(round_trip_3) } -BOOST_AUTO_TEST_CASE(round_trip_4) +BOOST_AUTO_TEST_CASE(array_float_none_valid) { // Float, none valid using Element = lue::FloatingPointElement<0>; @@ -83,7 +86,7 @@ BOOST_AUTO_TEST_CASE(round_trip_4) Element const fill_value{lue::policy::no_data_value}; Array array_written{lue::create_partitioned_array(array_shape, partition_shape, fill_value)}; - std::string const name{"lue_framework_io_raster_round_trip_4.tif"}; + std::string const name{"lue_framework_io_gdal_array_float_none_valid.tif"}; lue::to_gdal(array_written, name).wait(); @@ -93,3 +96,140 @@ BOOST_AUTO_TEST_CASE(round_trip_4) lue::test::check_arrays_are_equal(array_read, array_written); } + + +template +using Array = lue::PartitionedArray; + + +BOOST_AUTO_TEST_CASE(hyperslab) +{ + using namespace lue::value_policies; + + using ConditionArray = Array; + using IndexArray = Array; + using Shape = lue::ShapeT; + using Hyperslab = lue::Hyperslab<2>; + + Shape const array_shape{60, 40}; + Shape const partition_shape{10, 10}; + + ConditionArray const condition{ + lue::create_partitioned_array(array_shape, partition_shape, 1)}; + IndexArray const row_idxs_written = cell_index(condition, 0); + IndexArray const col_idxs_written = cell_index(condition, 1); + + std::string const row_idxs_name{"lue_framework_io_gdal_hyperslab_row_idxs.tif"}; + std::string const col_idxs_name{"lue_framework_io_gdal_hyperslab_col_idxs.tif"}; + + lue::to_gdal(row_idxs_written, row_idxs_name).wait(); + lue::to_gdal(col_idxs_written, col_idxs_name).wait(); + + Hyperslab const hyperslab{{30, 20}, {20, 10}}; + + IndexArray const row_idxs_read{ + lue::from_gdal(row_idxs_name, hyperslab, partition_shape)}; + IndexArray const col_idxs_read{ + lue::from_gdal(col_idxs_name, hyperslab, partition_shape)}; + + ConditionArray const hyperslab_condition{ + lue::create_partitioned_array(hyperslab.counts(), partition_shape, 1)}; + IndexArray const hyperslab_row_idxs_written = cell_index(hyperslab_condition, 0) + + static_cast(hyperslab.offsets()[0]); + IndexArray const hyperslab_col_idxs_written = cell_index(hyperslab_condition, 1) + + static_cast(hyperslab.offsets()[1]); + + lue::test::check_arrays_are_equal(row_idxs_read, hyperslab_row_idxs_written); + lue::test::check_arrays_are_equal(col_idxs_read, hyperslab_col_idxs_written); +} + + +BOOST_AUTO_TEST_CASE(valid_hyperslab) +{ + using namespace lue::value_policies; + + using ConditionArray = Array; + using IndexArray = Array; + using Shape = lue::ShapeT; + using Hyperslab = lue::Hyperslab<2>; + + Shape const array_shape{60, 40}; + Shape const partition_shape{10, 10}; + + ConditionArray const condition{ + lue::create_partitioned_array(array_shape, partition_shape, 1)}; + IndexArray const idxs_written = cell_index(condition, 0); + + std::string const idxs_name{"lue_framework_io_gdal_wrong_hyperslab_idxs.tif"}; + + lue::to_gdal(idxs_written, idxs_name).wait(); + + { + // Hyperslab corresponds with whole array + Hyperslab const hyperslab{{30, 20}, {60, 40}}; + BOOST_CHECK_NO_THROW(lue::from_gdal(idxs_name, hyperslab, partition_shape)); + } +} + + +BOOST_AUTO_TEST_CASE(incorrect_hyperslab) +{ + using namespace lue::value_policies; + + using ConditionArray = Array; + using IndexArray = Array; + using Shape = lue::ShapeT; + using Hyperslab = lue::Hyperslab<2>; + + Shape const array_shape{60, 40}; + Shape const partition_shape{10, 10}; + + ConditionArray const condition{ + lue::create_partitioned_array(array_shape, partition_shape, 1)}; + IndexArray const idxs_written = cell_index(condition, 0); + + std::string const idxs_name{"lue_framework_io_gdal_wrong_hyperslab_idxs.tif"}; + + lue::to_gdal(idxs_written, idxs_name).wait(); + + { + // Empty hyperslab + Hyperslab const hyperslab{{30, 20}, {0, 0}}; + BOOST_CHECK_EXCEPTION( + lue::from_gdal(idxs_name, hyperslab, partition_shape), + std::runtime_error, + [](auto const& exception) + { return std::string(exception.what()).find("shape must not be empty") != std::string::npos; }); + } + + { + // Empty hyperslab + Hyperslab const hyperslab{{30, 20}, {20, 0}}; + BOOST_CHECK_EXCEPTION( + lue::from_gdal(idxs_name, hyperslab, partition_shape), + std::runtime_error, + [](auto const& exception) + { return std::string(exception.what()).find("shape must not be empty") != std::string::npos; }); + } + + { + // Empty hyperslab + Hyperslab const hyperslab{{30, 20}, {0, 10}}; + BOOST_CHECK_EXCEPTION( + lue::from_gdal(idxs_name, hyperslab, partition_shape), + std::runtime_error, + [](auto const& exception) + { return std::string(exception.what()).find("shape must not be empty") != std::string::npos; }); + } + + { + // Hyperslab too large + Hyperslab const hyperslab{{55, 35}, {20, 20}}; + + BOOST_CHECK_EXCEPTION( + lue::from_gdal(idxs_name, hyperslab, partition_shape), + std::runtime_error, + [](auto const& exception) + { return std::string(exception.what()).find("extents beyond array") != std::string::npos; }); + } +} diff --git a/source/framework/python/CMakeLists.txt b/source/framework/python/CMakeLists.txt index ce95e608a..0e311274a 100644 --- a/source/framework/python/CMakeLists.txt +++ b/source/framework/python/CMakeLists.txt @@ -157,6 +157,7 @@ add_library(lue_py_framework SHARED source/core/create_partitioned_array.cpp ${generated_source_files} + source/core/hyperslab.cpp source/core/scalar.cpp source/core/wait.cpp diff --git a/source/framework/python/source/core/hyperslab.cpp b/source/framework/python/source/core/hyperslab.cpp new file mode 100644 index 000000000..4367b7c4d --- /dev/null +++ b/source/framework/python/source/core/hyperslab.cpp @@ -0,0 +1,48 @@ +#include "lue/framework/core/hyperslab.hpp" +#include +#include + + +using namespace pybind11::literals; + + +namespace lue::framework { + namespace { + + auto to_center(std::tuple const& center) + { + return Hyperslab<2>::Center{std::get<0>(center), std::get<1>(center)}; + }; + + auto to_shape(std::tuple const& shape) + { + return Hyperslab<2>::Shape{std::get<0>(shape), std::get<1>(shape)}; + }; + + } // Anonymous namespace + + + void bind_hyperslab(pybind11::module& module) + { + Rank const rank{2}; + using Hyperslab = lue::Hyperslab; + + pybind11::class_(module, "Hyperslab") + .def( + pybind11::init( + [](std::tuple const& center, std::tuple const& shape) + { return Hyperslab{to_center(center), to_shape(shape)}; }), + R"( + Construct a hyperslab given a @a center cell and a @a shape + + :param center: Indices of center cell + :param shape: Number of cells in each dimension +)", + "center"_a, + "shape"_a) + .def_property_readonly("offsets", &Hyperslab::offsets) + .def_property_readonly("counts", &Hyperslab::counts) + .def_property_readonly("strides", &Hyperslab::strides); + } + +} // namespace lue::framework diff --git a/source/framework/python/source/gdal/from_gdal.cpp b/source/framework/python/source/gdal/from_gdal.cpp index c2526c0e7..086c4667c 100644 --- a/source/framework/python/source/gdal/from_gdal.cpp +++ b/source/framework/python/source/gdal/from_gdal.cpp @@ -20,7 +20,127 @@ namespace lue::framework { } - auto from_gdal_py(std::string const& name, std::optional const& partition_shape) + template + auto from_gdal( + std::string const& name, + Hyperslab<2> const& hyperslab, + StaticShape<2> const& partition_shape) -> pybind11::object + { + return pybind11::cast(lue::from_gdal(name, hyperslab, partition_shape)); + } + + + auto hyperslab_from_gdal( + std::string const& name, + Hyperslab<2> const& hyperslab, + std::optional const& partition_shape) -> pybind11::object + { + gdal::DatasetPtr dataset{gdal::open_dataset(name, GDALAccess::GA_ReadOnly)}; + auto const raster_shape{gdal::shape(*dataset)}; + StaticShape<2> static_array_shape{raster_shape[0], raster_shape[1]}; + StaticShape<2> static_partition_shape{}; + + if (partition_shape) + { + DynamicShape const dynamic_partition_shape{tuple_to_shape(*partition_shape)}; + + if (static_array_shape.size() != dynamic_partition_shape.size()) + { + throw std::runtime_error(fmt::format( + "Rank of array shape and partition shape must be equal ({} != {})", + static_array_shape.size(), + dynamic_partition_shape.size())); + } + + static_partition_shape = dynamic_shape_to_static_shape<2>(dynamic_partition_shape); + } + else + { + static_partition_shape = default_partition_shape(hyperslab.counts()); + } + + pybind11::object result; + GDALDataType const data_type{gdal::data_type(*dataset)}; + + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_Byte) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } + +#if LUE_GDAL_SUPPORTS_8BIT_SIGNED_INTEGERS == 1 + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_Int8) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } +#endif + + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_UInt32) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } + + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_Int32) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } + +#if LUE_GDAL_SUPPORTS_64BIT_INTEGERS == 1 + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_UInt64) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } + + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_Int64) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } +#endif + + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_Float32) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } + + if constexpr (lue::arithmetic_element_supported) + { + if (data_type == GDT_Float64) + { + result = framework::from_gdal(name, hyperslab, static_partition_shape); + } + } + + if (!result) + { + throw std::runtime_error( + fmt::format("Unsupported datatype: {}", "TODO")); // as_string(data_type))); + } + + return result; + } + + + auto array_from_gdal(std::string const& name, std::optional const& partition_shape) -> pybind11::object { gdal::DatasetPtr dataset{gdal::open_dataset(name, GDALAccess::GA_ReadOnly)}; @@ -149,11 +269,29 @@ namespace lue::framework { { module.def( "from_gdal", - from_gdal_py, + array_from_gdal, + R"( + Create new array, filled with values read using GDAL + + :param str name: Name of dataset to read + :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 +)", + "name"_a, + pybind11::kw_only(), + "partition_shape"_a = std::optional{}, + pybind11::return_value_policy::move); + + module.def( + "from_gdal", + hyperslab_from_gdal, R"( Create new array, filled with values read using GDAL :param str name: Name of dataset to read + :param hyperslab: Hyperslab of cells to read :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. @@ -161,6 +299,7 @@ namespace lue::framework { )", "name"_a, pybind11::kw_only(), + "hyperslab"_a, "partition_shape"_a = std::optional{}, pybind11::return_value_policy::move); diff --git a/source/framework/python/source/submodule.cpp b/source/framework/python/source/submodule.cpp index 1a5d4de42..2fd32ea11 100644 --- a/source/framework/python/source/submodule.cpp +++ b/source/framework/python/source/submodule.cpp @@ -67,6 +67,7 @@ namespace lue::framework { void bind_hpx(pybind11::module& module); void bind_create_array(pybind11::module& module); + void bind_hyperslab(pybind11::module& module); void bind_wait_partitioned_array(pybind11::module& module); void bind_focal_operations(pybind11::module& module); void bind_global_operations(pybind11::module& module); @@ -142,6 +143,7 @@ namespace lue::framework { // Wrap high-level algorithms bind_create_array(submodule); + bind_hyperslab(submodule); bind_wait_partitioned_array(submodule); bind_focal_operations(submodule); bind_global_operations(submodule); diff --git a/source/framework/python/test/gdal_test.py b/source/framework/python/test/gdal_test.py index 28a36788e..1e57ab761 100644 --- a/source/framework/python/test/gdal_test.py +++ b/source/framework/python/test/gdal_test.py @@ -12,17 +12,16 @@ def tearDownModule(): lue_test.stop_hpx_runtime() -class ToNumPyTest(lue_test.TestCase): +class GdalTest(lue_test.TestCase): @lue_test.framework_test_case - def test_round_trip1(self): - # Use default partition shape + def test_array_default_partition_shape(self): array_shape = (600, 400) dtype = np.int32 fill_value = 5 array_written = lfr.create_array(array_shape, dtype, fill_value) - name = "to_gdal_round_trip1.tif" + name = "gdal_array_default_partition_shape.tif" written = lfr.to_gdal(array_written, name) written.wait() # Caveat! Don't read before the raster is written! @@ -34,8 +33,7 @@ def test_round_trip1(self): self.assertTrue(lfr.all(array_read == array_written).future.get()) @lue_test.framework_test_case - def test_round_trip2(self): - # Use specific partition shape + def test_array_custom_partition_shape(self): array_shape = (600, 400) partition_shape = (10, 10) @@ -45,7 +43,7 @@ def test_round_trip2(self): array_shape, dtype, fill_value, partition_shape=partition_shape ) - name = "to_gdal_round_trip2.tif" + name = "gdal_array_custom_partition_shape.tif" written = lfr.to_gdal(array_written, name) written.wait() # Caveat! Don't read before the raster is written! @@ -55,3 +53,41 @@ def test_round_trip2(self): self.assertEqual(array_read.dtype, array_written.dtype) self.assertEqual(array_read.shape, array_written.shape) self.assertTrue(lfr.all(array_read == array_written).future.get()) + + @lue_test.framework_test_case + def test_hyperslab(self): + + array_shape = (600, 400) + + condition = lfr.create_array(array_shape, lfr.boolean_element_type, 1) + row_idxs_written = lfr.cell_index(condition, 0) + col_idxs_written = lfr.cell_index(condition, 1) + + row_idxs_name = "gdal_hyperslab_row_idxs.tif" + written = lfr.to_gdal(row_idxs_written, row_idxs_name) + written.wait() + + col_idxs_name = "gdal_hyperslab_col_idxs.tif" + written = lfr.to_gdal(col_idxs_written, col_idxs_name) + written.wait() + + hyperslab_shape = (200, 100) + hyperslab = lfr.Hyperslab(center=(300, 200), shape=hyperslab_shape) + + row_idxs_read = lfr.from_gdal(row_idxs_name, hyperslab=hyperslab) + col_idxs_read = lfr.from_gdal(col_idxs_name, hyperslab=hyperslab) + + self.assertEqual(row_idxs_read.dtype, row_idxs_written.dtype) + self.assertEqual(col_idxs_read.dtype, col_idxs_written.dtype) + + self.assertEqual(row_idxs_read.shape, hyperslab_shape) + self.assertEqual(col_idxs_read.shape, hyperslab_shape) + + # Create arrays that contain the same values as the ones in the hyperslabs read from the raster + condition = lfr.create_array(hyperslab_shape, lfr.boolean_element_type, 1) + row_idxs_hyperslab = lfr.cell_index(condition, 0) + hyperslab.offsets[0] + col_idxs_hyperslab = lfr.cell_index(condition, 1) + hyperslab.offsets[1] + + # Verify the hyperslabs read contain values we expect + self.assertTrue(lfr.all(row_idxs_read == row_idxs_hyperslab).future.get()) + self.assertTrue(lfr.all(col_idxs_read == col_idxs_hyperslab).future.get())