diff --git a/algorithms/test/CMakeLists.txt b/algorithms/test/CMakeLists.txt
index 52948e1..6250e96 100644
--- a/algorithms/test/CMakeLists.txt
+++ b/algorithms/test/CMakeLists.txt
@@ -20,6 +20,7 @@ add_executable(gdxalgotest
     maximumtest.cpp
     nodatatest.cpp
     normalisetest.cpp
+    $<$<BOOL:${GDX_ENABLE_GEOMETRY}>:polygoncoveragetest.cpp>
     propdisttest.cpp
     rasterizetest.cpp
     rasterizelineantialiasedtest.cpp
diff --git a/algorithms/test/data/boundaries.gpkg b/algorithms/test/data/boundaries.gpkg
new file mode 100644
index 0000000..9124d7a
Binary files /dev/null and b/algorithms/test/data/boundaries.gpkg differ
diff --git a/algorithms/test/polygoncoveragetest.cpp b/algorithms/test/polygoncoveragetest.cpp
new file mode 100644
index 0000000..7d8d58f
--- /dev/null
+++ b/algorithms/test/polygoncoveragetest.cpp
@@ -0,0 +1,53 @@
+#include "gdx/algo/polygoncoverage.h"
+#include "gdx/test/testbase.h"
+#include "infra/crs.h"
+#include "infra/gdalalgo.h"
+#include "infra/gdalio.h"
+
+#include "testconfig.h"
+
+namespace gdx::test {
+
+using namespace inf;
+namespace gdal = inf::gdal;
+
+TEST_CASE("PolygonCoverage")
+{
+    auto boundaries = file::u8path(TEST_DATA_DIR) / "boundaries.gpkg";
+
+    auto ds = gdal::VectorDataSet::open(boundaries);
+    GeoMetadata outputExtent(120, 260, 11000.0, 140000.0, {1000.0, -1000.0}, std::numeric_limits<double>::quiet_NaN(), gdal::SpatialReference(crs::epsg::BelgianLambert72).export_to_wkt());
+
+    auto warpedMeta      = gdal::warp_metadata(outputExtent, crs::epsg::WGS84);
+    const auto coverages = create_polygon_coverages(warpedMeta, ds, gdx::BorderHandling::None, 1.0, {}, {}, "Code3", nullptr);
+
+    CHECK(coverages.size() == 3);
+
+    for (auto& coverage : coverages) {
+        auto totalCoverage = std::accumulate(coverages[0].cells.begin(), coverages[0].cells.end(), 0.0, [](double sum, const PolygonCellCoverage::CellInfo& cellInfo) {
+            return sum + cellInfo.coverage;
+        });
+
+        if (coverage.name == "BEB") {
+            CHECK(coverage.cells.size() == 145);
+
+            auto cellIter = std::find_if(coverage.cells.begin(), coverage.cells.end(), [](auto& c) { return c.computeGridCell == Cell(55, 147); });
+            CHECK(cellIter->cellCoverage == Approx(0.6037847694229548));
+
+        } else if (coverage.name == "BEF") {
+            CHECK(coverage.cells.size() == 10053);
+
+            auto cellIter = std::find_if(coverage.cells.begin(), coverage.cells.end(), [](auto& c) { return c.computeGridCell == Cell(55, 147); });
+            CHECK(cellIter->cellCoverage == Approx(0.3962152305751532));
+        } else if (coverage.name == "NL") {
+            CHECK(coverage.cells.size() == 28072);
+            // This cell does not overlap NL
+            CHECK(std::find_if(coverage.cells.begin(), coverage.cells.end(), [](auto& c) { return c.computeGridCell == Cell(55, 147); }) == coverage.cells.end());
+        } else {
+            CHECK_FALSE_MESSAGE("Unexpected polygon name", coverage.name);
+        }
+
+        CHECK_MESSAGE(totalCoverage == Approx(1.0), "Coverages for ", coverage.name, " do not add up to 1");
+    }
+}
+}
diff --git a/core/test/rasteriotest.cpp b/core/test/rasteriotest.cpp
index 22e6842..96e8f4d 100644
--- a/core/test/rasteriotest.cpp
+++ b/core/test/rasteriotest.cpp
@@ -72,7 +72,7 @@ TEST_CASE("read metadata")
 
 TEST_CASE("read metadata with projection")
 {
-    auto meta = inf::gdal::io::read_metadata(fs::u8path(TEST_DATA_DIR) / "../../../test/mapdata/landusebyte.tif");
+    auto meta = inf::gdal::io::read_metadata(file::u8path(TEST_DATA_DIR) / "../../../test/mapdata/landusebyte.tif");
     CHECK(meta.projection_frienly_name() == "EPSG:31370");
 }
 
@@ -286,7 +286,7 @@ template <typename T>
 void testRaster(std::string_view filename)
 {
     MaskedRaster<int32_t> referenceRaster;
-    referenceRaster.set_metadata(read_raster(fs::u8path(TEST_DATA_DIR) / std::string(filename), referenceRaster));
+    referenceRaster.set_metadata(read_raster(file::u8path(TEST_DATA_DIR) / std::string(filename), referenceRaster));
     write_raster(referenceRaster, "raster.asc");
 
     MaskedRaster<int32_t> writtenRaster;
@@ -376,18 +376,18 @@ TEST_CASE("write raster with nodata float")
 
     std::vector<std::string> expected({
         "ncols        5"s,
-            "nrows        3"s,
-            "xllcorner    1.000000000000"s,
-            "yllcorner    -10.000000000000"s,
-            "cellsize     4.000000000000"s,
+        "nrows        3"s,
+        "xllcorner    1.000000000000"s,
+        "yllcorner    -10.000000000000"s,
+        "cellsize     4.000000000000"s,
 #if GDAL_VERSION_NUM < GDAL_COMPUTE_VERSION(3, 8, 0)
-            "NODATA_value  -1"s,
+        "NODATA_value  -1"s,
 #else
-            "NODATA_value -1"s,
+        "NODATA_value -1"s,
 #endif
-            "0.0 1 2 3 4"s,
-            "5 6 7 8 9"s,
-            "4 3 2 -1 0"s,
+        "0.0 1 2 3 4"s,
+        "5 6 7 8 9"s,
+        "4 3 2 -1 0"s,
     });
 
     std::ifstream str("raster.asc");
@@ -612,7 +612,7 @@ TEST_CASE("Write raster with different data type")
 
     auto path = fs::temp_directory_path() / "rasterio" / "int32.tif";
 
-    auto ras = read_masked_raster<uint32_t>(fs::u8path(TEST_DATA_DIR) / "testraster.asc", extent);
+    auto ras = read_masked_raster<uint32_t>(file::u8path(TEST_DATA_DIR) / "testraster.asc", extent);
     gdx::write_raster(ras, path, typeid(int32_t));
 
     CHECK(inf::type_name(gdal::io::get_raster_type(path)) == inf::type_name(typeid(int32_t)));