Rename GeoTIFF API to xarray conventions (#1047)

open_geotiff replaces read_geotiff, to_geotiff replaces write_geotiff.
Adds .xrs.to_geotiff() accessor on DataArray and Dataset, and
.xrs.open_geotiff() on Dataset for spatially-windowed reads.

Author: brendancol

xrspatial/accessor.py

@@ -26,14 +26,14 @@ def __init__(self, obj):
     def plot(self, **kwargs):
         """Plot the DataArray, using an embedded TIFF colormap if present.
 
-        For palette/indexed-color GeoTIFFs (read via ``read_geotiff``),
+        For palette/indexed-color GeoTIFFs (read via ``open_geotiff``),
         the TIFF's color table is applied automatically with correct
         normalization.  For all other DataArrays, falls through to the
         standard ``da.plot()``.
 
         Usage::
 
-            da = read_geotiff('landcover.tif')
+            da = open_geotiff('landcover.tif')
             da.xrs.plot()  # palette colors used automatically
         """
         import numpy as np
@@ -460,6 +460,18 @@ def rasterize(self, geometries, **kwargs):
         from .rasterize import rasterize
         return rasterize(geometries, like=self._obj, **kwargs)
 
+    # ---- GeoTIFF I/O ----
+
+    def to_geotiff(self, path, **kwargs):
+        """Write this DataArray as a GeoTIFF.
+
+        Equivalent to ``to_geotiff(da, path, **kwargs)``.
+
+        See :func:`xrspatial.geotiff.to_geotiff` for full parameter docs.
+        """
+        from .geotiff import to_geotiff
+        return to_geotiff(self._obj, path, **kwargs)
+
 
 @xr.register_dataset_accessor("xrs")
 class XrsSpatialDatasetAccessor:
@@ -776,3 +788,75 @@ def rasterize(self, geometries, **kwargs):
             "Dataset has no 2D variable with 'y' and 'x' dimensions "
             "to use as rasterize template"
         )
+
+    # ---- GeoTIFF I/O ----
+
+    def to_geotiff(self, path, var=None, **kwargs):
+        """Write a Dataset variable as a GeoTIFF.
+
+        Parameters
+        ----------
+        path : str
+            Output file path.
+        var : str or None
+            Variable name to write.  If None, uses the first 2D variable
+            with y/x dimensions.
+        **kwargs
+            Passed to :func:`xrspatial.geotiff.to_geotiff`.
+        """
+        from .geotiff import to_geotiff
+        ds = self._obj
+        if var is not None:
+            return to_geotiff(ds[var], path, **kwargs)
+        for v in ds.data_vars:
+            da = ds[v]
+            if da.ndim >= 2 and 'y' in da.dims and 'x' in da.dims:
+                return to_geotiff(da, path, **kwargs)
+        raise ValueError(
+            "Dataset has no variable with 'y' and 'x' dimensions to write"
+        )
+
+    def open_geotiff(self, source, **kwargs):
+        """Read a GeoTIFF windowed to this Dataset's spatial extent.
+
+        Uses the Dataset's y/x coordinates to compute a pixel window,
+        then reads only that region from the file.
+
+        Parameters
+        ----------
+        source : str
+            File path to the GeoTIFF.
+        **kwargs
+            Passed to :func:`xrspatial.geotiff.open_geotiff` (except
+            ``window``, which is computed automatically).
+
+        Returns
+        -------
+        xr.DataArray
+            The windowed portion of the GeoTIFF.
+        """
+        from .geotiff import open_geotiff, _read_geo_info, _extent_to_window
+        ds = self._obj
+        if 'y' not in ds.coords or 'x' not in ds.coords:
+            raise ValueError(
+                "Dataset must have 'y' and 'x' coordinates to compute "
+                "a spatial window"
+            )
+        y = ds.coords['y'].values
+        x = ds.coords['x'].values
+        y_min, y_max = float(y.min()), float(y.max())
+        x_min, x_max = float(x.min()), float(x.max())
+
+        geo_info, file_h, file_w = _read_geo_info(source)
+        t = geo_info.transform
+
+        # Expand extent by half a pixel so we capture edge pixels
+        y_min -= abs(t.pixel_height) * 0.5
+        y_max += abs(t.pixel_height) * 0.5
+        x_min -= abs(t.pixel_width) * 0.5
+        x_max += abs(t.pixel_width) * 0.5
+
+        window = _extent_to_window(t, file_h, file_w,
+                                   y_min, y_max, x_min, x_max)
+        kwargs.pop('window', None)
+        return open_geotiff(source, window=window, **kwargs)

xrspatial/geotiff/__init__.py

@@ -4,12 +4,12 @@
 
 Public API
 ----------
-read_geotiff(source, ...)
+open_geotiff(source, ...)
     Read a GeoTIFF file to an xarray.DataArray.
-write_geotiff(data, path, ...)
+to_geotiff(data, path, ...)
     Write an xarray.DataArray as a GeoTIFF or COG.
-open_cog(url, ...)
-    Read a Cloud Optimized GeoTIFF from an HTTP URL.
+write_vrt(vrt_path, source_files, ...)
+    Generate a VRT mosaic XML from a list of GeoTIFF files.
 """
 from __future__ import annotations
 
@@ -20,7 +20,7 @@
 from ._reader import read_to_array
 from ._writer import write
 
-__all__ = ['read_geotiff', 'write_geotiff', 'write_vrt']
+__all__ = ['open_geotiff', 'to_geotiff', 'write_vrt']
 
 
 def _wkt_to_epsg(wkt_or_proj: str) -> int | None:
@@ -98,7 +98,55 @@ def _coords_to_transform(da: xr.DataArray) -> GeoTransform | None:
     )
 
 
-def read_geotiff(source: str, *, window=None,
+def _read_geo_info(source: str):
+    """Read only the geographic metadata and image dimensions from a GeoTIFF.
+
+    Returns (geo_info, height, width) without reading pixel data.
+    """
+    from ._geotags import extract_geo_info
+    from ._header import parse_all_ifds, parse_header
+
+    with open(source, 'rb') as f:
+        import mmap
+        data = mmap.mmap(f.fileno(), 0, access=mmap.ACCESS_READ)
+    try:
+        header = parse_header(data)
+        ifds = parse_all_ifds(data, header)
+        ifd = ifds[0]
+        geo_info = extract_geo_info(ifd, data, header.byte_order)
+        return geo_info, ifd.height, ifd.width
+    finally:
+        data.close()
+
+
+def _extent_to_window(transform, file_height, file_width,
+                      y_min, y_max, x_min, x_max):
+    """Convert geographic extent to pixel window (row_start, col_start, row_stop, col_stop).
+
+    Clamps to file bounds.
+    """
+    # Pixel coords from geographic coords
+    col_start = (x_min - transform.origin_x) / transform.pixel_width
+    col_stop = (x_max - transform.origin_x) / transform.pixel_width
+
+    row_start = (y_max - transform.origin_y) / transform.pixel_height
+    row_stop = (y_min - transform.origin_y) / transform.pixel_height
+
+    # pixel_height is typically negative, so row_start/row_stop may be swapped
+    if row_start > row_stop:
+        row_start, row_stop = row_stop, row_start
+    if col_start > col_stop:
+        col_start, col_stop = col_stop, col_start
+
+    row_start = max(0, int(np.floor(row_start)))
+    col_start = max(0, int(np.floor(col_start)))
+    row_stop = min(file_height, int(np.ceil(row_stop)))
+    col_stop = min(file_width, int(np.ceil(col_stop)))
+
+    return (row_start, col_start, row_stop, col_stop)
+
+
+def open_geotiff(source: str, *, window=None,
                  overview_level: int | None = None,
                  band: int | None = None,
                  name: str | None = None,
@@ -285,18 +333,18 @@ def _is_gpu_data(data) -> bool:
     return isinstance(data, _cupy_type)
 
 
-def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
-                  crs: int | str | None = None,
-                  nodata=None,
-                  compression: str = 'deflate',
-                  tiled: bool = True,
-                  tile_size: int = 256,
-                  predictor: bool = False,
-                  cog: bool = False,
-                  overview_levels: list[int] | None = None,
-                  overview_resampling: str = 'mean',
-                  bigtiff: bool | None = None,
-                  gpu: bool | None = None) -> None:
+def to_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
+               crs: int | str | None = None,
+               nodata=None,
+               compression: str = 'deflate',
+               tiled: bool = True,
+               tile_size: int = 256,
+               predictor: bool = False,
+               cog: bool = False,
+               overview_levels: list[int] | None = None,
+               overview_resampling: str = 'mean',
+               bigtiff: bool | None = None,
+               gpu: bool | None = None) -> None:
     """Write data as a GeoTIFF or Cloud Optimized GeoTIFF.
 
     Automatically dispatches to GPU compression when:
@@ -442,14 +490,6 @@ def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
     )
 
 
-def open_cog(url: str, **kwargs) -> xr.DataArray:
-    """Deprecated: use ``read_geotiff(url, ...)`` instead.
-
-    read_geotiff handles HTTP URLs, cloud URIs, and local files.
-    """
-    return read_geotiff(url, **kwargs)
-
-
 def read_geotiff_dask(source: str, *, chunks: int | tuple = 512,
                       overview_level: int | None = None,
                       name: str | None = None) -> xr.DataArray:

xrspatial/geotiff/tests/bench_vs_rioxarray.py

@@ -36,12 +36,12 @@ def _fmt_ms(seconds):
 def check_consistency(path):
     """Compare pixel values and geo metadata between the two readers."""
     import rioxarray  # noqa: F401
-    from xrspatial.geotiff import read_geotiff
+    from xrspatial.geotiff import open_geotiff
 
     rio_da = xr.open_dataarray(path, engine='rasterio')
     rio_arr = rio_da.squeeze('band').values.astype(np.float64)
 
-    our_da = read_geotiff(path)
+    our_da = open_geotiff(path)
     our_arr = our_da.values.astype(np.float64)
 
     # Shape
@@ -96,7 +96,7 @@ def check_consistency(path):
 def bench_read(path, runs=10):
     """Benchmark read performance."""
     import rioxarray  # noqa: F401
-    from xrspatial.geotiff import read_geotiff
+    from xrspatial.geotiff import open_geotiff
 
     def rio_read():
         da = xr.open_dataarray(path, engine='rasterio')
@@ -105,7 +105,7 @@ def rio_read():
         return da
 
     def our_read():
-        return read_geotiff(path)
+        return open_geotiff(path)
 
     rio_time, _ = _timer(rio_read, warmup=2, runs=runs)
     our_time, _ = _timer(our_read, warmup=2, runs=runs)
@@ -120,7 +120,7 @@ def our_read():
 def bench_write(shape=(512, 512), compression='deflate', runs=5):
     """Benchmark write performance."""
     import rioxarray  # noqa: F401
-    from xrspatial.geotiff import write_geotiff
+    from xrspatial.geotiff import to_geotiff
     from xrspatial.geotiff._geotags import GeoTransform
 
     rng = np.random.RandomState(42)
@@ -150,7 +150,7 @@ def rio_write():
 
     def our_write():
         p = os.path.join(tmpdir, 'our_out.tif')
-        write_geotiff(da_ours, p, compression=compression, tiled=False)
+        to_geotiff(da_ours, p, compression=compression, tiled=False)
         return os.path.getsize(p)
 
     rio_time, rio_size = _timer(rio_write, warmup=1, runs=runs)
@@ -166,7 +166,7 @@ def our_write():
 def bench_round_trip(shape=(256, 256), compression='deflate'):
     """Write with our module, read back with rioxarray, and vice versa."""
     import rioxarray  # noqa: F401
-    from xrspatial.geotiff import read_geotiff, write_geotiff
+    from xrspatial.geotiff import open_geotiff, to_geotiff
 
     rng = np.random.RandomState(99)
     arr = rng.rand(*shape).astype(np.float32)
@@ -179,7 +179,7 @@ def bench_round_trip(shape=(256, 256), compression='deflate'):
     our_path = os.path.join(tmpdir, 'ours.tif')
     da_ours = xr.DataArray(arr, dims=['y', 'x'], coords={'y': y, 'x': x},
                            attrs={'crs': 4326})
-    write_geotiff(da_ours, our_path, compression=compression, tiled=False)
+    to_geotiff(da_ours, our_path, compression=compression, tiled=False)
 
     rio_da = xr.open_dataarray(our_path, engine='rasterio')
     rio_arr = rio_da.squeeze('band').values if 'band' in rio_da.dims else rio_da.values
@@ -198,7 +198,7 @@ def bench_round_trip(shape=(256, 256), compression='deflate'):
     else:
         da_rio.rio.to_raster(rio_path)
 
-    our_da = read_geotiff(rio_path)
+    our_da = open_geotiff(rio_path)
     our_arr = our_da.values
 
     diff2 = float(np.nanmax(np.abs(arr - our_arr)))

xrspatial/geotiff/tests/test_cog.py

@@ -5,7 +5,7 @@
 import pytest
 import xarray as xr
 
-from xrspatial.geotiff import read_geotiff, write_geotiff
+from xrspatial.geotiff import open_geotiff, to_geotiff
 from xrspatial.geotiff._header import parse_header, parse_all_ifds
 from xrspatial.geotiff._writer import write
 from xrspatial.geotiff._geotags import GeoTransform, extract_geo_info
@@ -83,52 +83,52 @@ def test_read_write_round_trip(self, tmp_path):
         )
 
         path = str(tmp_path / 'round_trip.tif')
-        write_geotiff(da, path, compression='deflate', tiled=False)
+        to_geotiff(da, path, compression='deflate', tiled=False)
 
-        result = read_geotiff(path)
+        result = open_geotiff(path)
         np.testing.assert_array_almost_equal(result.values, data, decimal=5)
         assert result.attrs.get('crs') == 4326
 
-    def test_read_geotiff_name(self, tmp_path):
+    def test_open_geotiff_name(self, tmp_path):
         """DataArray name defaults to filename stem."""
         arr = np.zeros((4, 4), dtype=np.float32)
         path = str(tmp_path / 'myfile.tif')
         write(arr, path, compression='none', tiled=False)
 
-        da = read_geotiff(path)
+        da = open_geotiff(path)
         assert da.name == 'myfile'
 
-    def test_read_geotiff_custom_name(self, tmp_path):
+    def test_open_geotiff_custom_name(self, tmp_path):
         arr = np.zeros((4, 4), dtype=np.float32)
         path = str(tmp_path / 'test.tif')
         write(arr, path, compression='none', tiled=False)
 
-        da = read_geotiff(path, name='custom')
+        da = open_geotiff(path, name='custom')
         assert da.name == 'custom'
 
     def test_write_numpy_array(self, tmp_path):
-        """write_geotiff should accept raw numpy arrays too."""
+        """to_geotiff should accept raw numpy arrays too."""
         arr = np.arange(16, dtype=np.float32).reshape(4, 4)
         path = str(tmp_path / 'numpy.tif')
-        write_geotiff(arr, path, compression='none')
+        to_geotiff(arr, path, compression='none')
 
-        result = read_geotiff(path)
+        result = open_geotiff(path)
         np.testing.assert_array_equal(result.values, arr)
 
     def test_write_3d_rgb(self, tmp_path):
         """3D arrays (height, width, bands) should write multi-band."""
         arr = np.zeros((4, 4, 3), dtype=np.uint8)
         arr[:, :, 0] = 255  # red channel
         path = str(tmp_path / 'rgb.tif')
-        write_geotiff(arr, path, compression='none')
+        to_geotiff(arr, path, compression='none')
 
-        result = read_geotiff(path)
+        result = open_geotiff(path)
         np.testing.assert_array_equal(result.values, arr)
 
     def test_write_rejects_4d(self, tmp_path):
         arr = np.zeros((2, 3, 4, 4), dtype=np.float32)
         with pytest.raises(ValueError, match="Expected 2D or 3D"):
-            write_geotiff(arr, str(tmp_path / 'bad.tif'))
+            to_geotiff(arr, str(tmp_path / 'bad.tif'))
 
 
 def read_to_array_local(path):