Merge pull request #122 from quantifyearth/mwd-sparse-geotiff

Add support for sparse geotiffs

Author: mdales

CHANGES.md

@@ -3,6 +3,7 @@
 ### Added
 
 * Added a `to_graphviz` call for expressions to help with debugging.
+* Support sparse GeoTIFF creation.
 
 ### Fixed
 

pyproject.toml

@@ -55,6 +55,7 @@ dev = [
     "mike",
     "mkdocs-gen-files",
     "black",
+    "tifffile",
 ]
 
 [project.urls]

tests/unit/test_sparse_files.py

@@ -0,0 +1,95 @@
+import tempfile
+from pathlib import Path
+
+import numpy as np
+import pytest
+import tifffile
+import yirgacheffe as yg
+
+def test_sprase_needs_nodata() -> None:
+    data = np.zeros((5000, 5000))
+    with tempfile.TemporaryDirectory() as tempdir:
+        filename = Path(tempdir) / "test.tif"
+        with yg.from_array(data, (0, 0), yg.MapProjection("EPSG:4326", 10.0, -10.0)) as layer:
+            with pytest.raises(ValueError):
+                layer.to_geotiff(filename, sparse=True)
+
+
+@pytest.mark.parametrize("dtype,zero", [
+    (np.int8, 0),
+    (np.float32, 0.0),
+])
+def test_all_sparse(dtype: type, zero: int | float) -> None:
+    data = np.zeros((5000, 5000)).astype(dtype) # type: ignore
+    with tempfile.TemporaryDirectory() as tempdir:
+        filename = Path(tempdir) / "test.tif"
+        with yg.from_array(data, (0, 0), yg.MapProjection("EPSG:4326", 10.0, -10.0)) as layer:
+            layer.to_geotiff(filename, nodata=zero, sparse=True)
+        with yg.read_raster(filename) as layer:
+            # no data is nan, so convert that to 0, as otherwise
+            # this test will miss things
+            assert layer.nan_to_num().sum() == 0
+
+        # now look into the tiff structure to see that all blocks are tiff
+        with tifffile.TiffFile(filename) as tif:
+            page = tif.pages[0] # type: ignore
+            offsets = page.tags.get('StripOffsets').value # type: ignore
+            byte_counts = page.tags.get('StripByteCounts').value # type: ignore
+
+            for offset, count in zip(offsets, byte_counts):
+                assert offset == 0
+                assert count == 0
+
+
+@pytest.mark.parametrize("dtype,zero", [
+    (np.int8, 0),
+    (np.float32, 0.0),
+])
+def test_none_sparse(dtype: type, zero: int | float) -> None:
+    data = np.ones((5000, 5000)).astype(dtype)  # type: ignore
+    with tempfile.TemporaryDirectory() as tempdir:
+        filename = Path(tempdir) / "test.tif"
+        with yg.from_array(data, (0, 0), yg.MapProjection("EPSG:4326", 10.0, -10.0)) as layer:
+            layer.to_geotiff(filename, nodata=zero, sparse=True)
+        with yg.read_raster(filename) as layer:
+            assert layer.sum() == 5000 * 5000
+
+        # now look into the tiff structure to see that all blocks are tiff
+        with tifffile.TiffFile(filename) as tif:
+            page = tif.pages[0] # type: ignore
+            offsets = page.tags.get('StripOffsets').value # type: ignore
+            byte_counts = page.tags.get('StripByteCounts').value # type: ignore
+
+            for offset, count in zip(offsets, byte_counts):
+                assert offset != 0
+                assert count != 0
+
+
+@pytest.mark.parametrize("dtype,zero", [
+    (np.int8, 0),
+    (np.float32, 0.0),
+])
+def test_mixed_sparse(dtype: type, zero: int | float) -> None:
+    data = np.array([[i % 2] * 5000 for i in range(5000)]).astype(dtype)  # type: ignore
+    with tempfile.TemporaryDirectory() as tempdir:
+        filename = Path(tempdir) / "test.tif"
+        with yg.from_array(data, (0, 0), yg.MapProjection("EPSG:4326", 10.0, -10.0)) as layer:
+            layer.to_geotiff(filename, nodata=zero, sparse=True)
+        with yg.read_raster(filename) as layer:
+            # no data is nan, so convert that to 0, as otherwise
+            # this test will miss things
+            assert layer.nan_to_num().sum() == 5000 * 5000 / 2
+
+        # now look into the tiff structure to see that all blocks are tiff
+        with tifffile.TiffFile(filename) as tif:
+            page = tif.pages[0] # type: ignore
+            offsets = page.tags.get('StripOffsets').value # type: ignore
+            byte_counts = page.tags.get('StripByteCounts').value # type: ignore
+
+            for i, (offset, count) in enumerate(zip(offsets, byte_counts)):
+                if i % 2:
+                    assert offset != 0
+                    assert count != 0
+                else:
+                    assert offset == 0
+                    assert count == 0

yirgacheffe/_operators/__init__.py

@@ -312,6 +312,7 @@ def to_geotiff(
         parallelism: int | bool | None = None,
         callback: Callable[[float], None] | None = None,
         nodata: float | int | None = None,
+        sparse: bool = False,
     ) -> float | None:
         """Saves a calculation to a raster file, optionally also returning the sum of pixels.
 
@@ -323,7 +324,8 @@ def to_geotiff(
                 yirgacheffe will pick a sensible value.
             callback: If passed, this callback will be called periodically with a progress update for the saving,
                 with a value between 0.0 and 1.0.
-            nodata: Nominate a value to be stored as nodata in the result
+            nodata: Nominate a value to be stored as nodata in the result.
+            sparse: If True then save a sparse GeoTIFF as per GDAL's extension to the GeoTIFF standard.
 
         Returns:
             Either returns None, or the sum of the pixels in the resulting raster if `and_sum` was specified.
@@ -334,6 +336,7 @@ def to_geotiff(
             parallelism=parallelism,
             callback=callback,
             nodata=nodata,
+            sparse=sparse,
         )
 
     def sum(self):
@@ -1328,8 +1331,12 @@ def to_geotiff(
         parallelism: int | bool | None = None,
         callback: Callable[[float], None] | None = None,
         nodata: float | int | None = None,
+        sparse: bool = False,
     ) -> float | None:
 
+        if sparse and nodata is None:
+            raise ValueError("Nodata value must be provided for sparse GeoTIFFs")
+
         # We want to write to a tempfile before we move the result into place, but we can't use
         # the actual $TMPDIR as that might be on a different device, and so we use a file next to where
         # the final file will be, so we just need to rename the file at the end, not move it. But for special cases,
@@ -1351,7 +1358,8 @@ def to_geotiff(
             with RasterLayer.empty_raster_layer_like(
                 self,
                 filename=inner_filename, # type: ignore
-                nodata=nodata
+                nodata=nodata,
+                sparse=sparse,
             ) as layer:
                 if parallelism is None:
                     result = self.save(layer, and_sum=and_sum, callback=callback)

yirgacheffe/layers/rasters.py

@@ -33,7 +33,8 @@ def empty_raster_layer(
         nodata: float | int | None = None,
         nbits: int | None = None,
         threads: int | None = None,
-        bands: int=1
+        bands: int=1,
+        sparse: bool=False,
     ) -> RasterLayer:
         abs_xstep, abs_ystep = abs(scale.xstep), abs(scale.ystep)
 
@@ -48,6 +49,8 @@ def empty_raster_layer(
             options.append(f"NUM_THREADS={threads}")
         if nbits is not None:
             options.append(f"NBITS={nbits}")
+        if sparse:
+            options.append("SPARSE_OK=YES")
 
         map_projection = MapProjection(projection, scale.xstep, scale.ystep)
         # If the area is projected, use that, otherwise we need to project it
@@ -100,7 +103,8 @@ def empty_raster_layer_like(
         nodata: float | int | None = None,
         nbits: int | None = None,
         threads: int | None = None,
-        bands: int=1
+        bands: int=1,
+        sparse: bool=False,
     ) -> RasterLayer:
         if area is None:
             area = layer.area
@@ -135,6 +139,8 @@ def empty_raster_layer_like(
             options.append(f"NUM_THREADS={threads}")
         if nbits is not None:
             options.append(f"NBITS={nbits}")
+        if sparse:
+            options.append("SPARSE_OK=YES")
 
         if filename:
             driver = gdal.GetDriverByName('GTiff')