Rename GeoTIFF API to xarray conventions (#1047) (#1061)

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.
Includes accessor tests and user guide notebook.

Author: brendancol

README.md

@@ -141,25 +141,29 @@ Native GeoTIFF and Cloud Optimized GeoTIFF reader/writer. No GDAL required.
 
 | Name | Description | NumPy | Dask | CuPy GPU | Dask+CuPy GPU | Cloud |
 |:-----|:------------|:-----:|:----:|:--------:|:-------------:|:-----:|
-| [read_geotiff](xrspatial/geotiff/__init__.py) | Read GeoTIFF / COG / VRT | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
-| [write_geotiff](xrspatial/geotiff/__init__.py) | Write DataArray as GeoTIFF / COG | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
+| [open_geotiff](xrspatial/geotiff/__init__.py) | Read GeoTIFF / COG / VRT | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
+| [to_geotiff](xrspatial/geotiff/__init__.py) | Write DataArray as GeoTIFF / COG | ✅️ | ✅️ | ✅️ | ✅️ | ✅️ |
 | [write_vrt](xrspatial/geotiff/__init__.py) | Generate VRT mosaic from GeoTIFFs | ✅️ | | | | |
 
-`read_geotiff` and `write_geotiff` auto-dispatch to the correct backend:
+`open_geotiff` and `to_geotiff` auto-dispatch to the correct backend:
 
 ```python
-read_geotiff('dem.tif')                              # NumPy
-read_geotiff('dem.tif', chunks=512)                  # Dask
-read_geotiff('dem.tif', gpu=True)                    # CuPy (nvCOMP + GDS)
-read_geotiff('dem.tif', gpu=True, chunks=512)        # Dask + CuPy
-read_geotiff('https://example.com/cog.tif')          # HTTP COG
-read_geotiff('s3://bucket/dem.tif')                  # Cloud (S3/GCS/Azure)
-read_geotiff('mosaic.vrt')                           # VRT mosaic (auto-detected)
-
-write_geotiff(cupy_array, 'out.tif')                 # auto-detects GPU
-write_geotiff(data, 'out.tif', gpu=True)             # force GPU compress
-write_geotiff(data, 'ortho.tif', compression='jpeg') # JPEG for orthophotos
+open_geotiff('dem.tif')                              # NumPy
+open_geotiff('dem.tif', chunks=512)                  # Dask
+open_geotiff('dem.tif', gpu=True)                    # CuPy (nvCOMP + GDS)
+open_geotiff('dem.tif', gpu=True, chunks=512)        # Dask + CuPy
+open_geotiff('https://example.com/cog.tif')          # HTTP COG
+open_geotiff('s3://bucket/dem.tif')                  # Cloud (S3/GCS/Azure)
+open_geotiff('mosaic.vrt')                           # VRT mosaic (auto-detected)
+
+to_geotiff(cupy_array, 'out.tif')                    # auto-detects GPU
+to_geotiff(data, 'out.tif', gpu=True)                # force GPU compress
+to_geotiff(data, 'ortho.tif', compression='jpeg')    # JPEG for orthophotos
 write_vrt('mosaic.vrt', ['tile1.tif', 'tile2.tif'])  # generate VRT
+
+# Accessor methods
+da.xrs.to_geotiff('out.tif', compression='lzw')     # write from DataArray
+ds.xrs.open_geotiff('large_dem.tif')                 # read windowed to Dataset extent
 ```
 
 **Compression codecs:** Deflate, LZW (Numba JIT), ZSTD, PackBits, JPEG (Pillow), uncompressed
@@ -493,10 +497,10 @@ Importing `xrspatial` registers an `.xrs` accessor on DataArrays and Datasets, g
 
 ```python
 import xrspatial
-from xrspatial.geotiff import read_geotiff
+from xrspatial.geotiff import open_geotiff
 
 # Read a GeoTIFF (no GDAL required)
-elevation = read_geotiff('dem.tif')
+elevation = open_geotiff('dem.tif')
 
 # Surface analysis — call operations directly on the DataArray
 slope = elevation.xrs.slope()

docs/source/user_guide/multispectral.ipynb

@@ -41,7 +41,7 @@
    },
    "outputs": [],
    "source": [
-    "import datashader as ds\nfrom datashader.colors import Elevation\nimport datashader.transfer_functions as tf\nfrom datashader.transfer_functions import shade\nfrom datashader.transfer_functions import stack\nfrom datashader.transfer_functions import dynspread\nfrom datashader.transfer_functions import set_background\nfrom datashader.transfer_functions import Images, Image\nfrom datashader.utils import orient_array\nimport numpy as np\nimport xarray as xr\nfrom xrspatial.geotiff import read_geotiff"
+    "import datashader as ds\nfrom datashader.colors import Elevation\nimport datashader.transfer_functions as tf\nfrom datashader.transfer_functions import shade\nfrom datashader.transfer_functions import stack\nfrom datashader.transfer_functions import dynspread\nfrom datashader.transfer_functions import set_background\nfrom datashader.transfer_functions import Images, Image\nfrom datashader.utils import orient_array\nimport numpy as np\nimport xarray as xr\nfrom xrspatial.geotiff import open_geotiff"
    ]
   },
   {
@@ -132,7 +132,7 @@
     }
    ],
    "source": [
-    "SCENE_ID = \"LC80030172015001LGN00\"\nEXTS = {\n    \"blue\": \"B2\",\n    \"green\": \"B3\",\n    \"red\": \"B4\",\n    \"nir\": \"B5\",\n}\n\ncvs = ds.Canvas(plot_width=1024, plot_height=1024)\nlayers = {}\nfor name, ext in EXTS.items():\n    layer = read_geotiff(f\"../../../xrspatial-examples/data/{SCENE_ID}_{ext}.tiff\", band=0)\n    layer.name = name\n    layer = cvs.raster(layer, agg=\"mean\")\n    layer.data = orient_array(layer)\n    layers[name] = layer\nlayers"
+    "SCENE_ID = \"LC80030172015001LGN00\"\nEXTS = {\n    \"blue\": \"B2\",\n    \"green\": \"B3\",\n    \"red\": \"B4\",\n    \"nir\": \"B5\",\n}\n\ncvs = ds.Canvas(plot_width=1024, plot_height=1024)\nlayers = {}\nfor name, ext in EXTS.items():\n    layer = open_geotiff(f\"../../../xrspatial-examples/data/{SCENE_ID}_{ext}.tiff\", band=0)\n    layer.name = name\n    layer = cvs.raster(layer, agg=\"mean\")\n    layer.data = orient_array(layer)\n    layers[name] = layer\nlayers"
    ]
   },
   {

examples/user_guide/25_GLCM_Texture.ipynb

@@ -282,7 +282,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import os\nfrom xrspatial.geotiff import read_geotiff\n\n\nCOG_URL = (\n    'https://sentinel-cogs.s3.us-west-2.amazonaws.com/'\n    'sentinel-s2-l2a-cogs/10/S/EG/2023/9/'\n    'S2B_10SEG_20230921_0_L2A/B08.tif'\n)\n\ntry:\n    nir_da = read_geotiff(COG_URL, band=0, window=(2100, 5300, 2600, 5800))\n    nir = nir_da.values.astype(np.float64)\n    print(f'Downloaded NIR band: {nir.shape}, range {nir.min():.0f} to {nir.max():.0f}')\nexcept Exception as exc:\n    print(f'Remote read failed ({exc}), using synthetic fallback')\n    rng_sat = np.random.default_rng(99)\n    nir = np.zeros((500, 500), dtype=np.float64)\n    nir[:, 250:] = rng_sat.normal(80, 10, (500, 250)).clip(20, 200)\n    nir[:, :250] = rng_sat.normal(1800, 400, (500, 250)).clip(300, 4000)\n\nsatellite = xr.DataArray(nir, dims=['y', 'x'],\n                         coords={'y': np.arange(nir.shape[0], dtype=float),\n                                 'x': np.arange(nir.shape[1], dtype=float)})\n\nfig, ax = plt.subplots(figsize=(7, 7))\nsatellite.plot.imshow(ax=ax, cmap='gray', vmax=float(np.percentile(nir, 98)),\n                      add_colorbar=False)\nax.set_title('Sentinel-2 NIR band')\nax.set_axis_off()\nplt.tight_layout()"
+    "import os\nfrom xrspatial.geotiff import open_geotiff\n\n\nCOG_URL = (\n    'https://sentinel-cogs.s3.us-west-2.amazonaws.com/'\n    'sentinel-s2-l2a-cogs/10/S/EG/2023/9/'\n    'S2B_10SEG_20230921_0_L2A/B08.tif'\n)\n\ntry:\n    nir_da = open_geotiff(COG_URL, band=0, window=(2100, 5300, 2600, 5800))\n    nir = nir_da.values.astype(np.float64)\n    print(f'Downloaded NIR band: {nir.shape}, range {nir.min():.0f} to {nir.max():.0f}')\nexcept Exception as exc:\n    print(f'Remote read failed ({exc}), using synthetic fallback')\n    rng_sat = np.random.default_rng(99)\n    nir = np.zeros((500, 500), dtype=np.float64)\n    nir[:, 250:] = rng_sat.normal(80, 10, (500, 250)).clip(20, 200)\n    nir[:, :250] = rng_sat.normal(1800, 400, (500, 250)).clip(300, 4000)\n\nsatellite = xr.DataArray(nir, dims=['y', 'x'],\n                         coords={'y': np.arange(nir.shape[0], dtype=float),\n                                 'x': np.arange(nir.shape[1], dtype=float)})\n\nfig, ax = plt.subplots(figsize=(7, 7))\nsatellite.plot.imshow(ax=ax, cmap='gray', vmax=float(np.percentile(nir, 98)),\n                      add_colorbar=False)\nax.set_title('Sentinel-2 NIR band')\nax.set_axis_off()\nplt.tight_layout()"
    ]
   },
   {

examples/user_guide/35_GeoTIFF_IO.ipynb

@@ -35,7 +35,7 @@
    },
    "outputs": [],
    "source": [
-    "import pandas\nimport matplotlib.pyplot as plt\nimport geopandas as gpd\n\nimport xarray as xr\nimport numpy as np\nimport cupy\nfrom xrspatial.geotiff import read_geotiff\n\nimport xrspatial"
+    "import pandas\nimport matplotlib.pyplot as plt\nimport geopandas as gpd\n\nimport xarray as xr\nimport numpy as np\nimport cupy\nfrom xrspatial.geotiff import open_geotiff\n\nimport xrspatial"
    ]
   },
   {
@@ -66,7 +66,7 @@
    },
    "outputs": [],
    "source": [
-    "file_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\n\nraster = read_geotiff(file_name, band=0)\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
+    "file_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\n\nraster = open_geotiff(file_name, band=0)\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
    ]
   },
   {

examples/user_guide/images/geotiff_io_preview.png

@@ -47,7 +47,7 @@
    },
    "outputs": [],
    "source": [
-    "import xarray as xr\nfrom xrspatial.geotiff import read_geotiff\nimport xrspatial\n\nfile_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\nraster = read_geotiff(file_name, band=0)\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
+    "import xarray as xr\nfrom xrspatial.geotiff import open_geotiff\nimport xrspatial\n\nfile_name = '../xrspatial-examples/data/colorado_merge_3arc_resamp.tif'\nraster = open_geotiff(file_name, band=0)\nraster.name = 'Colorado Elevation Raster'\n\nxmin, xmax = raster.x.data.min(), raster.x.data.max()\nymin, ymax = raster.y.data.min(), raster.y.data.max()\n\nxmin, xmax, ymin, ymax"
    ]
   },
   {

examples/viewshed_gpu.ipynb

@@ -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
@@ -482,6 +482,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:
@@ -820,3 +832,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)

examples/xarray-spatial_classification-methods.ipynb

@@ -4,12 +4,10 @@
 
 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.
 """
 from __future__ import annotations
 
@@ -20,7 +18,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 +96,53 @@ 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.
+    """
+    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
+
+    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,7 +329,7 @@ def _is_gpu_data(data) -> bool:
     return isinstance(data, _cupy_type)
 
 
-def write_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
+def to_geotiff(data: xr.DataArray | np.ndarray, path: str, *,
                   crs: int | str | None = None,
                   nodata=None,
                   compression: str = 'deflate',
@@ -445,14 +489,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: