ADD georeferenced GeoTIFF export via to_raster()

docs/spectral-indices.md

@@ -76,6 +76,42 @@ SAVICalculator(
 
 Sentinel-2 processing baseline 04.00 and later carries a `BOA_ADD_OFFSET` of $-1000$, hence the separate preset. If your product is already reflectance, the defaults are correct and nothing needs passing.
 
+## Output: Picture or Data Product
+
+Every tool can write its result two ways, and the difference matters.
+
+| | `execute()` | `to_raster()` |
+|---|---|---|
+| Format | PNG | GeoTIFF |
+| Values | quantized to 256 levels per channel | full precision (`float32` / `int32`) |
+| Pixel grid | resampled by `dpi` and `bbox_inches` | identical to the input |
+| CRS / transform | none | copied from the source band |
+| Use | inspection, reporting, figures | GIS overlay, zonal statistics, differencing |
+
+`execute()` renders through matplotlib, with margins, colorbar and watermark baked in, so the result is a **picture of** the computed array rather than the array. A pixel that was `0.6237` cannot be recovered from it, which rules out thresholding, zonal statistics over mapped units, and multitemporal differencing. Use it for figures.
+
+`to_raster()` writes the array itself:
+
+```Python
+from fezrs import NDVICalculator
+
+calculator = NDVICalculator(nir_path="nir.tif", red_path="red.tif")
+calculator.process()
+calculator.to_raster("./exports/ndvi.tif")
+```
+
+The output carries the CRS and affine transform of the source band, so it lands in the right place when opened over other layers in QGIS or ArcGIS. Multi-component results, such as PCA's `(6, height, width)`, are written as multi-band rasters.
+
+**Defaults**, chosen for scene-scale raster products:
+
+- `float32` for continuous results — well beyond reflectance precision, half the size of `float64` — and `int32` for integer label maps such as classifications.
+- `nodata = NaN` for floating point output.
+- `tiled=True`, `compress="deflate"` with `predictor=3` for float data, and `BIGTIFF="IF_SAFER"`.
+
+Override any of them with `dtype=`, `nodata=` and `compress=`.
+
+**A source without a CRS raises.** Writing an identity transform instead would produce a file that looks georeferenced while placing the scene at the coordinate origin. If the inputs carry no spatial referencing, `execute()` is the appropriate output.
+
 ## Mathematical & Scientific Formulations
 
 ### `NDVICalculator` (Normalized Difference Vegetation Index)

fezrs/base.py

@@ -1,5 +1,8 @@
 # Import packages and libraries
 from abc import ABC
+
+import numpy as np
+import rasterio as rio
 from PIL import Image
 from pathlib import Path
 from uuid import uuid4
@@ -140,6 +143,102 @@ def _export_file(
         plt.close(fig)
         return filename
 
+    def to_raster(
+        self,
+        output_path: BandPathType,
+        reference_band: str | None = None,
+        dtype: str | None = None,
+        nodata=None,
+        compress: str = "deflate",
+    ):
+        """
+        Write the computed result as a georeferenced GeoTIFF.
+
+        ``execute()`` renders through matplotlib and saves a PNG: values are
+        quantized to 256 levels per channel, the pixel grid is resampled by dpi
+        and bbox_inches, and CRS and transform are discarded. That output is a
+        picture of the result, not the result. This method writes the array
+        itself, carrying the CRS and affine transform of the source band, so it
+        can be overlaid in a GIS, intersected with mapped units, differenced
+        against another date, or used for zonal statistics.
+
+        Multi-component outputs, such as PCA's ``(6, height, width)``, are
+        written as multi-band rasters.
+
+        Args:
+            output_path: Destination ``.tif`` path. Parent directories are
+                created.
+            reference_band: Band whose spatial referencing to copy. Defaults to
+                the first supplied band.
+            dtype: Output dtype. Defaults to ``float32`` for continuous results
+                and ``int32`` for integer label maps such as classifications.
+                float32 is well beyond reflectance precision and half the size
+                of float64.
+            nodata: Nodata value. Defaults to NaN for floating point output.
+            compress: GeoTIFF compression. ``deflate`` with a horizontal
+                differencing predictor is the usual choice for float rasters.
+
+        Returns:
+            str: Path to the written raster.
+
+        Raises:
+            ValueError: If nothing has been computed, or the source carried no
+                spatial referencing to propagate.
+        """
+        if self._output is None:
+            raise ValueError("Data not computed.")
+
+        profile = self.files_handler.get_raster_profile(reference_band)
+
+        if profile is None or profile["crs"] is None:
+            raise ValueError(
+                "The source band carries no CRS, so the result cannot be "
+                "written as a georeferenced raster. Reading the inputs from "
+                "GeoTIFFs that declare a CRS and transform will fix this. Use "
+                "execute() if a plain image is what you want."
+            )
+
+        array = np.asarray(self._output)
+        if array.ndim == 2:
+            array = array[np.newaxis, :, :]
+        elif array.ndim != 3:
+            raise ValueError(
+                f"Cannot write an array with {array.ndim} dimensions as a raster."
+            )
+
+        if dtype is None:
+            dtype = (
+                "int32" if np.issubdtype(array.dtype, np.integer) else "float32"
+            )
+
+        if nodata is None and not np.issubdtype(np.dtype(dtype), np.integer):
+            nodata = float("nan")
+
+        output_path = Path(output_path)
+        output_path.parent.mkdir(parents=True, exist_ok=True)
+
+        creation = {
+            "driver": "GTiff",
+            "height": array.shape[1],
+            "width": array.shape[2],
+            "count": array.shape[0],
+            "dtype": dtype,
+            "crs": profile["crs"],
+            "transform": profile["transform"],
+            "nodata": nodata,
+            "tiled": True,
+            "compress": compress,
+            "BIGTIFF": "IF_SAFER",
+        }
+        if compress == "deflate" and not np.issubdtype(np.dtype(dtype), np.integer):
+            # Horizontal differencing for floating point data.
+            creation["predictor"] = 3
+
+        with rio.open(output_path, "w", **creation) as destination:
+            destination.write(array.astype(dtype))
+
+        return str(output_path)
+
     def execute(
         self,
         output_path: BandPathType,

fezrs/utils/file_handler.py

@@ -96,6 +96,27 @@ def _rasterio_image_tifs(path: str):
     return rio.open(path)
 
 
+def _raster_profile(path: str) -> Dict:
+    """
+    Read the spatial referencing of a raster without loading its pixels.
+
+    Args:
+        path (str): Path to a raster file.
+
+    Returns:
+        Dict: CRS, affine transform, nodata value, dtype and shape.
+    """
+    with rio.open(path) as source:
+        return {
+            "crs": source.crs,
+            "transform": source.transform,
+            "nodata": source.nodata,
+            "dtype": source.dtypes[0],
+            "height": source.height,
+            "width": source.width,
+        }
+
+
 class FileHandler:
     """
     FileHandler is a utility class for managing and processing geospatial image files.
@@ -225,6 +246,42 @@ def get_bands(self, requested_bands: Optional[List[BandNameType]] = None):
             if self.bands.get(band) is not None
         }
 
+    def get_raster_profile(
+        self, band: Optional[BandNameType] = None
+    ) -> Optional[Dict]:
+        """
+        Retrieve the spatial referencing of an input band.
+
+        Bands are loaded for computation through scikit-image, which discards
+        CRS, transform and nodata. This reads that metadata back via rasterio so
+        a result can be written out as a georeferenced raster rather than only
+        as a picture of one.
+
+        Args:
+            band (Optional[BandNameType]): Band to describe. Defaults to the
+                first band that was supplied.
+
+        Returns:
+            Optional[Dict]: Profile mapping, or None when no source is available.
+        """
+        if band is None:
+            candidates = [
+                name
+                for name, path in self.band_paths.items()
+                if path is not None
+            ]
+            if not candidates and self.tif_paths:
+                return _raster_profile(str(self.tif_paths[0]))
+            if not candidates:
+                return None
+            band = candidates[0]
+
+        path = self.band_paths.get(band)
+        if path is None or not os.path.exists(path):
+            return None
+
+        return _raster_profile(str(path))
+
     def get_metadata_bands(
         self, requested_bands: Optional[list[BandNameType]] = None
     ) -> Dict[str, Dict]:

paper/paper.md

@@ -108,7 +108,7 @@ FEZrs was built to address this fragmentation by consolidating frequently used r
 
 FEZrs is organized around a shared abstract base class (`BaseTool`) that standardizes the lifecycle of every analysis module: input validation, processing, visualization customization, and export. Each calculator (for example, `NDVICalculator`, `GaussianCalculator`, or `KMeansCalculator`) inherits this interface, accepts band file paths through a common file-handling layer, and exposes a small set of methods such as `execute` and `histogram_export`. This pattern keeps the public API uniform across spectral indices, filters, enhancement tools, change-detection utilities, and machine-learning modules, reducing cognitive overhead when composing multi-step workflows.
 
-The package is modular by domain: (`spectral_indices`, `filters`, `image_enhancement`, `change_detection`, `clustering`, `glcm`, `pca`, `svm`, and others), while shared utilities handle band-path typing, file I/O, and histogram support. Raster inputs are typically multi-band geospatial imagery. Outputs are written as figures and processed products suitable for inspection, reporting, and further analysis. The design deliberately favors composition of independent calculators over a single monolithic pipeline object, so researchers can select only the methods required for a given study while still benefiting from consistent validation and export behavior.
+The package is modular by domain: (`spectral_indices`, `filters`, `image_enhancement`, `change_detection`, `clustering`, `glcm`, `pca`, `svm`, and others), while shared utilities handle band-path typing, file I/O, and histogram support. Raster inputs are typically multi-band geospatial imagery. Results can be exported two ways: as rendered figures via `execute()`, for inspection and reporting, and as georeferenced GeoTIFF rasters via `to_raster()`, which preserve full numerical precision along with the coordinate reference system and affine transform of the source imagery, and are therefore suitable for GIS overlay, zonal statistics, and multitemporal analysis. The design deliberately favors composition of independent calculators over a single monolithic pipeline object, so researchers can select only the methods required for a given study while still benefiting from consistent validation and export behavior.
 
 # Research Impact Statement