GeoTensor

The GeoTensor class is a subclass of np.ndarray that stores the spatial affine transform, the coordinate reference system (crs) and no data value (called fill_value_default). When a GeoTensor is sliced, its transform attribute is shifted accordingly. Additionally its transform is also shifted if the GeoTensor is padded. GeoTensors are restricted to be 2D, 3D or 4D arrays and their two last dimensions are assumed to be the y and x spatial axis.

GeoTensor implements the GeoData protocol. This makes it fully compatible with all read_* methods in the library, allowing for operations like reprojection, subseting from spatial coordinates, mosaicking or vectorization.

For a detailed guide on working with the NumPy API in GeoTensors, see the GeoTensor NumPy API tutorial which covers slicing, operations, masking, and more practical examples.

As a subclass of np.ndarray, operations with GeoTensor objects work similar than operations with np.ndarrays. However, there are some restrictions that we have implemented to keep consistency with the GeoTensor concept. If you need to use the numpy implementation you can access the bare numpy object with the .values attribute. Below there's a list with restrictions on numpy operations:

1. Slicing a GeoTensor is more restrictive than a numpy array. It only allows to slice with lists, numbers or slices. In particular the spatial dimensions can only be sliced with slices. Slicing for inplace modification is not restricted (i.e. you can slice with boolean arrays to modify certain values of the object). See isel and getitem methods.

2. Binary operations (such as add add, multiply mul, ==, | etc) check, for GeoTensor inputs, if they have the same_extent; that is, same transform crs and spatial dimensions (width and height).

3. squeeze, expand_dims and transpose make sure spatial dimensions (last two axes) are not modified and kept at the end of the array.

4. concatenate and stack make sure all operated GeoTensors have same_extent and shape. concatenate does not allow to concatenate on the spatial dims.

5. Reductions (such as np.mean or np.all) return GeoTensor object if the spatial dimensions are preserved and np.ndarray or scalars otherwise. This is handled by the array_ufunc method.

Additional Features

• Masking utilities: Methods like validmask() and invalidmask() create boolean masks based on the fill_value_default.
• Window-based access: read_from_window() and write_from_window() provide window-based operations using rasterio's window system.
• File I/O: Class methods load_file() and load_bytes() for easily loading GeoTensors from files or memory.
• Footprint extraction: Methods to extract the valid data footprint as vector geometries (footprint() and valid_footprint()).
• Coordinate transformation: The meshgrid() method creates coordinate arrays for the spatial dimensions.
• Resizing: The resize() method changes the spatial resolution while maintaining geospatial information correct (i.e. changing the spatial resolution of the transform).
• Metadata storage: GeoTensor includes an attrs dictionary for storing additional metadata like tags and band descriptions.

API Reference

GeoTensor: A NumPy ndarray subclass with geospatial metadata.

This module provides the core data structure for georeader - a numpy array subclass that carries geospatial metadata (CRS, transform, fill value) through all operations. GeoTensor enables numpy-style array operations while preserving georeferencing.

Memory Layout & Dimensions

GeoTensor supports 2D, 3D, and 4D arrays. The last two dimensions are always spatial::

┌─────────────────────────────────────────────────────────────────────────┐
│                     GEOTENSOR DIMENSION CONVENTIONS                     │
├─────────────────────────────────────────────────────────────────────────┤
│                                                                         │
│  2D: (H, W)           Single-band raster (e.g., DEM, mask)              │
│                       shape = (height, width)                           │
│                                                                         │
│  3D: (C, H, W)        Multi-band raster (e.g., RGB, multispectral)      │
│                       shape = (channels, height, width)                 │
│                                                                         │
│  4D: (T, C, H, W)     Time-series cube (e.g., satellite time stack)     │
│                       shape = (time, channels, height, width)           │
│                                                                         │
│  Dimension names:  dims = ("y", "x") or ("band", "y", "x") or           │
│                          ("time", "band", "y", "x")                     │
│                                                                         │
│  Note: "y" decreases downward (row index), "x" increases rightward      │
└─────────────────────────────────────────────────────────────────────────┘

Coordinate Systems & Transform

GeoTensor uses an affine transform to map pixel coordinates to geographic coordinates::

┌─────────────────────────────────────────────────────────────────────────┐
│              AFFINE TRANSFORM: PIXEL ↔ GEOGRAPHIC COORDINATES           │
├─────────────────────────────────────────────────────────────────────────┤
│                                                                         │
│  Pixel Space (row, col)              Geographic Space (x, y)            │
│  ┌───┬───┬───┬───┬───┐              ┌─────────────────────────┐         │
│  │0,0│0,1│0,2│0,3│...│              │OriginX ──────────────►  │         │
│  ├───┼───┼───┼───┼───┤     ═══►     │OriginY                  │         │
│  │1,0│1,1│   │   │   │   Transform  │   │                     │         │
│  ├───┼───┼───┼───┼───┤              │   ▼                     │         │
│  │2,0│   │   │   │   │              │        (CRS units)      │         │
│  └───┴───┴───┴───┴───┘              └─────────────────────────┘         │
│                                                                         │
│  Affine Transform:  | a  b  c |    x_geo = a * col + b * row + c        │
│                     | d  e  f |    y_geo = d * col + e * row + f        │
│                                                                         │
│  Typical (north-up):  a = pixel_width (positive)                        │
│                       e = -pixel_height (negative, y decreases down)    │
│                       c = origin_x (upper-left corner)                  │
│                       f = origin_y (upper-left corner)                  │
│                       b, d = 0 (no rotation/shear)                      │
│                                                                         │
│  Resolution:  res = (|a|, |e|) in CRS units (e.g., meters, degrees)     │
└─────────────────────────────────────────────────────────────────────────┘

NumPy Subclass Behavior

GeoTensor inherits from np.ndarray, so all numpy operations work natively::

┌─────────────────────────────────────────────────────────────────────────┐
│                    NUMPY OPERATIONS PRESERVE GEOSPATIAL INFO            │
├─────────────────────────────────────────────────────────────────────────┤
│                                                                         │
│  Arithmetic:     gt + 1, gt * 2, gt1 + gt2 (same extent required)       │
│  Comparison:     gt > 0, gt == other                                    │
│  Ufuncs:         np.sqrt(gt), np.log(gt), np.clip(gt, 0, 1)             │
│  Aggregation:    gt.mean(), gt.sum(axis=0)  # Returns scalar/array      │
│  Slicing:        gt[0], gt[:, 10:20, 10:20]  # Updates transform        │
│                                                                         │
│  Important: Operations that change spatial dimensions (slicing)         │
│  automatically update the transform to maintain georeferencing!         │
│                                                                         │
│  Example:                                                               │
│    gt_subset = gt[:, 100:200, 100:200]                                  │
│    # gt_subset.transform origin shifted by (100*res_x, 100*res_y)       │
└─────────────────────────────────────────────────────────────────────────┘

Key Attributes

=============== ============================================================ Attribute Description =============== ============================================================ values View of underlying numpy array data transform rasterio.Affine mapping pixels to geographic coordinates crs Coordinate reference system (EPSG code or WKT) fill_value_default Value used for out-of-bounds reads (default: 0) bounds Geographic bounds (minx, miny, maxx, maxy) in CRS units res Pixel resolution as (x_res, y_res) tuple shape Array shape (always ends with height, width) dims Dimension names tuple ("band", "y", "x") etc. attrs Optional metadata dictionary =============== ============================================================

Quick Start

Creating a GeoTensor::

import numpy as np
import rasterio
from georeader.geotensor import GeoTensor

# Create from scratch
data = np.random.rand(3, 100, 100)  # 3 bands, 100x100 pixels
transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)  # 10m resolution
gt = GeoTensor(data, transform=transform, crs="EPSG:32630")

# Access properties
print(gt.bounds)      # (500000.0, 4649000.0, 501000.0, 4650000.0)
print(gt.res)         # (10.0, 10.0)
print(gt.shape)       # (3, 100, 100)

Array operations::

# Arithmetic (preserves geospatial info)
gt_scaled = gt * 10000
gt_offset = gt + 0.1

# Band math (NDVI example)
red, nir = gt[0], gt[1]
ndvi = (nir - red) / (nir + red + 1e-10)

# Spatial slicing (automatically updates transform)
subset = gt[:, 50:100, 50:100]  # 50x50 subset
print(subset.bounds)  # Different from gt.bounds!

Module Functions

=================== ============================================================ Function Description =================== ============================================================ GeoTensor Core class - numpy array with geospatial metadata _vsi_path Convert cloud URLs to GDAL VSI paths get_rio_options_path Configure rasterio environment options =================== ============================================================

See Also

georeader.read : Functions to subset GeoTensors from polygons or windows. georeader.save : Functions to save GeoTensors to GeoTIFF files. georeader.rasterio_reader : RasterioReader class for lazy file access. georeader.abstract_reader : GeoData protocol that GeoTensor implements.

References

• Rasterio Affine transforms: https://rasterio.readthedocs.io/en/latest/topics/transforms.html
• NumPy subclassing: https://numpy.org/doc/stable/user/basics.subclassing.html

GeoTensor

Bases: ndarray

A numpy ndarray subclass with geospatial metadata.

GeoTensor wraps numpy arrays with geospatial information including coordinate reference system (CRS), affine transform, and fill value. It supports 2D, 3D, or 4D tensors where the last two dimensions are always spatial (y, x).

Use :meth:__new__ to create instances with values, transform, and CRS.

Attributes:

Name	Type	Description
values	NDArray	numpy or torch tensor
transform	Affine	affine geospatial transform
crs	Any	coordinate reference system
fill_value_default	Optional[Union[int, float]]	Value to fill when reading out of bounds. Could be None. Defaults to 0.
shape	Tuple	shape of the tensor
res	Tuple[float, float]	resolution of the tensor
dtype		data type of the tensor
height	int	height of the tensor
width	int	width of the tensor
count	int	number of bands in the tensor
bounds	Tuple[float, float, float, float]	bounds of the tensor
dims	Tuple[str]	names of the dimensions
attrs	Dict[str, Any]	dictionary with the attributes of the GeoTensor

Examples:

>>> import numpy as np
>>> transform = rasterio.Affine(1, 0, 0, 0, -1, 0)
>>> crs = "EPSG:4326"
>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)

Source code in georeader/geotensor.py

class GeoTensor(np.ndarray):
    """
    A numpy ndarray subclass with geospatial metadata.

    GeoTensor wraps numpy arrays with geospatial information including coordinate
    reference system (CRS), affine transform, and fill value. It supports 2D, 3D,
    or 4D tensors where the last two dimensions are always spatial (y, x).

    Use :meth:`__new__` to create instances with values, transform, and CRS.

    Attributes:
        values (NDArray): numpy or torch tensor
        transform (rasterio.Affine): affine geospatial transform
        crs (Any): coordinate reference system
        fill_value_default (Optional[Union[int, float]], optional): Value to fill when
            reading out of bounds. Could be None. Defaults to 0.
        shape (Tuple): shape of the tensor
        res (Tuple[float, float]): resolution of the tensor
        dtype: data type of the tensor
        height (int): height of the tensor
        width (int): width of the tensor
        count (int): number of bands in the tensor
        bounds (Tuple[float, float, float, float]): bounds of the tensor
        dims (Tuple[str]): names of the dimensions
        attrs (Dict[str, Any]): dictionary with the attributes of the GeoTensor

    Examples:
        >>> import numpy as np
        >>> transform = rasterio.Affine(1, 0, 0, 0, -1, 0)
        >>> crs = "EPSG:4326"
        >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)

    """

    def __new__(
        cls,
        values: NDArray,
        transform: rasterio.Affine,
        crs: Any,
        fill_value_default: Optional[Union[int, float]] = 0,
        attrs: Optional[Dict[str, Any]] = None,
    ):
        """
        This class is a wrapper around a numpy or torch tensor with geospatial information.

        Args:
            values (NDArray): numpy or torch tensor
            transform (rasterio.Affine): affine geospatial transform
            crs (Any): coordinate reference system
            fill_value_default (Optional[Union[int, float]], optional): Value to fill when
                reading out of bounds. Could be None. Defaults to 0.
            attrs (Optional[Dict[str, Any]], optional): dictionary with the attributes of the GeoTensor
                Defaults to None.

        Raises:
            ValueError: when the shape of the tensor is not 2d, 3d or 4d.
        """
        obj = np.asarray(values).view(cls)

        obj.transform = transform
        obj.crs = crs
        obj.fill_value_default = fill_value_default
        shape = obj.shape
        if (len(shape) < 2) or (len(shape) > 4):
            raise ValueError(f"Expected 2d-4d array found {shape}")

        obj.attrs = attrs if attrs is not None else {}

        return obj

    def __array_finalize__(self, obj: Optional[Union[np.ndarray, Self]]) -> None:
        """
        Initialize attributes when a new GeoTensor is created from an existing array.

        This method is called whenever a new array object is created from an existing array
        (e.g., through slicing, view casting, or copy operations).

        Args:
            obj (Optional[np.ndarray]): The array object from which the new array is created.
                                       Can be None if the array is being created from scratch.
        """
        if obj is None:
            return

        if hasattr(obj, "transform"):
            self.transform: rasterio.Affine = getattr(obj, "transform", None)
        if hasattr(obj, "crs"):
            self.crs = getattr(obj, "crs", None)
        if hasattr(obj, "fill_value_default"):
            self.fill_value_default = getattr(obj, "fill_value_default", None)
        if hasattr(obj, "attrs"):
            self.attrs = getattr(obj, "attrs", None)

    def __array_ufunc__(self, ufunc, method, *inputs, **kwargs) -> Union["GeoTensor", NDArray]:
        """
        Handle NumPy universal functions applied to this GeoTensor.

        This method is called when a NumPy universal function (ufunc) is applied to the GeoTensor.
        It converts GeoTensor inputs to NumPy arrays, applies the ufunc, and converts array results
        back to GeoTensor objects with the same geospatial metadata.

        Args:
            ufunc (np.ufunc): The NumPy universal function being applied.
            method (str): The method of the ufunc ('__call__', 'reduce', etc.).
            inputs (tuple): The input arrays to the ufunc.
            kwargs (dict): Additional keyword arguments to the ufunc.

        Returns:
            Union[GeoTensor, NDArray]: If the result is an array, returns a new GeoTensor with the same
                geospatial attributes. Otherwise, returns the original result.
        """
        # Normal processing for most operations
        inputs_arr = tuple(
            x.view(np.ndarray) if isinstance(x, GeoTensor) else x for x in inputs
        )
        # Handle 'out' argument if present
        out = kwargs.pop("out", None)
        out_arrays = None

        if out:
            # Convert GeoTensor outputs to regular arrays
            if isinstance(out, tuple):
                out_arrays = tuple(
                    o.view(np.ndarray) if isinstance(o, GeoTensor) else o for o in out
                )
            else:
                out_arrays = (
                    (out.view(np.ndarray),) if isinstance(out, GeoTensor) else (out,)
                )
            kwargs["out"] = out_arrays

        # Delegate to numpy's implementation
        result = super().__array_ufunc__(ufunc, method, *inputs_arr, **kwargs)

        cast_to_geotensor = self._preserved_spatial(method, **kwargs)

        # Propagate metadata to output arrays
        if out_arrays:
            for o_orig, o_new in zip(
                out if isinstance(out, tuple) else [out], out_arrays
            ):
                if (
                    cast_to_geotensor
                    and isinstance(o_orig, GeoTensor)
                    and isinstance(o_new, np.ndarray)
                ):
                    o_new = o_orig.array_as_geotensor(o_new)

        # Normal ufunc processing for other cases
        if cast_to_geotensor:
            return self.array_as_geotensor(result)

        return result

    def _preserved_spatial(self, method: str, **kwargs) -> bool:
        """Special handling for reduction operations (sum, mean, max, etc.)"""
        # Extract reduction axis (default is flattening)
        if method != "reduce":
            return True  # No reduction, preserve spatial dims

        axis = kwargs.get("axis", None)

        # Check if reduction preserves spatial structure (last 2 dims untouched)
        if axis is not None:
            if isinstance(axis, int):
                preserve_spatial = axis not in [-1, -2, self.ndim - 1, self.ndim - 2]
            else:  # tuple of axes
                preserve_spatial = all(
                    ax not in [-1, -2, self.ndim - 1, self.ndim - 2] for ax in axis
                )
        else:
            preserve_spatial = False  # Full reduction eliminates all dims

        # For full reductions or spatial dim reductions, return plain array/scalar
        return preserve_spatial

    def array_as_geotensor(
        self,
        result: Union[np.ndarray, Self],
        fill_value_default: Optional[numbers.Number] = None,
    ) -> Self:
        """
        Convert a NumPy array result back to a GeoTensor.

        Args:
            result (Union[np.ndarray, Self]): Any NumPy array or GeoTensor.
            fill_value_default: fill value for the returned GeoTensor.

        Returns:
            Self: A new GeoTensor with the same geospatial attributes as the original.
        """

        # Propagate metadata for array results
        if isinstance(result, np.ndarray):
            if result.shape[-2:] != self.shape[-2:]:
                raise ValueError("Operation altered spatial dimensions!")

            if fill_value_default is None:
                fill_value_default = self.fill_value_default

            result = GeoTensor(
                result,
                transform=self.transform,
                crs=self.crs,
                fill_value_default=fill_value_default,
                attrs=self.attrs,
            )

        return result

    def __array__(self, dtype: Optional[np.dtype] = None) -> np.ndarray:
        """
        Convert the GeoTensor to a standard NumPy array.

        This method is called by np.asarray() and most NumPy functions to get
        the underlying NumPy array representation of this object.

        Args:
            dtype (Optional[np.dtype]): The desired data type for the returned array.
                                       If None, the array's current dtype is preserved.

        Returns:
            np.ndarray: A NumPy array view of this GeoTensor.
        """
        return np.asarray(self.view(np.ndarray), dtype=dtype)

    @property
    def values(self):
        """Return a view of the array (memory shared with original)"""
        return self.view(np.ndarray)

    @property
    def dims(self) -> Tuple[str]:
        # TODO allow different ordering of dimensions?
        shape = self.shape
        if len(shape) == 2:
            dims = ("y", "x")
        elif len(shape) == 3:
            dims = ("band", "y", "x")
        elif len(shape) == 4:
            dims = ("time", "band", "y", "x")
        else:
            raise ValueError(f"Unexpected 2d-4d array found {shape}")

        return dims

    def to_json(self) -> Dict[str, Any]:
        return {
            "values": self.values.tolist(),
            "transform": [
                self.transform.a,
                self.transform.b,
                self.transform.c,
                self.transform.d,
                self.transform.e,
                self.transform.f,
            ],
            "crs": str(self.crs),
            "fill_value_default": self.fill_value_default,
        }

    @classmethod
    def from_json(cls, json: Dict[str, Any]) -> Self:
        return cls(
            np.array(json["values"]),
            rasterio.Affine(*json["transform"]),
            json["crs"],
            json["fill_value_default"],
        )

    @property
    def res(self) -> Tuple[float, float]:
        return window_utils.res(self.transform)

    @property
    def dtype(self):
        return self.values.dtype

    @property
    def height(self) -> int:
        return self.shape[-2]

    @property
    def width(self) -> int:
        return self.shape[-1]

    @property
    def count(self) -> int:
        return self.shape[-3]

    @property
    def bounds(self) -> Tuple[float, float, float, float]:
        return window_bounds(
            rasterio.windows.Window(
                row_off=0, col_off=0, height=self.height, width=self.width
            ),
            self.transform,
        )

    def set_dtype(self, dtype: np.dtype) -> None:
        """Set the dtype of the GeoTensor in-place.

        .. deprecated::
            Use astype() instead for a cleaner numpy-compatible interface.
            This method modifies the array in place which can be confusing.

        Args:
            dtype (np.dtype): Target data type for the array.

        Raises:
            RuntimeError: If the dtype cannot be changed in-place (due to numpy subclass limitations).
        """
        warnings.warn(
            "set_dtype() is deprecated. Use astype() for a cleaner interface.",
            DeprecationWarning,
            stacklevel=2
        )
        original_dtype = self.dtype
        requested_dtype = np.dtype(dtype)

        # If requesting the same dtype, nothing to do
        if original_dtype == requested_dtype:
            return

        # Convert using numpy's astype
        result = np.asarray(self).astype(dtype)
        # Since we can't truly change dtype in-place when sizes differ,
        # we update the internal data by using the ndarray resize mechanism
        # This modifies the underlying buffer
        np.ndarray.resize(self, result.shape, refcheck=False)
        self[:] = result

        # Check if dtype actually changed when a change was requested
        if self.dtype != requested_dtype:
            raise RuntimeError(
                f"set_dtype() failed to change dtype from {original_dtype} to {requested_dtype}. "
                f"This is due to numpy subclass limitations. Use astype() instead to create "
                f"a new GeoTensor with the desired dtype."
            )

    # Not needed due to ufunc implementation?
    # def astype(self, dtype) -> Self:
    #     return GeoTensor(
    #         self.values.astype(dtype), self.transform, self.crs, self.fill_value_default
    #     )

    def meshgrid(self, dst_crs: Optional[Any] = None) -> Tuple[NDArray, NDArray]:
        """
        Create a meshgrid of spatial dimensions of the GeoTensor.

        Args:
            dst_crs (Optional[Any], optional): output coordinate reference system. Defaults to None.

        Returns:
            Tuple[NDArray, NDArray]: 2D arrays of xs and ys coordinates.
        """
        from georeader import griddata

        return griddata.meshgrid(
            self.transform,
            self.width,
            self.height,
            source_crs=self.crs,
            dst_crs=dst_crs,
        )

    def load(self) -> Self:
        return self

    def __copy__(self) -> Self:
        return GeoTensor(
            self.values.copy(), self.transform, self.crs, self.fill_value_default
        )

    def copy(self) -> Self:
        return self.__copy__()

    def same_extent(self, other: Self, precision: float = 1e-3) -> bool:
        """
        Check if two GeoTensors have the same georeferencing (crs, transform and spatial dimensions).

        Args:
            other (GeoTensor | GeoData): GeoTensor to compare with. Other GeoData object can be passed (it requires crs, transform and shape attributes)
            precision (float, optional): precision to compare the transform. Defaults to 1e-3.

        Returns:
            bool: True if both GeoTensors have the same georeferencing.
        """
        return (
            self.transform.almost_equals(other.transform, precision=precision)
            and window_utils.compare_crs(self.crs, other.crs)
            and (self.shape[-2:] == other.shape[-2:])
        )

    def __add__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Add a value or array to this GeoTensor element-wise.

        Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
        When adding two GeoTensors, they must have the same spatial extent
        (matching transform, CRS, and spatial dimensions).

        Broadcasting Rules:
        - Scalar: Added to every pixel
        - Array: Must be broadcastable to self.values shape
        - GeoTensor: Must have identical georeferencing (same_extent)

        Args:
            other (Union[numbers.Number, np.ndarray, GeoTensor]): Value to add. Can be:
                - Scalar (int, float): Added to all pixels
                - np.ndarray: Must be broadcastable with self.values
                - GeoTensor: Must have same spatial extent

        Returns:
            GeoTensor: New GeoTensor with same transform, CRS, and fill_value_default.
                Shape matches self.shape (or broadcast result for arrays).

        Raises:
            ValueError: If other is a GeoTensor and georeferencing doesn't match.

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>>
            >>> # Create sample GeoTensor (3 bands, 100x100 pixels)
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> data = np.random.rand(3, 100, 100)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
            >>>
            >>> # Add scalar to all pixels
            >>> gt_offset = gt + 0.1
            >>> print(gt_offset.shape)  # (3, 100, 100)
            >>>
            >>> # Add per-band offset using broadcasting
            >>> band_offsets = np.array([0.1, 0.2, 0.3])[:, None, None]  # Shape: (3, 1, 1)
            >>> gt_adjusted = gt + band_offsets
            >>>
            >>> # Add two GeoTensors (must have same extent)
            >>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
            >>> gt_sum = gt + gt2
            >>>
            >>> # Error case: mismatched georeferencing
            >>> gt_different = GeoTensor(data, rasterio.Affine(20, 0, 0, 0, -20, 0), crs="EPSG:4326")
            >>> # gt + gt_different  # Raises ValueError

        Note:
            - Result inherits transform, CRS, and fill_value_default from self
            - For GeoTensor addition, use `read.read_reproject_like(other, self)`
              to align georeferencing before adding

        See Also:
            - `same_extent`: Check if two GeoTensors have matching georeferencing
            - `read.read_reproject_like`: Reproject one GeoTensor to match another
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for addition. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values + other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def __sub__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Subtract a value or array from this GeoTensor element-wise.

        Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
        When subtracting two GeoTensors, they must have the same spatial extent
        (matching transform, CRS, and spatial dimensions).

        Args:
            other (Union[numbers.Number, np.ndarray, GeoTensor]): Value to subtract. Can be:
                - Scalar (int, float): Subtracted from all pixels
                - np.ndarray: Must be broadcastable with self.values
                - GeoTensor: Must have same spatial extent

        Returns:
            GeoTensor: New GeoTensor with same transform, CRS, and fill_value_default.

        Raises:
            ValueError: If other is a GeoTensor and georeferencing doesn't match.

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>>
            >>> # Create sample GeoTensor
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> data = np.random.rand(3, 100, 100)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
            >>>
            >>> # Subtract scalar (e.g., remove baseline offset)
            >>> gt_adjusted = gt - 0.1
            >>>
            >>> # Change detection: subtract two GeoTensors
            >>> gt_before = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
            >>> gt_after = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
            >>> change = gt_after - gt_before  # Positive = increase, negative = decrease
            >>>
            >>> # Center data by subtracting mean per band
            >>> band_means = gt.values.mean(axis=(1, 2), keepdims=True)
            >>> gt_centered = gt - band_means

        See Also:
            - `same_extent`: Check if two GeoTensors have matching georeferencing
            - `read.read_reproject_like`: Reproject one GeoTensor to match another
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for substraction. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values - other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def __mul__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Multiply this GeoTensor by a value or array element-wise.

        Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
        Common uses include scaling, applying masks, and band math operations.

        Args:
            other (Union[numbers.Number, np.ndarray, GeoTensor]): Value to multiply. Can be:
                - Scalar (int, float): Multiplies all pixels
                - np.ndarray: Must be broadcastable with self.values
                - GeoTensor: Must have same spatial extent

        Returns:
            GeoTensor: New GeoTensor with result of multiplication.

        Raises:
            ValueError: If other is a GeoTensor and georeferencing doesn't match.

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>>
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> data = np.random.rand(3, 100, 100)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
            >>>
            >>> # Scale all values (e.g., unit conversion)
            >>> gt_scaled = gt * 10000  # Reflectance [0-1] to [0-10000]
            >>>
            >>> # Apply per-band gain coefficients
            >>> gains = np.array([1.1, 1.0, 0.95])[:, None, None]  # Shape: (3, 1, 1)
            >>> gt_calibrated = gt * gains
            >>>
            >>> # Apply binary mask (mask out invalid pixels)
            >>> mask = np.random.rand(100, 100) > 0.5  # Boolean mask
            >>> gt_masked = gt * mask  # Broadcasts mask to all bands
            >>>
            >>> # Element-wise product of two rasters
            >>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
            >>> gt_product = gt * gt2

        Note:
            - For masking, consider using `gt[mask] = fill_value` for in-place updates
            - Result inherits georeferencing from self

        See Also:
            - `__truediv__`: Division operation
            - `__setitem__`: In-place value assignment with masks
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for multiplication. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values * other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def __truediv__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Divide this GeoTensor by a value or array element-wise.

        Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
        Common uses include normalization, ratio calculations, and index computation.

        Args:
            other (Union[numbers.Number, np.ndarray, GeoTensor]): Divisor. Can be:
                - Scalar (int, float): Divides all pixels
                - np.ndarray: Must be broadcastable with self.values
                - GeoTensor: Must have same spatial extent

        Returns:
            GeoTensor: New GeoTensor with result of division.

        Raises:
            ValueError: If other is a GeoTensor and georeferencing doesn't match.

        Note:
            Division by zero produces inf or nan (numpy behavior).
            Consider adding small epsilon for numerical stability: ``gt / (other + 1e-10)``

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>>
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> data = np.random.rand(3, 100, 100)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
            >>>
            >>> # Normalize to [0, 1] range
            >>> gt_norm = gt / gt.values.max()
            >>>
            >>> # Per-band normalization
            >>> band_maxes = gt.values.max(axis=(1, 2))[:, None, None] + 1e-10
            >>> gt_normalized = gt / band_maxes
            >>>
            >>> # Compute NDVI-like ratio: (NIR - Red) / (NIR + Red)
            >>> # Assuming band 0 = Red, band 1 = NIR
            >>> red = gt.values[0]
            >>> nir = gt.values[1]
            >>> ndvi_values = (nir - red) / (nir + red + 1e-10)  # Add epsilon
            >>> ndvi = GeoTensor(ndvi_values, gt.transform, gt.crs)
            >>>
            >>> # Ratio of two rasters
            >>> gt2 = GeoTensor(np.random.rand(3, 100, 100) + 0.1, transform, crs="EPSG:32630")
            >>> ratio = gt / gt2

        See Also:
            - `__mul__`: Multiplication operation
            - `clip`: Clip values to valid range after division
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for division. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values / other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def __and__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Perform bitwise AND operation between two GeoTensors. The georeferencing must match.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): GeoTensor or array-like to AND with.

        Raises:
            ValueError: if the georeferencing does not match when other is a GeoTensor.

        Returns:
            GeoTensor: GeoTensor with the result of the bitwise AND operation.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for bitwise AND. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values & other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def __or__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Perform bitwise OR operation between two GeoTensors. The georeferencing must match.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): GeoTensor or array-like to OR with.

        Raises:
            ValueError: if the georeferencing does not match when other is a GeoTensor.

        Returns:
            GeoTensor: GeoTensor with the result of the bitwise OR operation.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for bitwise OR. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values | other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def __eq__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Element-wise equality comparison between GeoTensors or with a scalar/array.
        The georeferencing must match if comparing with another GeoTensor.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

        Raises:
            ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

        Returns:
            GeoTensor: GeoTensor with boolean values indicating equality.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for comparison. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values == other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    def __ne__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Element-wise inequality comparison between GeoTensors or with a scalar/array.
        The georeferencing must match if comparing with another GeoTensor.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

        Raises:
            ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

        Returns:
            GeoTensor: GeoTensor with boolean values indicating inequality.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for comparison. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values != other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    def __lt__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Element-wise less than comparison between GeoTensors or with a scalar/array.
        The georeferencing must match if comparing with another GeoTensor.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

        Raises:
            ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

        Returns:
            GeoTensor: GeoTensor with boolean values indicating less than relationship.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for comparison. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values < other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    def __le__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Element-wise less than or equal comparison between GeoTensors or with a scalar/array.
        The georeferencing must match if comparing with another GeoTensor.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

        Raises:
            ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

        Returns:
            GeoTensor: GeoTensor with boolean values indicating less than or equal relationship.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for comparison. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values <= other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    def __gt__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Element-wise greater than comparison between GeoTensors or with a scalar/array.
        The georeferencing must match if comparing with another GeoTensor.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

        Raises:
            ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

        Returns:
            GeoTensor: GeoTensor with boolean values indicating greater than relationship.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for comparison. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values > other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    def __ge__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
        """
        Element-wise greater than or equal comparison between GeoTensors or with a scalar/array.
        The georeferencing must match if comparing with another GeoTensor.

        Args:
            other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

        Raises:
            ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

        Returns:
            GeoTensor: GeoTensor with boolean values indicating greater than or equal relationship.
        """
        if isinstance(other, GeoTensor):
            if self.same_extent(other):
                other = other.values
            else:
                raise ValueError(
                    "GeoTensor georref must match for comparison. "
                    "Use `read.read_reproject_like(other, self)` to "
                    "to reproject `other` to `self` georreferencing."
                )

        result_values = self.values >= other

        return GeoTensor(
            result_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    def squeeze(self, axis=None) -> Self:
        """
        Remove single-dimensional entries from the shape of the GeoTensor values.
        It does not squeeze the spatial dimensions (last two dimensions).

        Returns:
            GeoTensor: GeoTensor with the squeezed values.
        """
        if axis is None:
            axis = tuple(range(self.values.ndim - 2))
        else:
            if isinstance(axis, int):
                axis = (axis,)
            # Check if spatial dimesions will be squeezed
            if self.width == 1:
                if any(a in (-1, self.values.ndim - 1) for a in axis):
                    raise ValueError(
                        "Cannot squeeze spatial dimensions. "
                        "Use `squeeze(axis=0)` to squeeze the first dimension."
                    )
            elif self.height == 1:
                if any(a in (-2, self.values.ndim - 2) for a in axis):
                    raise ValueError(
                        "Cannot squeeze spatial dimensions. "
                        "Use `squeeze(axis=0)` to squeeze the first dimension."
                    )

        # squeeze all but last two dimensions
        squeezed_values = np.squeeze(self.values, axis=axis)

        return GeoTensor(
            squeezed_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def expand_dims(self, axis: Union[int, tuple]) -> Self:
        """
        Expand the dimensions of the GeoTensor values while preserving the spatial dimensions.

        This method ensures that no dimensions are added after or in between the spatial dimensions
        (which are always the last two dimensions).

        Args:
            axis (Union[int, tuple]): Position or positions where new axes should be inserted.
                Must be less than the number of dimensions minus 2 (to preserve spatial dims).
                Positions are counted from the first dimension.

        Returns:
            GeoTensor: GeoTensor with the expanded values.

        Raises:
            ValueError: If trying to add dimensions at or after the spatial dimensions.

        Examples:
            >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
            >>> # Add a new dimension at axis 0
            >>> gt_expanded = gt.expand_dims(0)
            >>> assert gt_expanded.shape == (1, 3, 100, 100)
        """
        ndim = len(self.shape)

        # Check if axis is valid (not interfering with spatial dimensions)
        if isinstance(axis, int):
            if axis >= ndim - 2 or axis < -ndim:
                raise ValueError(
                    f"Cannot add dimension at or after spatial dimensions. "
                    f"Axis must be < {ndim - 2} or >= {-ndim}, got {axis}"
                )
            # Convert negative axis to positive
            if axis < 0:
                axis = ndim + axis
        else:  # tuple of axes
            for ax in axis:
                if ax >= ndim - 2 or ax < -ndim:
                    raise ValueError(
                        f"Cannot add dimension at or after spatial dimensions. "
                        f"All axes must be < {ndim - 2} or >= {-ndim}, got {ax}"
                    )

        # Use numpy expand_dims to add the new dimensions
        expanded_values = np.expand_dims(self.values, axis=axis)

        return GeoTensor(
            expanded_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def clip(self, a_min: Optional[np.array], a_max: Optional[np.array]) -> Self:
        """
        Clip the GeoTensor values between the GeoTensor min and max values.

        Args:
            a_min (float): Minimum value.
            a_max (float): Maximum value.

        Returns:
            GeoTensor: GeoTensor with the clipped values.
        """
        clipped_values = np.clip(self.values, a_min, a_max)
        return GeoTensor(
            clipped_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def __getitem__(self, key: Any) -> "GeoTensor":
        """
        Get the values of the GeoTensor using the given key.

        Args:
            key (Any): Key to index the GeoTensor (int, slice, tuple of slices, etc.).

        Returns:
            GeoTensor: GeoTensor with the selected values.
        """
        if not isinstance(key, tuple):
            key = (key,)

        sel_dict = {}
        for i, k in enumerate(self.dims):
            if i < len(key):
                if key[i] is None:
                    raise NotImplementedError(
                        f"Adding axis is not permitted to GeoTensors. Use `expand_dims`"
                    )
                elif isinstance(key[i], type(...)):
                    raise NotImplementedError(
                        f"Using elipsis is not permitted with GeoTensors. Use `values` attribute"
                    )
                sel_dict[k] = key[i]
            else:
                sel_dict[k] = slice(None)

        return self.isel(sel_dict)

    def isel(self, sel: Dict[str, Union[slice, list, int]]) -> Self:
        """
        Select data using dimension-based indexing (xarray-style interface).

        Index into the GeoTensor using named dimensions rather than positional
        indices. Automatically updates the geotransform when slicing spatial
        dimensions (x, y) to maintain correct georeferencing.

        This method provides xarray-compatible syntax while preserving full
        geospatial metadata. Spatial slices shift the transform origin and
        optionally rescale resolution if step > 1.

        Dimension Mapping::

            Standard dimension names for 3D GeoTensor (C, H, W):
            ┌──────────┬────────────────┬────────────────────┐
            │ Dim name │ Array axis      │ Description        │
            ├──────────┼────────────────┼────────────────────┤
            │ "band"   │ axis 0 (C)      │ Spectral bands     │
            │ "y"      │ axis 1 (H)      │ Rows (north-south) │
            │ "x"      │ axis 2 (W)      │ Cols (east-west)   │
            └──────────┴────────────────┴────────────────────┘

        Args:
            sel (Dict[str, Union[slice, list, int]]): Selection dictionary mapping
                dimension names to slice objects, lists of indices, or integers.
                For spatial dimensions ("x", "y"), only slice objects are allowed.

                Common patterns:
                - `{"x": slice(10, 110)}`: Columns 10-109 (100 cols)
                - `{"y": slice(20, 120)}`: Rows 20-119 (100 rows)
                - `{"band": 0}`: Select first band (reduces dimensionality)
                - `{"band": [0, 2, 3]}`: Select bands by list
                - `{"x": slice(0, 100, 2)}`: Every other column (step=2)

        Returns:
            GeoTensor: Sliced GeoTensor with:
                - Updated transform: origin shifted to slice start
                - Updated transform: resolution scaled if step > 1
                - Same CRS and fill_value_default

        Raises:
            NotImplementedError: If dimension name not in self.dims.
            NotImplementedError: If spatial dimension uses non-slice indexing.

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>>
            >>> # Create 3-band 100x100 GeoTensor
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> data = np.random.rand(3, 100, 100).astype(np.float32)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
            >>> print(f"Original: {gt.shape}, origin: ({gt.transform.c}, {gt.transform.f})")
            Original: (3, 100, 100), origin: (500000, 4650000)
            >>>
            >>> # Slice spatial subset (maintains bands)
            >>> subset = gt.isel({"x": slice(20, 80), "y": slice(10, 60)})
            >>> print(f"Subset: {subset.shape}")
            Subset: (3, 50, 60)
            >>> # Transform shifted: new origin at pixel (20, 10)
            >>> print(f"New origin: ({subset.transform.c}, {subset.transform.f})")
            New origin: (500200.0, 4649900.0)
            >>>
            >>> # Downsample with step (every 2nd pixel)
            >>> downsampled = gt.isel({"x": slice(0, 100, 2), "y": slice(0, 100, 2)})
            >>> print(f"Downsampled: {downsampled.shape}")
            Downsampled: (3, 50, 50)
            >>> # Resolution doubled (10m -> 20m)
            >>> print(f"Resolution: {downsampled.res}")
            Resolution: (20.0, 20.0)
            >>>
            >>> # Select single band (reduces to 2D)
            >>> red_band = gt.isel({"band": 0})
            >>> print(f"Red band: {red_band.shape}")
            Red band: (100, 100)
            >>>
            >>> # Combined: subset + band selection
            >>> roi = gt.isel({
            ...     "band": [1, 2],  # Green and blue
            ...     "x": slice(25, 75),
            ...     "y": slice(25, 75)
            ... })
            >>> print(f"ROI: {roi.shape}")
            ROI: (2, 50, 50)

        See Also:
            - `__getitem__`: Positional indexing with tuple syntax: `gt[0, 10:60, 20:80]`
            - `read_from_window`: Slice using rasterio Window objects
            - `pad`: Add padding to spatial dimensions
        """
        for k in sel:
            if k not in self.dims:
                raise NotImplementedError(f"Axis {k} not in {self.dims}")

        slice_tuple = self._slice_tuple(sel)

        # Compùte the window to shift the transform
        slices_window = []
        for k in ["y", "x"]:
            if k in sel:
                if not isinstance(sel[k], slice):
                    raise NotImplementedError(
                        f"Only slice selection supported for x, y dims, found {sel[k]}"
                    )
                slices_window.append(sel[k])
            else:
                size = self.width if (k == "x") else self.height
                slices_window.append(slice(0, size))

        window_current = rasterio.windows.Window.from_slices(
            *slices_window, boundless=False, height=self.height, width=self.width
        )

        transform_current = rasterio.windows.transform(
            window_current, transform=self.transform
        )

        # Scale the spatial transform if the step of the slices > 1
        step_rows = slices_window[0].step
        step_cols = slices_window[1].step
        if step_rows is None:
            step_rows = 1

        if step_cols is None:
            step_cols = 1

        if (step_rows != 1) or (step_cols != 1):
            transform_current = transform_current * Affine.scale(step_cols, step_rows)

        return GeoTensor(
            self.values[slice_tuple],
            transform_current,
            self.crs,
            self.fill_value_default,
            attrs=self.attrs,
        )

    def _slice_tuple(self, sel: Dict[str, Union[slice, list, int]]) -> tuple:
        slice_list = []
        for k in self.dims:
            if k in sel:
                slice_list.append(sel[k])
            else:
                slice_list.append(slice(None))
        return tuple(slice_list)

    def footprint(self, crs: Optional[str] = None) -> Polygon:
        """Returns the footprint of the GeoTensor as a Polygon.

        Args:
            crs (Optional[str], optional): Target coordinate reference system for the footprint.
                If None, returns footprint in the GeoTensor's native CRS. Common formats:
                "EPSG:4326" (WGS84), "EPSG:32630" (UTM), CRS object, or WKT string.
                Use "EPSG:4326" for compatibility with web mapping and GeoJSON.
                Defaults to None.

        Returns:
            Polygon: Shapely Polygon with exactly 5 vertices (4 corners + closing point)
                representing the rectangular footprint. Coordinates are in the specified
                CRS (or native CRS if crs=None).

        Examples:
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>> import numpy as np
            >>>
            >>> # Example 1: Get footprint in native CRS
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)  # UTM transform
            >>> data = np.random.rand(3, 100, 100)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
            >>> footprint_utm = gt.footprint()
            >>> print(f"Bounds (UTM): {footprint_utm.bounds}")
            >>> # (xmin, ymin, xmax, ymax) in UTM meters: (500000, 4649000, 501000, 4650000)

            >>> # Example 2: Transform footprint to WGS84 for web mapping
            >>> footprint_wgs84 = gt.footprint(crs="EPSG:4326")
            >>> print(f"Bounds (WGS84): {footprint_wgs84.bounds}")
            >>> # (lon_min, lat_min, lon_max, lat_max) in degrees
            >>> # Can be used directly in Leaflet, Google Maps, etc.

            >>> # Example 3: Check if rasters overlap
            >>> gt1 = GeoTensor.load_file('image1.tif')
            >>> gt2 = GeoTensor.load_file('image2.tif')
            >>> # Get footprints in common CRS
            >>> fp1 = gt1.footprint(crs="EPSG:4326")
            >>> fp2 = gt2.footprint(crs="EPSG:4326")
            >>> if fp1.intersects(fp2):
            ...     print("Rasters overlap!")
            ...     overlap_area = fp1.intersection(fp2).area
            ...     print(f"Overlap area: {overlap_area} square degrees")

            >>> # Example 4: Export footprint as GeoJSON
            >>> from shapely.geometry import mapping
            >>> import json
            >>> footprint = gt.footprint(crs="EPSG:4326")
            >>> geojson = {
            ...     "type": "Feature",
            ...     "geometry": mapping(footprint),
            ...     "properties": {"name": "Raster extent"}
            ... }
            >>> with open('footprint.geojson', 'w') as f:
            ...     json.dump(geojson, f)

            >>> # Example 5: Check if point is within raster extent
            >>> from shapely.geometry import Point
            >>> point_of_interest = Point(-3.7038, 40.4168)  # Madrid coordinates
            >>> footprint = gt.footprint(crs="EPSG:4326")
            >>> if footprint.contains(point_of_interest):
            ...     print("Point is within raster extent")

            >>> # Example 6: Calculate raster area in square kilometers
            >>> footprint_utm = gt.footprint()  # In UTM (meters)
            >>> area_sqm = footprint_utm.area
            >>> area_sqkm = area_sqm / 1_000_000
            >>> print(f"Raster covers {area_sqkm:.2f} km²")

        Note:
            - Footprint represents the full raster extent, including nodata regions
            - For actual data coverage (excluding nodata), use valid_footprint()
            - Polygon always has rectangular shape (4 corners) even for rotated rasters
            - CRS transformation is performed if target CRS differs from native CRS
            - Footprint is cached-free (computed on-demand from transform and shape)
        """
        pol = window_utils.window_polygon(
            rasterio.windows.Window(
                row_off=0, col_off=0, height=self.shape[-2], width=self.shape[-1]
            ),
            self.transform,
        )
        if (crs is None) or window_utils.compare_crs(self.crs, crs):
            return pol

        return window_utils.polygon_to_crs(pol, self.crs, crs)

    def valid_footprint(
        self, crs: Optional[str] = None, method: str = "all"
    ) -> Union[MultiPolygon, Polygon]:
        """
        vectorizes the valid values of the GeoTensor and returns the footprint as a Polygon.

        Args:
            crs (Optional[str], optional): Coordinate reference system. Defaults to None.
            method (str, optional): "all" or "any" to aggregate the channels of the image. Defaults to "all".

        Returns:
            Polygon or MultiPolygon: footprint of the GeoTensor.
        """
        valid_values = self.values != self.fill_value_default
        if len(valid_values.shape) > 2:
            if method == "all":
                valid_values = np.all(
                    valid_values,
                    axis=tuple(np.arange(0, len(valid_values.shape) - 2).tolist()),
                )
            elif method == "any":
                valid_values = np.any(
                    valid_values,
                    axis=tuple(np.arange(0, len(valid_values.shape) - 2).tolist()),
                )
            else:
                raise NotImplementedError(
                    f"Method {method} to aggregate channels not implemented"
                )

        from georeader import vectorize

        polygons = vectorize.get_polygons(valid_values, transform=self.transform)
        if len(polygons) == 0:
            raise ValueError("GeoTensor has no valid values")
        elif len(polygons) == 1:
            pol = polygons[0]
        else:
            pol = MultiPolygon(polygons)
        if crs is None:
            return pol

        return window_utils.polygon_to_crs(pol, self.crs, crs)

    def __repr__(self) -> str:
        return f""" 
         Transform: {self.transform}
         Shape: {self.shape}
         Resolution: {self.res}
         Bounds: {self.bounds}
         CRS: {self.crs}
         fill_value_default: {self.fill_value_default}
        """

    def pad(
        self,
        pad_width: Union[Dict[str, Tuple[int, int]], List[Tuple[int, int]]],
        mode: str = "constant",
        **kwargs: Any,
    ) -> "GeoTensor":
        """
        Add padding around the GeoTensor, updating georeferencing.

        Expands the spatial dimensions of the GeoTensor by adding pixels on each
        side, while automatically updating the geotransform to maintain correct
        geographic positioning. The new pixels are filled according to the
        specified mode.

        Padding is essential for:
        - CNN inference requiring fixed input sizes
        - Avoiding edge artifacts in convolution operations
        - Creating uniform tile sizes for batch processing

        Padding behavior::

            pad_width = {"x": (2, 3), "y": (1, 4)}

            Original (5×5):           Padded (10×10):
            ┌─────────────┐           ┌ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ┐
            │ █ █ █ █ █ │           │   ← 1 row top           │
            │ █ █ █ █ █ │    →      │ ┌─────────────┐          │
            │ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
            │ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
            │ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
            └─────────────┘           │ │ █ █ █ █ █ │          │
                                      │ │ █ █ █ █ █ │          │
                                      │ └─────────────┘          │
                                      │   ← 4 rows bottom       │
                                      └ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ┘
                                        ↑           ↑
                                        2 cols      3 cols
                                        left        right

        Args:
            pad_width (Union[Dict[str, Tuple[int, int]], List[Tuple[int, int]]]):
                Padding specification. Can be:

                - Dict: Dimension names to (before, after) padding::

                    {"x": (left, right), "y": (top, bottom)}
                    {"x": (10, 10), "y": (5, 5), "band": (0, 0)}

                - List/Tuple: Padding per dimension in order of self.dims::

                    # For 3D GeoTensor with dims (band, y, x):
                    [(band_before, band_after), (y_before, y_after), (x_before, x_after)]

            mode (str, optional): Padding mode (see numpy.pad). Options:

                - "constant": Fill with constant value (default)
                - "edge": Replicate edge pixels
                - "reflect": Mirror reflection at edges
                - "wrap": Wrap around (toroidal)

                Defaults to "constant".

            **kwargs: Additional arguments passed to np.pad:

                - constant_values: Fill value for "constant" mode.
                  Defaults to self.fill_value_default.

        Returns:
            GeoTensor: Padded GeoTensor with:
                - Updated transform: origin shifted by padding amount
                - Shape: increased by sum of padding in each dimension
                - Same CRS and fill_value_default

        Raises:
            ValueError: If list-style pad_width has wrong length.

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>>
            >>> # Create 3-band 100x100 GeoTensor
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> data = np.random.rand(3, 100, 100).astype(np.float32)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630", fill_value_default=0)
            >>> print(f"Original: {gt.shape}")
            Original: (3, 100, 100)
            >>>
            >>> # Pad spatial dimensions (dict style)
            >>> padded = gt.pad({"x": (32, 32), "y": (32, 32)})
            >>> print(f"Padded: {padded.shape}")
            Padded: (3, 164, 164)
            >>> # Transform shifted: origin moved by -32 pixels (320m)
            >>> print(f"New origin: ({padded.transform.c}, {padded.transform.f})")
            New origin: (499680.0, 4650320.0)
            >>>
            >>> # Pad with list style (must match dims order)
            >>> padded = gt.pad([(0, 0), (16, 16), (16, 16)])
            >>> print(f"Padded: {padded.shape}")
            Padded: (3, 132, 132)
            >>>
            >>> # Asymmetric padding for edge tiles
            >>> edge_padded = gt.pad({"x": (0, 50), "y": (0, 25)})
            >>> print(f"Edge padded: {edge_padded.shape}")
            Edge padded: (3, 125, 150)
            >>>
            >>> # Edge replication instead of constant fill
            >>> padded_edge = gt.pad({"x": (10, 10), "y": (10, 10)}, mode="edge")

        Note:
            - Transform is updated so padded region has correct georeferencing
            - For CNN inference, also consider `read_from_window` with boundless=True
            - Use pad_width={("band": (0, 0), "x": ..., "y": ...} to be explicit

        See Also:
            - `pad_array`: Lower-level method with dict-only interface
            - `read_from_window`: Read with implicit padding for out-of-bounds
            - `window_utils.pad_window`: Expand window coordinates
        """
        if isinstance(pad_width, list) or isinstance(pad_width, tuple):
            if len(pad_width) != len(self.dims):
                raise ValueError(
                    f"Expected {len(self.dims)} pad widths found {len(pad_width)}"
                )
            pad_width_dict = {}
            for i, k in enumerate(self.dims):
                pad_width_dict[k] = pad_width[i]
        else:
            pad_width_dict = pad_width
        return self.pad_array(pad_width_dict, mode=mode, **kwargs)

    def pad_array(
        self,
        pad_width: Dict[str, Tuple[int, int]],
        mode: str = "constant",
        constant_values: Optional[Any] = None,
    ) -> Self:
        """
        Pad the GeoTensor.

        Args:
            pad_width (_type_, optional):  dictionary with Tuple to pad for each dimension
                `{"x": (pad_x_0, pad_x_1), "y": (pad_y_0, pad_y_1)}`.
            mode (str, optional): pad mode (see np.pad or torch.nn.functional.pad). Defaults to "constant".
            constant_values (Any, optional): _description_. Defaults to `self.fill_value_default`.

        Returns:
            GeoTensor: padded GeoTensor.

        Examples:
            >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
            >>> gt.pad_array({"x": (10, 10), "y": (10, 10)})
            >>> assert gt.shape == (3, 120, 120)
        """
        if constant_values is None and mode == "constant":
            if self.fill_value_default is None:
                raise ValueError(
                    f"Mode constant either requires constant_values passed or fill_value_default not None in current GeoTensor"
                )
            constant_values = self.fill_value_default

        pad_list_np = []
        for k in self.dims:
            if k in pad_width:
                pad_list_np.append(pad_width[k])
            else:
                pad_list_np.append((0, 0))

        kwargs_extra = {}
        if mode == "constant":
            kwargs_extra["constant_values"] = constant_values
        values_new = np.pad(self.values, tuple(pad_list_np), mode=mode, **kwargs_extra)

        # Compute the new transform
        slices_window = []
        for k in ["y", "x"]:
            size = self.width if (k == "x") else self.height
            if k in pad_width:
                slices_window.append(slice(-pad_width[k][0], size + pad_width[k][1]))
            else:
                slices_window.append(slice(0, size))

        window_current = rasterio.windows.Window.from_slices(
            *slices_window, boundless=True
        )
        transform_current = rasterio.windows.transform(
            window_current, transform=self.transform
        )
        return GeoTensor(
            values_new,
            transform_current,
            self.crs,
            self.fill_value_default,
            attrs=self.attrs,
        )

    def resize(
        self,
        output_shape: Optional[Tuple[int, int]] = None,
        resolution_dst: Optional[Tuple[float, float]] = None,
        anti_aliasing: bool = True,
        anti_aliasing_sigma: Optional[Union[float, np.ndarray]] = None,
        interpolation: Optional[str] = "bilinear",
        mode_pad: str = "constant",
    ) -> Self:
        """
        Resize the geotensor to match a certain size output_shape. This function works with GeoTensors of 2D, 3D and 4D.
        The geoinformation of the output tensor is changed accordingly.

        Args:
            output_shape: output spatial shape if None resolution_dst must be provided. If not provided,
                the output shape is computed from the resolution_dst rounding to the closest integer.
            resolution_dst: output resolution if None output_shape must be provided.
            anti_aliasing: Whether to apply a Gaussian filter to smooth the image prior to downsampling
            anti_aliasing_sigma:  anti_aliasing_sigma : {float}, optional
                Standard deviation for Gaussian filtering used when anti-aliasing.
                By default, this value is chosen as (s - 1) / 2 where s is the
                downsampling factor, where s > 1
            interpolation: Algorithm used for resizing: 'nearest' | 'bilinear' | 'bicubic'
            mode_pad: mode pad for resize function

        Returns:
             resized GeoTensor

        Examples:
            >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
            >>> resized = gt.resize((50, 50))
            >>> assert resized.shape == (3, 50, 50)
            >>> assert resized.res == (2*gt.res[0], 2*gt.res[1])
        """
        input_shape = self.shape
        spatial_shape = input_shape[-2:]
        resolution_or = self.res

        if output_shape is None:
            assert (
                resolution_dst is not None
            ), f"Can't have output_shape and resolution_dst as None"
            output_shape = int(
                round(spatial_shape[0] * resolution_or[0] / resolution_dst[0])
            ), int(round(spatial_shape[1] * resolution_or[1] / resolution_dst[1]))
        else:
            assert (
                resolution_dst is None
            ), f"Both output_shape and resolution_dst can't be provided"
            assert (
                len(output_shape) == 2
            ), f"Expected output shape to be the spatial dimensions found: {output_shape}"
            resolution_dst = (
                spatial_shape[0] * resolution_or[0] / output_shape[0],
                spatial_shape[1] * resolution_or[1] / output_shape[1],
            )

        # Compute output transform
        transform_scale = rasterio.Affine.scale(
            resolution_dst[0] / resolution_or[0], resolution_dst[1] / resolution_or[1]
        )
        transform = self.transform * transform_scale

        from skimage.transform import resize

        # https://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.resize
        output_tensor = np.ndarray(input_shape[:-2] + output_shape, dtype=self.dtype)
        if len(input_shape) == 4:
            for i, j in product(range(0, input_shape[0]), range(0, input_shape[1])):
                if (
                    (not anti_aliasing)
                    or (anti_aliasing_sigma is None)
                    or isinstance(anti_aliasing_sigma, numbers.Number)
                ):
                    anti_aliasing_sigma_iter = anti_aliasing_sigma
                else:
                    anti_aliasing_sigma_iter = anti_aliasing_sigma[i, j]
                output_tensor[i, j] = resize(
                    self.values[i, j],
                    output_shape,
                    order=ORDERS[interpolation],
                    anti_aliasing=anti_aliasing,
                    preserve_range=False,
                    cval=self.fill_value_default,
                    mode=mode_pad,
                    anti_aliasing_sigma=anti_aliasing_sigma_iter,
                )
        elif len(input_shape) == 3:
            for i in range(0, input_shape[0]):
                if (
                    (not anti_aliasing)
                    or (anti_aliasing_sigma is None)
                    or isinstance(anti_aliasing_sigma, numbers.Number)
                ):
                    anti_aliasing_sigma_iter = anti_aliasing_sigma
                else:
                    anti_aliasing_sigma_iter = anti_aliasing_sigma[i]
                output_tensor[i] = resize(
                    self.values[i],
                    output_shape,
                    order=ORDERS[interpolation],
                    anti_aliasing=anti_aliasing,
                    preserve_range=False,
                    cval=self.fill_value_default,
                    mode=mode_pad,
                    anti_aliasing_sigma=anti_aliasing_sigma_iter,
                )
        else:
            output_tensor[...] = resize(
                self.values,
                output_shape,
                order=ORDERS[interpolation],
                anti_aliasing=anti_aliasing,
                preserve_range=False,
                cval=self.fill_value_default,
                mode=mode_pad,
                anti_aliasing_sigma=anti_aliasing_sigma,
            )

        return GeoTensor(
            output_tensor,
            transform=transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def transpose(self, axes=None) -> Self:
        """
        Permute the dimensions of the GeoTensor while keeping the spatial dimensions at the end.

        Args:
            axes (tuple, optional): If specified, it must be a tuple or list of axes. The last two
                values must be the original spatial dimensions indices (ndim-2, ndim-1).
                If None, the non-spatial dimensions are reversed while spatial dimensions remain at the end.

        Returns:
            GeoTensor: A view of the array with dimensions transposed.

        Raises:
            ValueError: If the spatial dimensions are moved from their last positions.

        Examples:
            >>> gt = GeoTensor(np.random.rand(3, 4, 100, 200), transform, crs)
            >>> # Shape is (3, 4, 100, 200)
            >>> gt_t = gt.transpose()
            >>> # Shape is now (4, 3, 100, 200)
            >>>
            >>> # You can also specify axes explicitly:
            >>> gt_t = gt.transpose((1, 0, 2, 3))  # Valid: spatial dims remain at end
            >>> # But this would raise an error:
            >>> # gt.transpose((0, 2, 1, 3))  # Invalid: spatial dims must stay at end
        """
        ndim = len(self.shape)

        if ndim <= 2:
            # Nothing meaningful to transpose for arrays with only spatial dimensions
            return self.copy()

        # Original spatial dimensions indices
        y_dim = ndim - 2
        x_dim = ndim - 1

        if axes is None:
            # Reverse all dimensions except the spatial ones which stay at the end
            non_spatial_axes = list(range(ndim - 2))
            non_spatial_axes.reverse()
            axes = tuple(non_spatial_axes + [y_dim, x_dim])
        else:
            # Convert to tuple if necessary
            axes = tuple(axes)

            # Check if axes has the right length
            if len(axes) != ndim:
                raise ValueError(
                    f"axes should contain {ndim} dimensions, got {len(axes)}"
                )

            # Check if the last two values in axes are the spatial dimensions
            if axes[-2:] != (y_dim, x_dim):
                raise ValueError(
                    "Cannot change the position of spatial dimensions. "
                    f"The last two axes must be {y_dim} and {x_dim}."
                )

        # Perform the transpose
        transposed_values = np.transpose(self.values, axes)

        return GeoTensor(
            transposed_values,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )

    def validmask(self) -> Self:
        """
        Returns a mask of the valid values of the GeoTensor. The mask is a boolean array
        with the same shape as the GeoTensor values, where True indicates valid values and
        False indicates invalid values.
        The mask is created by comparing the values of the GeoTensor with the `self.fill_value_default`.

        Returns:
            Self: GeoTensor with the valid boolean mask.
        """
        if self.fill_value_default is None:
            return GeoTensor(
                np.ones(self.shape, dtype=bool),
                transform=self.transform,
                crs=self.crs,
                fill_value_default=self.fill_value_default,
                attrs=self.attrs,
            )
        return GeoTensor(
            values=self.values != self.fill_value_default,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    def invalidmask(self) -> Self:
        """
        Returns a mask of the invalid values of the GeoTensor. The mask is a boolean array
        with the same shape as the GeoTensor values, where True indicates invalid values and
        False indicates valid values.
        The mask is created by comparing the values of the GeoTensor with the `self.fill_value_default`.

        Returns:
            Self: GeoTensor with the invalid boolean mask.
        """
        if self.fill_value_default is None:
            return GeoTensor(
                np.zeros(self.shape, dtype=bool),
                transform=self.transform,
                crs=self.crs,
                fill_value_default=self.fill_value_default,
            )
        return GeoTensor(
            values=self.values == self.fill_value_default,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=False,
            attrs=self.attrs,
        )

    @classmethod
    def load_file(
        cls,
        path: str,
        fs: Optional[Any] = None,
        load_tags: bool = False,
        load_descriptions: bool = False,
        rio_env_options: Optional[Dict[str, str]] = None,
    ) -> Self:
        """
        Load a GeoTensor from a file. It uses rasterio to read the data. This function
        loads all the data in memory. For a lazy loading reader use `georeader.rasterio_reader`.

        Args:
            path (str): Path to the file.
            fs (Optional[Any], optional): fsspec.Filesystem object. Defaults to None.
            load_descriptions (bool, optional): If True, load the description of the image. Defaults to False.
            load_tags (bool, optional): If True, load the tags of the image. Defaults to False.
            rio_env_options (Optional[Dict[str, str]], optional): Rasterio environment options. Defaults to None.

        Returns:
            GeoTensor: GeoTensor object with the loaded data.
        """

        if fs is not None:
            if load_descriptions:
                raise NotImplementedError(
                    """Description loading not supported with `fsspec`. This is because
            the `descriptions` attribute cannote be loaded from a byte stream. This is a limitation of `rasterio`.
            The issue is related to how `rasterio.io.MemoryFile` handles band descriptions
            compared to direct file access. This is a known limitation when working
            with in-memory file representations in GDAL (which `rasterio` uses under
            the hood). If you need to load descriptions, you can use `georeader.rasterio_reader`
            class."""
                )

            with fs.open(path, "rb") as fh:
                return cls.load_bytes(
                    fh.read(), load_tags=load_tags, rio_env_options=rio_env_options
                )

        tags = None
        descriptions = None
        rio_env_options = (
            RIO_ENV_OPTIONS_DEFAULT if rio_env_options is None else rio_env_options
        )
        with rasterio.Env(**get_rio_options_path(rio_env_options, path)):
            with rasterio.open(path) as src:
                data = src.read()
                transform = src.transform
                crs = src.crs
                fill_value_default = src.nodata
                if load_tags:
                    tags = src.tags()
                if load_descriptions:
                    descriptions = tuple(src.descriptions)

        attrs = {}
        if tags is not None:
            attrs["tags"] = tags

        if descriptions is not None:
            attrs["descriptions"] = descriptions

        return cls(
            data, transform, crs, fill_value_default=fill_value_default, attrs=attrs
        )

    @classmethod
    def load_bytes(
        cls,
        bytes_read: Union[bytes, bytearray, memoryview],
        load_tags: bool = False,
        rio_env_options: Optional[Dict[str, str]] = None,
    ) -> Self:
        """
        Load a GeoTensor from a byte stream. It uses rasterio to read the data.


        Args:
            bytes_read (Union[bytes, bytearray, memoryview]): Byte stream to read.
            load_tags (bool, optional): if True, load the tags of the image. Defaults to False.
            rio_env_options (Optional[Dict[str, str]], optional): Rasterio environment options. Defaults to None.

        Returns:
            __class__: GeoTensor object with the loaded data.

        Note:
            The `descriptions` attribute cannote be loaded from a byte stream. This is a limitation of `rasterio`.
            The issue is related to how `rasterio.io.MemoryFile` handles band descriptions
            compared to direct file access. This is a known limitation when working
            with in-memory file representations in GDAL (which `rasterio` uses under
            the hood). If you need to load descriptions, you should use `georeader.rasterio_reader`
            class.
        """
        import rasterio.io

        tags = None
        rio_env_options = (
            RIO_ENV_OPTIONS_DEFAULT if rio_env_options is None else rio_env_options
        )
        with rasterio.Env(**rio_env_options):
            with rasterio.io.MemoryFile(bytes_read) as mem:
                with mem.open() as src:
                    data = src.read()
                    transform = src.transform
                    crs = src.crs
                    fill_value_default = src.nodata
                    if load_tags:
                        tags = src.tags()

        attrs = {}
        if tags is not None:
            attrs["tags"] = tags

        return cls(
            data, transform, crs, fill_value_default=fill_value_default, attrs=attrs
        )

    def write_from_window(self, data: np.ndarray, window: rasterio.windows.Window):
        """
        Writes array to GeoTensor values object at the given window position. If window surpasses the bounds of this
        object it crops the data to fit the object.

        Args:
            data: Tensor to write. Expected: spatial dimensions `window.width`, `window.height`. Rest: same as `self`
            window: Window object that specifies the spatial location to write the data

        Examples:
            >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
            >>> data = np.random.rand(3, 50, 50)
            >>> window = rasterio.windows.Window(col_off=7, row_off=9, width=50, height=50)
            >>> gt.write_from_window(data, window)

        """
        window_data = rasterio.windows.Window(
            col_off=0, row_off=0, width=self.width, height=self.height
        )
        if not rasterio.windows.intersect(window, window_data):
            return

        assert data.shape[-2:] == (
            window.height,
            window.width,
        ), f"window {window} has different shape than data {data.shape}"
        assert (
            data.shape[:-2] == self.shape[:-2]
        ), f"Dimension of data in non-spatial channels found {data.shape} expected: {self.shape}"

        slice_dict, pad_width = window_utils.get_slice_pad(window_data, window)
        slice_list = self._slice_tuple(slice_dict)
        # need_pad = any(p != 0 for p in pad_width["x"] + pad_width["y"])

        slice_data_spatial_x = slice(
            pad_width["x"][0], None if pad_width["x"][1] == 0 else -pad_width["x"][1]
        )
        slice_data_spatial_y = slice(
            pad_width["y"][0], None if pad_width["y"][1] == 0 else -pad_width["y"][1]
        )
        slice_data = self._slice_tuple(
            {"x": slice_data_spatial_x, "y": slice_data_spatial_y}
        )
        self.values[slice_list] = data[slice_data]

    def read_from_window(
        self, window: rasterio.windows.Window, boundless: bool = True
    ) -> Self:
        """
        Extract a spatial subset using a rasterio Window.

        Reads data from the GeoTensor corresponding to the specified pixel window,
        with optional boundless reading that pads out-of-bounds regions with the
        fill value. This method provides the same interface as rasterio's windowed
        reading, enabling seamless transitions between lazy file readers and
        in-memory GeoTensors.

        Boundless vs Bounded Reading::

            Window extends beyond data:     boundless=True:      boundless=False:
            ┌───────────────┐                 ┌─────────┐           ┌─────┐
            │ ┌─────────┐   │                 │▒▒▒▒▒▒▒▒▒│           │     │
            │ │  DATA   │   │                 │▒▒▒███▒▒▒│           │ ███ │
            │ │         │   │  ─────────►    │▒▒▒███▒▒▒│           │ ███ │
            │ └─────────┘   │                 │▒▒▒▒▒▒▒▒▒│           └─────┘
            └───────────────┘                 └─────────┘           Returns only
                 Window                   Padded with            intersection
                 request                  fill_value

        Args:
            window (rasterio.windows.Window): Pixel window defining the region
                to extract. Format: Window(col_off, row_off, width, height).
                May have negative offsets or extend beyond data bounds if
                boundless=True.
            boundless (bool, optional): If True, allows reading windows that
                extend beyond data boundaries. Out-of-bounds regions are filled
                with `self.fill_value_default`. If False, returns only the
                intersection with the data extent. Defaults to True.

        Returns:
            GeoTensor: Extracted data with:
                - Shape: (C, window.height, window.width) if boundless=True,
                  or smaller if boundless=False and window exceeds bounds
                - Transform: Updated to window origin
                - Same CRS and fill_value_default

        Raises:
            rasterio.windows.WindowError: If window does not intersect the data
                at all (no overlap).

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> import rasterio.windows
            >>> from georeader import GeoTensor
            >>>
            >>> # Create test GeoTensor
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> data = np.arange(300).reshape(3, 10, 10).astype(np.float32)
            >>> gt = GeoTensor(data, transform, crs="EPSG:32630", fill_value_default=-1)
            >>> print(f"Data shape: {gt.shape}")
            Data shape: (3, 10, 10)
            >>>
            >>> # Read window fully within data
            >>> window = rasterio.windows.Window(2, 3, 5, 4)  # col_off, row_off, width, height
            >>> subset = gt.read_from_window(window)
            >>> print(f"Subset shape: {subset.shape}")
            Subset shape: (3, 4, 5)
            >>>
            >>> # Boundless read: window extends beyond data
            >>> window_edge = rasterio.windows.Window(7, 7, 6, 6)  # Extends 3 pixels past edge
            >>> padded = gt.read_from_window(window_edge, boundless=True)
            >>> print(f"Padded shape: {padded.shape}")
            Padded shape: (3, 6, 6)
            >>> # Corners filled with fill_value_default (-1)
            >>> print(f"Has padding: {(padded.values == -1).any()}")
            Has padding: True
            >>>
            >>> # Bounded read: only get intersection
            >>> clipped = gt.read_from_window(window_edge, boundless=False)
            >>> print(f"Clipped shape: {clipped.shape}")
            Clipped shape: (3, 3, 3)
            >>>
            >>> # CNN tiling: process image in overlapping windows
            >>> tile_size = 4
            >>> stride = 2
            >>> for row in range(0, gt.height, stride):
            ...     for col in range(0, gt.width, stride):
            ...         tile_window = rasterio.windows.Window(col, row, tile_size, tile_size)
            ...         tile = gt.read_from_window(tile_window, boundless=True)
            ...         # Process tile with CNN...

        Note:
            - Boundless reading is essential for processing edge tiles in CNNs
            - Transform is always updated to match the output data origin
            - Use `get_slice_pad` directly for manual slice/pad control
            - Compatible with RasterioReader interface for polymorphic code

        See Also:
            - `isel`: Dictionary-based dimension indexing
            - `window_utils.get_slice_pad`: Compute slice and padding parameters
            - `RasterioReader.read_from_window`: Lazy file-based equivalent
        """

        window_data = rasterio.windows.Window(
            col_off=0, row_off=0, width=self.width, height=self.height
        )
        if boundless:
            slice_dict, pad_width = window_utils.get_slice_pad(window_data, window)
            need_pad = any(p != 0 for p in pad_width["x"] + pad_width["y"])
            X_sliced = self.isel(slice_dict)
            if need_pad:
                X_sliced = X_sliced.pad(
                    pad_width=pad_width,
                    mode="constant",
                    constant_values=self.fill_value_default,
                )
            return X_sliced
        else:
            window_read = rasterio.windows.intersection(window, window_data)
            slice_y, slice_x = window_read.toslices()
            slice_dict = {"x": slice_x, "y": slice_y}
            slices_ = self._slice_tuple(slice_dict)
            transform_current = rasterio.windows.transform(
                window_read, transform=self.transform
            )
            return GeoTensor(
                self.values[slices_],
                transform_current,
                self.crs,
                self.fill_value_default,
                attrs=self.attrs,
            )

    @classmethod
    def stack(cls, geotensors: List[Self]) -> Self:
        """
        Stack multiple GeoTensors along a new leading dimension.

        Creates a new GeoTensor with an additional dimension at the front,
        containing all input GeoTensors. All inputs must have identical
        georeferencing (transform, CRS) and shape.

        This is useful for:
        - Creating time-series stacks from multiple acquisitions
        - Batching for CNN inference
        - Multi-temporal analysis

        Stacking behavior::

            Input: 3 GeoTensors each (3, H, W)

            gt1: (3, 100, 100)  ─┐
            gt2: (3, 100, 100)  ─┼──► stacked: (3, 3, 100, 100)
            gt3: (3, 100, 100)  ─┘     ↑
                                    new time/batch dim

        Args:
            geotensors (List[GeoTensor]): List of GeoTensors to stack. Requirements:

                - Non-empty list
                - All must have identical shape
                - All must have same transform (same_extent)
                - All must have same CRS
                - All must have same fill_value_default

        Returns:
            GeoTensor: Stacked GeoTensor with:

                - Shape: (len(geotensors),) + original_shape
                - Transform: Same as input GeoTensors
                - CRS: Same as input GeoTensors
                - fill_value_default: Same as input GeoTensors
                - attrs: Copied from first GeoTensor

        Raises:
            AssertionError: If list is empty.
            AssertionError: If GeoTensors have different extents.
            AssertionError: If GeoTensors have different shapes.
            AssertionError: If GeoTensors have different fill_value_default.

        Examples:
            >>> import numpy as np
            >>> import rasterio
            >>> from georeader import GeoTensor
            >>>
            >>> # Create multiple GeoTensors (e.g., different dates)
            >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
            >>> gt_jan = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
            >>> gt_feb = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
            >>> gt_mar = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
            >>>
            >>> # Stack into time series
            >>> time_series = GeoTensor.stack([gt_jan, gt_feb, gt_mar])
            >>> print(f"Time series shape: {time_series.shape}")
            Time series shape: (3, 3, 100, 100)
            >>> # dims: (time, band, y, x)
            >>>
            >>> # Access individual time steps
            >>> jan_data = time_series[0]  # Returns GeoTensor with shape (3, 100, 100)
            >>>
            >>> # Compute temporal median (across time dimension)
            >>> median = np.nanmedian(time_series.values, axis=0)
            >>> # median shape: (3, 100, 100)
            >>>
            >>> # Stack single GeoTensor (adds dimension)
            >>> single_stacked = GeoTensor.stack([gt_jan])
            >>> print(f"Single stacked: {single_stacked.shape}")
            Single stacked: (1, 3, 100, 100)

        Note:
            - All GeoTensors must have identical georeferencing
            - Use `concatenate` to join along existing dimension instead
            - For memory efficiency with large stacks, consider lazy loading
            - The first GeoTensor's attrs are preserved

        See Also:
            - `concatenate`: Join along existing axis (no new dimension)
            - `same_extent`: Check if GeoTensors have matching georeferencing
            - `numpy.stack`: Underlying operation for values
        """
        assert len(geotensors) > 0, "Empty list provided can't concat"

        if len(geotensors) == 1:
            gt = geotensors[0]
            return GeoTensor(
                gt.values[np.newaxis],
                transform=gt.transform,
                crs=gt.crs,
                fill_value_default=gt.fill_value_default,
                attrs=gt.attrs,
            )

        first_geotensor = geotensors[0]
        array_out = np.zeros(
            (len(geotensors),) + first_geotensor.shape, dtype=first_geotensor.dtype
        )
        array_out[0] = first_geotensor.values

        for i, geo in enumerate(geotensors[1:]):
            assert geo.same_extent(first_geotensor), f"Different size in concat {i+1}"
            assert (
                geo.shape == first_geotensor.shape
            ), f"Different shape in concat {i+1}"
            assert (
                geo.fill_value_default == first_geotensor.fill_value_default
            ), "Different fill_value_default in concat"
            array_out[i + 1] = geo.values

        return cls(
            array_out,
            transform=first_geotensor.transform,
            crs=first_geotensor.crs,
            fill_value_default=first_geotensor.fill_value_default,
            attrs=first_geotensor.attrs,
        )

    @classmethod
    def concatenate(cls, geotensors: List[Self], axis: int = 0) -> Self:
        """
        Concatenates a list of geotensors along a given axis, assert that all of them has same shape, transform and crs.

        Args:
            geotensors: list of geotensors to concat. All with same shape, transform and crs.
            axis: axis to concatenate. Must be less than the number of dimensions of the geotensors minus 2.
                default is 0.

        Returns:
            geotensor with extra dim at the front: (len(geotensors),) + shape

        Examples:
            >>> gt1 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
            >>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
            >>> gt3 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
            >>> gt = concatenate([gt1, gt2, gt3], axis=0)
            >>> assert gt.shape == (9, 100, 100)
        """
        assert len(geotensors) > 0, "Empty list provided can't concat"

        if len(geotensors) == 1:
            return geotensors[0].copy()

        first_geotensor = geotensors[0]

        # Assert the axis is NOT an spatial axis
        assert (
            axis < len(first_geotensor.shape) - 2
        ), f"Can't concatenate along spatial axis"

        for i, geo in enumerate(geotensors[1:]):
            assert geo.same_extent(first_geotensor), f"Different extent in concat {i+1}"
            assert (
                geo.shape == first_geotensor.shape
            ), f"Different shape in concat {i+1}"
            assert (
                geo.fill_value_default == first_geotensor.fill_value_default
            ), "Different fill_value_default in concat"

        array_out = np.concatenate([gt.values for gt in geotensors], axis=axis)

        return cls(
            array_out,
            transform=first_geotensor.transform,
            crs=first_geotensor.crs,
            fill_value_default=first_geotensor.fill_value_default,
            attrs=first_geotensor.attrs,
        )

values property

Return a view of the array (memory shared with original)

__add__(other)

Add a value or array to this GeoTensor element-wise.

Supports broadcasting with scalars, numpy arrays, or other GeoTensors. When adding two GeoTensors, they must have the same spatial extent (matching transform, CRS, and spatial dimensions).

Broadcasting Rules: - Scalar: Added to every pixel - Array: Must be broadcastable to self.values shape - GeoTensor: Must have identical georeferencing (same_extent)

Parameters:

Name	Type	Description	Default
other	Union[Number, ndarray, GeoTensor]	Value to add. Can be: - Scalar (int, float): Added to all pixels - np.ndarray: Must be broadcastable with self.values - GeoTensor: Must have same spatial extent	required

Returns:

Name	Type	Description
GeoTensor	Self	New GeoTensor with same transform, CRS, and fill_value_default. Shape matches self.shape (or broadcast result for arrays).

Raises:

Type	Description
ValueError	If other is a GeoTensor and georeferencing doesn't match.

Examples:

>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>>
>>> # Create sample GeoTensor (3 bands, 100x100 pixels)
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>>
>>> # Add scalar to all pixels
>>> gt_offset = gt + 0.1
>>> print(gt_offset.shape)  # (3, 100, 100)
>>>
>>> # Add per-band offset using broadcasting
>>> band_offsets = np.array([0.1, 0.2, 0.3])[:, None, None]  # Shape: (3, 1, 1)
>>> gt_adjusted = gt + band_offsets
>>>
>>> # Add two GeoTensors (must have same extent)
>>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
>>> gt_sum = gt + gt2
>>>
>>> # Error case: mismatched georeferencing
>>> gt_different = GeoTensor(data, rasterio.Affine(20, 0, 0, 0, -20, 0), crs="EPSG:4326")
>>> # gt + gt_different  # Raises ValueError

Note

• Result inherits transform, CRS, and fill_value_default from self
• For GeoTensor addition, use read.read_reproject_like(other, self) to align georeferencing before adding

See Also

• same_extent: Check if two GeoTensors have matching georeferencing
• read.read_reproject_like: Reproject one GeoTensor to match another

Source code in georeader/geotensor.py

def __add__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Add a value or array to this GeoTensor element-wise.

    Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
    When adding two GeoTensors, they must have the same spatial extent
    (matching transform, CRS, and spatial dimensions).

    Broadcasting Rules:
    - Scalar: Added to every pixel
    - Array: Must be broadcastable to self.values shape
    - GeoTensor: Must have identical georeferencing (same_extent)

    Args:
        other (Union[numbers.Number, np.ndarray, GeoTensor]): Value to add. Can be:
            - Scalar (int, float): Added to all pixels
            - np.ndarray: Must be broadcastable with self.values
            - GeoTensor: Must have same spatial extent

    Returns:
        GeoTensor: New GeoTensor with same transform, CRS, and fill_value_default.
            Shape matches self.shape (or broadcast result for arrays).

    Raises:
        ValueError: If other is a GeoTensor and georeferencing doesn't match.

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>>
        >>> # Create sample GeoTensor (3 bands, 100x100 pixels)
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> data = np.random.rand(3, 100, 100)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
        >>>
        >>> # Add scalar to all pixels
        >>> gt_offset = gt + 0.1
        >>> print(gt_offset.shape)  # (3, 100, 100)
        >>>
        >>> # Add per-band offset using broadcasting
        >>> band_offsets = np.array([0.1, 0.2, 0.3])[:, None, None]  # Shape: (3, 1, 1)
        >>> gt_adjusted = gt + band_offsets
        >>>
        >>> # Add two GeoTensors (must have same extent)
        >>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
        >>> gt_sum = gt + gt2
        >>>
        >>> # Error case: mismatched georeferencing
        >>> gt_different = GeoTensor(data, rasterio.Affine(20, 0, 0, 0, -20, 0), crs="EPSG:4326")
        >>> # gt + gt_different  # Raises ValueError

    Note:
        - Result inherits transform, CRS, and fill_value_default from self
        - For GeoTensor addition, use `read.read_reproject_like(other, self)`
          to align georeferencing before adding

    See Also:
        - `same_extent`: Check if two GeoTensors have matching georeferencing
        - `read.read_reproject_like`: Reproject one GeoTensor to match another
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for addition. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values + other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

__and__(other)

Perform bitwise AND operation between two GeoTensors. The georeferencing must match.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	GeoTensor or array-like to AND with.	required

Raises:

Type	Description
ValueError	if the georeferencing does not match when other is a GeoTensor.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with the result of the bitwise AND operation.

Source code in georeader/geotensor.py

def __and__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Perform bitwise AND operation between two GeoTensors. The georeferencing must match.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): GeoTensor or array-like to AND with.

    Raises:
        ValueError: if the georeferencing does not match when other is a GeoTensor.

    Returns:
        GeoTensor: GeoTensor with the result of the bitwise AND operation.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for bitwise AND. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values & other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

__array__(dtype=None)

Convert the GeoTensor to a standard NumPy array.

This method is called by np.asarray() and most NumPy functions to get the underlying NumPy array representation of this object.

Parameters:

Name	Type	Description	Default
dtype	Optional[dtype]	The desired data type for the returned array. If None, the array's current dtype is preserved.	None

Returns:

Type	Description
ndarray	np.ndarray: A NumPy array view of this GeoTensor.

Source code in georeader/geotensor.py

def __array__(self, dtype: Optional[np.dtype] = None) -> np.ndarray:
    """
    Convert the GeoTensor to a standard NumPy array.

    This method is called by np.asarray() and most NumPy functions to get
    the underlying NumPy array representation of this object.

    Args:
        dtype (Optional[np.dtype]): The desired data type for the returned array.
                                   If None, the array's current dtype is preserved.

    Returns:
        np.ndarray: A NumPy array view of this GeoTensor.
    """
    return np.asarray(self.view(np.ndarray), dtype=dtype)

__array_finalize__(obj)

Initialize attributes when a new GeoTensor is created from an existing array.

This method is called whenever a new array object is created from an existing array (e.g., through slicing, view casting, or copy operations).

Parameters:

Name	Type	Description	Default
obj	Optional[ndarray]	The array object from which the new array is created. Can be None if the array is being created from scratch.	required

Source code in georeader/geotensor.py

def __array_finalize__(self, obj: Optional[Union[np.ndarray, Self]]) -> None:
    """
    Initialize attributes when a new GeoTensor is created from an existing array.

    This method is called whenever a new array object is created from an existing array
    (e.g., through slicing, view casting, or copy operations).

    Args:
        obj (Optional[np.ndarray]): The array object from which the new array is created.
                                   Can be None if the array is being created from scratch.
    """
    if obj is None:
        return

    if hasattr(obj, "transform"):
        self.transform: rasterio.Affine = getattr(obj, "transform", None)
    if hasattr(obj, "crs"):
        self.crs = getattr(obj, "crs", None)
    if hasattr(obj, "fill_value_default"):
        self.fill_value_default = getattr(obj, "fill_value_default", None)
    if hasattr(obj, "attrs"):
        self.attrs = getattr(obj, "attrs", None)

__array_ufunc__(ufunc, method, *inputs, **kwargs)

Handle NumPy universal functions applied to this GeoTensor.

This method is called when a NumPy universal function (ufunc) is applied to the GeoTensor. It converts GeoTensor inputs to NumPy arrays, applies the ufunc, and converts array results back to GeoTensor objects with the same geospatial metadata.

Parameters:

Name	Type	Description	Default
ufunc	ufunc	The NumPy universal function being applied.	required
method	str	The method of the ufunc ('call', 'reduce', etc.).	required
inputs	tuple	The input arrays to the ufunc.	()
kwargs	dict	Additional keyword arguments to the ufunc.	{}

Returns:

Type	Description
Union[GeoTensor, NDArray]	Union[GeoTensor, NDArray]: If the result is an array, returns a new GeoTensor with the same geospatial attributes. Otherwise, returns the original result.

Source code in georeader/geotensor.py

def __array_ufunc__(self, ufunc, method, *inputs, **kwargs) -> Union["GeoTensor", NDArray]:
    """
    Handle NumPy universal functions applied to this GeoTensor.

    This method is called when a NumPy universal function (ufunc) is applied to the GeoTensor.
    It converts GeoTensor inputs to NumPy arrays, applies the ufunc, and converts array results
    back to GeoTensor objects with the same geospatial metadata.

    Args:
        ufunc (np.ufunc): The NumPy universal function being applied.
        method (str): The method of the ufunc ('__call__', 'reduce', etc.).
        inputs (tuple): The input arrays to the ufunc.
        kwargs (dict): Additional keyword arguments to the ufunc.

    Returns:
        Union[GeoTensor, NDArray]: If the result is an array, returns a new GeoTensor with the same
            geospatial attributes. Otherwise, returns the original result.
    """
    # Normal processing for most operations
    inputs_arr = tuple(
        x.view(np.ndarray) if isinstance(x, GeoTensor) else x for x in inputs
    )
    # Handle 'out' argument if present
    out = kwargs.pop("out", None)
    out_arrays = None

    if out:
        # Convert GeoTensor outputs to regular arrays
        if isinstance(out, tuple):
            out_arrays = tuple(
                o.view(np.ndarray) if isinstance(o, GeoTensor) else o for o in out
            )
        else:
            out_arrays = (
                (out.view(np.ndarray),) if isinstance(out, GeoTensor) else (out,)
            )
        kwargs["out"] = out_arrays

    # Delegate to numpy's implementation
    result = super().__array_ufunc__(ufunc, method, *inputs_arr, **kwargs)

    cast_to_geotensor = self._preserved_spatial(method, **kwargs)

    # Propagate metadata to output arrays
    if out_arrays:
        for o_orig, o_new in zip(
            out if isinstance(out, tuple) else [out], out_arrays
        ):
            if (
                cast_to_geotensor
                and isinstance(o_orig, GeoTensor)
                and isinstance(o_new, np.ndarray)
            ):
                o_new = o_orig.array_as_geotensor(o_new)

    # Normal ufunc processing for other cases
    if cast_to_geotensor:
        return self.array_as_geotensor(result)

    return result

__eq__(other)

Element-wise equality comparison between GeoTensors or with a scalar/array. The georeferencing must match if comparing with another GeoTensor.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	Value to compare with.	required

Raises:

Type	Description
ValueError	If comparing with a GeoTensor and the georeferencing doesn't match.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with boolean values indicating equality.

Source code in georeader/geotensor.py

def __eq__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Element-wise equality comparison between GeoTensors or with a scalar/array.
    The georeferencing must match if comparing with another GeoTensor.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

    Raises:
        ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

    Returns:
        GeoTensor: GeoTensor with boolean values indicating equality.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for comparison. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values == other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

__ge__(other)

Element-wise greater than or equal comparison between GeoTensors or with a scalar/array. The georeferencing must match if comparing with another GeoTensor.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	Value to compare with.	required

Raises:

Type	Description
ValueError	If comparing with a GeoTensor and the georeferencing doesn't match.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with boolean values indicating greater than or equal relationship.

Source code in georeader/geotensor.py

def __ge__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Element-wise greater than or equal comparison between GeoTensors or with a scalar/array.
    The georeferencing must match if comparing with another GeoTensor.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

    Raises:
        ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

    Returns:
        GeoTensor: GeoTensor with boolean values indicating greater than or equal relationship.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for comparison. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values >= other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

__getitem__(key)

Get the values of the GeoTensor using the given key.

Parameters:

Name	Type	Description	Default
key	Any	Key to index the GeoTensor (int, slice, tuple of slices, etc.).	required

Returns:

Name	Type	Description
GeoTensor	GeoTensor	GeoTensor with the selected values.

Source code in georeader/geotensor.py

def __getitem__(self, key: Any) -> "GeoTensor":
    """
    Get the values of the GeoTensor using the given key.

    Args:
        key (Any): Key to index the GeoTensor (int, slice, tuple of slices, etc.).

    Returns:
        GeoTensor: GeoTensor with the selected values.
    """
    if not isinstance(key, tuple):
        key = (key,)

    sel_dict = {}
    for i, k in enumerate(self.dims):
        if i < len(key):
            if key[i] is None:
                raise NotImplementedError(
                    f"Adding axis is not permitted to GeoTensors. Use `expand_dims`"
                )
            elif isinstance(key[i], type(...)):
                raise NotImplementedError(
                    f"Using elipsis is not permitted with GeoTensors. Use `values` attribute"
                )
            sel_dict[k] = key[i]
        else:
            sel_dict[k] = slice(None)

    return self.isel(sel_dict)

__gt__(other)

Element-wise greater than comparison between GeoTensors or with a scalar/array. The georeferencing must match if comparing with another GeoTensor.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	Value to compare with.	required

Raises:

Type	Description
ValueError	If comparing with a GeoTensor and the georeferencing doesn't match.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with boolean values indicating greater than relationship.

Source code in georeader/geotensor.py

def __gt__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Element-wise greater than comparison between GeoTensors or with a scalar/array.
    The georeferencing must match if comparing with another GeoTensor.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

    Raises:
        ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

    Returns:
        GeoTensor: GeoTensor with boolean values indicating greater than relationship.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for comparison. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values > other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

__le__(other)

Element-wise less than or equal comparison between GeoTensors or with a scalar/array. The georeferencing must match if comparing with another GeoTensor.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	Value to compare with.	required

Raises:

Type	Description
ValueError	If comparing with a GeoTensor and the georeferencing doesn't match.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with boolean values indicating less than or equal relationship.

Source code in georeader/geotensor.py

def __le__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Element-wise less than or equal comparison between GeoTensors or with a scalar/array.
    The georeferencing must match if comparing with another GeoTensor.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

    Raises:
        ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

    Returns:
        GeoTensor: GeoTensor with boolean values indicating less than or equal relationship.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for comparison. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values <= other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

__lt__(other)

Element-wise less than comparison between GeoTensors or with a scalar/array. The georeferencing must match if comparing with another GeoTensor.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	Value to compare with.	required

Raises:

Type	Description
ValueError	If comparing with a GeoTensor and the georeferencing doesn't match.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with boolean values indicating less than relationship.

Source code in georeader/geotensor.py

def __lt__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Element-wise less than comparison between GeoTensors or with a scalar/array.
    The georeferencing must match if comparing with another GeoTensor.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

    Raises:
        ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

    Returns:
        GeoTensor: GeoTensor with boolean values indicating less than relationship.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for comparison. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values < other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

__mul__(other)

Multiply this GeoTensor by a value or array element-wise.

Supports broadcasting with scalars, numpy arrays, or other GeoTensors. Common uses include scaling, applying masks, and band math operations.

Parameters:

Name	Type	Description	Default
other	Union[Number, ndarray, GeoTensor]	Value to multiply. Can be: - Scalar (int, float): Multiplies all pixels - np.ndarray: Must be broadcastable with self.values - GeoTensor: Must have same spatial extent	required

Returns:

Name	Type	Description
GeoTensor	Self	New GeoTensor with result of multiplication.

Raises:

Type	Description
ValueError	If other is a GeoTensor and georeferencing doesn't match.

Examples:

>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>>
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>>
>>> # Scale all values (e.g., unit conversion)
>>> gt_scaled = gt * 10000  # Reflectance [0-1] to [0-10000]
>>>
>>> # Apply per-band gain coefficients
>>> gains = np.array([1.1, 1.0, 0.95])[:, None, None]  # Shape: (3, 1, 1)
>>> gt_calibrated = gt * gains
>>>
>>> # Apply binary mask (mask out invalid pixels)
>>> mask = np.random.rand(100, 100) > 0.5  # Boolean mask
>>> gt_masked = gt * mask  # Broadcasts mask to all bands
>>>
>>> # Element-wise product of two rasters
>>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
>>> gt_product = gt * gt2

Note

• For masking, consider using gt[mask] = fill_value for in-place updates
• Result inherits georeferencing from self

See Also

• __truediv__: Division operation
• __setitem__: In-place value assignment with masks

Source code in georeader/geotensor.py

def __mul__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Multiply this GeoTensor by a value or array element-wise.

    Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
    Common uses include scaling, applying masks, and band math operations.

    Args:
        other (Union[numbers.Number, np.ndarray, GeoTensor]): Value to multiply. Can be:
            - Scalar (int, float): Multiplies all pixels
            - np.ndarray: Must be broadcastable with self.values
            - GeoTensor: Must have same spatial extent

    Returns:
        GeoTensor: New GeoTensor with result of multiplication.

    Raises:
        ValueError: If other is a GeoTensor and georeferencing doesn't match.

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>>
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> data = np.random.rand(3, 100, 100)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
        >>>
        >>> # Scale all values (e.g., unit conversion)
        >>> gt_scaled = gt * 10000  # Reflectance [0-1] to [0-10000]
        >>>
        >>> # Apply per-band gain coefficients
        >>> gains = np.array([1.1, 1.0, 0.95])[:, None, None]  # Shape: (3, 1, 1)
        >>> gt_calibrated = gt * gains
        >>>
        >>> # Apply binary mask (mask out invalid pixels)
        >>> mask = np.random.rand(100, 100) > 0.5  # Boolean mask
        >>> gt_masked = gt * mask  # Broadcasts mask to all bands
        >>>
        >>> # Element-wise product of two rasters
        >>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
        >>> gt_product = gt * gt2

    Note:
        - For masking, consider using `gt[mask] = fill_value` for in-place updates
        - Result inherits georeferencing from self

    See Also:
        - `__truediv__`: Division operation
        - `__setitem__`: In-place value assignment with masks
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for multiplication. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values * other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

__ne__(other)

Element-wise inequality comparison between GeoTensors or with a scalar/array. The georeferencing must match if comparing with another GeoTensor.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	Value to compare with.	required

Raises:

Type	Description
ValueError	If comparing with a GeoTensor and the georeferencing doesn't match.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with boolean values indicating inequality.

Source code in georeader/geotensor.py

def __ne__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Element-wise inequality comparison between GeoTensors or with a scalar/array.
    The georeferencing must match if comparing with another GeoTensor.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): Value to compare with.

    Raises:
        ValueError: If comparing with a GeoTensor and the georeferencing doesn't match.

    Returns:
        GeoTensor: GeoTensor with boolean values indicating inequality.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for comparison. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values != other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

__new__(values, transform, crs, fill_value_default=0, attrs=None)

This class is a wrapper around a numpy or torch tensor with geospatial information.

Parameters:

Name	Type	Description	Default
values	NDArray	numpy or torch tensor	required
transform	Affine	affine geospatial transform	required
crs	Any	coordinate reference system	required
fill_value_default	Optional[Union[int, float]]	Value to fill when reading out of bounds. Could be None. Defaults to 0.	0
attrs	Optional[Dict[str, Any]]	dictionary with the attributes of the GeoTensor Defaults to None.	None

Raises:

Type	Description
ValueError	when the shape of the tensor is not 2d, 3d or 4d.

Source code in georeader/geotensor.py

def __new__(
    cls,
    values: NDArray,
    transform: rasterio.Affine,
    crs: Any,
    fill_value_default: Optional[Union[int, float]] = 0,
    attrs: Optional[Dict[str, Any]] = None,
):
    """
    This class is a wrapper around a numpy or torch tensor with geospatial information.

    Args:
        values (NDArray): numpy or torch tensor
        transform (rasterio.Affine): affine geospatial transform
        crs (Any): coordinate reference system
        fill_value_default (Optional[Union[int, float]], optional): Value to fill when
            reading out of bounds. Could be None. Defaults to 0.
        attrs (Optional[Dict[str, Any]], optional): dictionary with the attributes of the GeoTensor
            Defaults to None.

    Raises:
        ValueError: when the shape of the tensor is not 2d, 3d or 4d.
    """
    obj = np.asarray(values).view(cls)

    obj.transform = transform
    obj.crs = crs
    obj.fill_value_default = fill_value_default
    shape = obj.shape
    if (len(shape) < 2) or (len(shape) > 4):
        raise ValueError(f"Expected 2d-4d array found {shape}")

    obj.attrs = attrs if attrs is not None else {}

    return obj

__or__(other)

Perform bitwise OR operation between two GeoTensors. The georeferencing must match.

Parameters:

Name	Type	Description	Default
other	Union[Number, NDArray, GeoTensor]	GeoTensor or array-like to OR with.	required

Raises:

Type	Description
ValueError	if the georeferencing does not match when other is a GeoTensor.

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with the result of the bitwise OR operation.

Source code in georeader/geotensor.py

def __or__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Perform bitwise OR operation between two GeoTensors. The georeferencing must match.

    Args:
        other (Union[numbers.Number, NDArray, GeoTensor]): GeoTensor or array-like to OR with.

    Raises:
        ValueError: if the georeferencing does not match when other is a GeoTensor.

    Returns:
        GeoTensor: GeoTensor with the result of the bitwise OR operation.
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for bitwise OR. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values | other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

__sub__(other)

Subtract a value or array from this GeoTensor element-wise.

Supports broadcasting with scalars, numpy arrays, or other GeoTensors. When subtracting two GeoTensors, they must have the same spatial extent (matching transform, CRS, and spatial dimensions).

Parameters:

Name	Type	Description	Default
other	Union[Number, ndarray, GeoTensor]	Value to subtract. Can be: - Scalar (int, float): Subtracted from all pixels - np.ndarray: Must be broadcastable with self.values - GeoTensor: Must have same spatial extent	required

Returns:

Name	Type	Description
GeoTensor	Self	New GeoTensor with same transform, CRS, and fill_value_default.

Raises:

Type	Description
ValueError	If other is a GeoTensor and georeferencing doesn't match.

Examples:

>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>>
>>> # Create sample GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>>
>>> # Subtract scalar (e.g., remove baseline offset)
>>> gt_adjusted = gt - 0.1
>>>
>>> # Change detection: subtract two GeoTensors
>>> gt_before = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
>>> gt_after = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
>>> change = gt_after - gt_before  # Positive = increase, negative = decrease
>>>
>>> # Center data by subtracting mean per band
>>> band_means = gt.values.mean(axis=(1, 2), keepdims=True)
>>> gt_centered = gt - band_means

See Also

• same_extent: Check if two GeoTensors have matching georeferencing
• read.read_reproject_like: Reproject one GeoTensor to match another

Source code in georeader/geotensor.py

def __sub__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Subtract a value or array from this GeoTensor element-wise.

    Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
    When subtracting two GeoTensors, they must have the same spatial extent
    (matching transform, CRS, and spatial dimensions).

    Args:
        other (Union[numbers.Number, np.ndarray, GeoTensor]): Value to subtract. Can be:
            - Scalar (int, float): Subtracted from all pixels
            - np.ndarray: Must be broadcastable with self.values
            - GeoTensor: Must have same spatial extent

    Returns:
        GeoTensor: New GeoTensor with same transform, CRS, and fill_value_default.

    Raises:
        ValueError: If other is a GeoTensor and georeferencing doesn't match.

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>>
        >>> # Create sample GeoTensor
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> data = np.random.rand(3, 100, 100)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
        >>>
        >>> # Subtract scalar (e.g., remove baseline offset)
        >>> gt_adjusted = gt - 0.1
        >>>
        >>> # Change detection: subtract two GeoTensors
        >>> gt_before = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
        >>> gt_after = GeoTensor(np.random.rand(3, 100, 100), transform, crs="EPSG:32630")
        >>> change = gt_after - gt_before  # Positive = increase, negative = decrease
        >>>
        >>> # Center data by subtracting mean per band
        >>> band_means = gt.values.mean(axis=(1, 2), keepdims=True)
        >>> gt_centered = gt - band_means

    See Also:
        - `same_extent`: Check if two GeoTensors have matching georeferencing
        - `read.read_reproject_like`: Reproject one GeoTensor to match another
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for substraction. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values - other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

__truediv__(other)

Divide this GeoTensor by a value or array element-wise.

Supports broadcasting with scalars, numpy arrays, or other GeoTensors. Common uses include normalization, ratio calculations, and index computation.

Parameters:

Name	Type	Description	Default
other	Union[Number, ndarray, GeoTensor]	Divisor. Can be: - Scalar (int, float): Divides all pixels - np.ndarray: Must be broadcastable with self.values - GeoTensor: Must have same spatial extent	required

Returns:

Name	Type	Description
GeoTensor	Self	New GeoTensor with result of division.

Raises:

Type	Description
ValueError	If other is a GeoTensor and georeferencing doesn't match.

Note

Division by zero produces inf or nan (numpy behavior). Consider adding small epsilon for numerical stability: gt / (other + 1e-10)

Examples:

>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>>
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>>
>>> # Normalize to [0, 1] range
>>> gt_norm = gt / gt.values.max()
>>>
>>> # Per-band normalization
>>> band_maxes = gt.values.max(axis=(1, 2))[:, None, None] + 1e-10
>>> gt_normalized = gt / band_maxes
>>>
>>> # Compute NDVI-like ratio: (NIR - Red) / (NIR + Red)
>>> # Assuming band 0 = Red, band 1 = NIR
>>> red = gt.values[0]
>>> nir = gt.values[1]
>>> ndvi_values = (nir - red) / (nir + red + 1e-10)  # Add epsilon
>>> ndvi = GeoTensor(ndvi_values, gt.transform, gt.crs)
>>>
>>> # Ratio of two rasters
>>> gt2 = GeoTensor(np.random.rand(3, 100, 100) + 0.1, transform, crs="EPSG:32630")
>>> ratio = gt / gt2

See Also

• __mul__: Multiplication operation
• clip: Clip values to valid range after division

Source code in georeader/geotensor.py

def __truediv__(self, other: Union[numbers.Number, NDArray, Self]) -> Self:
    """
    Divide this GeoTensor by a value or array element-wise.

    Supports broadcasting with scalars, numpy arrays, or other GeoTensors.
    Common uses include normalization, ratio calculations, and index computation.

    Args:
        other (Union[numbers.Number, np.ndarray, GeoTensor]): Divisor. Can be:
            - Scalar (int, float): Divides all pixels
            - np.ndarray: Must be broadcastable with self.values
            - GeoTensor: Must have same spatial extent

    Returns:
        GeoTensor: New GeoTensor with result of division.

    Raises:
        ValueError: If other is a GeoTensor and georeferencing doesn't match.

    Note:
        Division by zero produces inf or nan (numpy behavior).
        Consider adding small epsilon for numerical stability: ``gt / (other + 1e-10)``

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>>
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> data = np.random.rand(3, 100, 100)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
        >>>
        >>> # Normalize to [0, 1] range
        >>> gt_norm = gt / gt.values.max()
        >>>
        >>> # Per-band normalization
        >>> band_maxes = gt.values.max(axis=(1, 2))[:, None, None] + 1e-10
        >>> gt_normalized = gt / band_maxes
        >>>
        >>> # Compute NDVI-like ratio: (NIR - Red) / (NIR + Red)
        >>> # Assuming band 0 = Red, band 1 = NIR
        >>> red = gt.values[0]
        >>> nir = gt.values[1]
        >>> ndvi_values = (nir - red) / (nir + red + 1e-10)  # Add epsilon
        >>> ndvi = GeoTensor(ndvi_values, gt.transform, gt.crs)
        >>>
        >>> # Ratio of two rasters
        >>> gt2 = GeoTensor(np.random.rand(3, 100, 100) + 0.1, transform, crs="EPSG:32630")
        >>> ratio = gt / gt2

    See Also:
        - `__mul__`: Multiplication operation
        - `clip`: Clip values to valid range after division
    """
    if isinstance(other, GeoTensor):
        if self.same_extent(other):
            other = other.values
        else:
            raise ValueError(
                "GeoTensor georref must match for division. "
                "Use `read.read_reproject_like(other, self)` to "
                "to reproject `other` to `self` georreferencing."
            )

    result_values = self.values / other

    return GeoTensor(
        result_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

array_as_geotensor(result, fill_value_default=None)

Convert a NumPy array result back to a GeoTensor.

Parameters:

Name	Type	Description	Default
result	Union[ndarray, Self]	Any NumPy array or GeoTensor.	required
fill_value_default	Optional[Number]	fill value for the returned GeoTensor.	None

Returns:

Name	Type	Description
Self	Self	A new GeoTensor with the same geospatial attributes as the original.

Source code in georeader/geotensor.py

def array_as_geotensor(
    self,
    result: Union[np.ndarray, Self],
    fill_value_default: Optional[numbers.Number] = None,
) -> Self:
    """
    Convert a NumPy array result back to a GeoTensor.

    Args:
        result (Union[np.ndarray, Self]): Any NumPy array or GeoTensor.
        fill_value_default: fill value for the returned GeoTensor.

    Returns:
        Self: A new GeoTensor with the same geospatial attributes as the original.
    """

    # Propagate metadata for array results
    if isinstance(result, np.ndarray):
        if result.shape[-2:] != self.shape[-2:]:
            raise ValueError("Operation altered spatial dimensions!")

        if fill_value_default is None:
            fill_value_default = self.fill_value_default

        result = GeoTensor(
            result,
            transform=self.transform,
            crs=self.crs,
            fill_value_default=fill_value_default,
            attrs=self.attrs,
        )

    return result

clip(a_min, a_max)

Clip the GeoTensor values between the GeoTensor min and max values.

Parameters:

Name	Type	Description	Default
a_min	float	Minimum value.	required
a_max	float	Maximum value.	required

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with the clipped values.

Source code in georeader/geotensor.py

def clip(self, a_min: Optional[np.array], a_max: Optional[np.array]) -> Self:
    """
    Clip the GeoTensor values between the GeoTensor min and max values.

    Args:
        a_min (float): Minimum value.
        a_max (float): Maximum value.

    Returns:
        GeoTensor: GeoTensor with the clipped values.
    """
    clipped_values = np.clip(self.values, a_min, a_max)
    return GeoTensor(
        clipped_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

concatenate(geotensors, axis=0) classmethod

Concatenates a list of geotensors along a given axis, assert that all of them has same shape, transform and crs.

Parameters:

Name	Type	Description	Default
geotensors	List[Self]	list of geotensors to concat. All with same shape, transform and crs.	required
axis	int	axis to concatenate. Must be less than the number of dimensions of the geotensors minus 2. default is 0.	0

Returns:

Type	Description
Self	geotensor with extra dim at the front: (len(geotensors),) + shape

Examples:

>>> gt1 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt3 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt = concatenate([gt1, gt2, gt3], axis=0)
>>> assert gt.shape == (9, 100, 100)

Source code in georeader/geotensor.py

@classmethod
def concatenate(cls, geotensors: List[Self], axis: int = 0) -> Self:
    """
    Concatenates a list of geotensors along a given axis, assert that all of them has same shape, transform and crs.

    Args:
        geotensors: list of geotensors to concat. All with same shape, transform and crs.
        axis: axis to concatenate. Must be less than the number of dimensions of the geotensors minus 2.
            default is 0.

    Returns:
        geotensor with extra dim at the front: (len(geotensors),) + shape

    Examples:
        >>> gt1 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
        >>> gt2 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
        >>> gt3 = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
        >>> gt = concatenate([gt1, gt2, gt3], axis=0)
        >>> assert gt.shape == (9, 100, 100)
    """
    assert len(geotensors) > 0, "Empty list provided can't concat"

    if len(geotensors) == 1:
        return geotensors[0].copy()

    first_geotensor = geotensors[0]

    # Assert the axis is NOT an spatial axis
    assert (
        axis < len(first_geotensor.shape) - 2
    ), f"Can't concatenate along spatial axis"

    for i, geo in enumerate(geotensors[1:]):
        assert geo.same_extent(first_geotensor), f"Different extent in concat {i+1}"
        assert (
            geo.shape == first_geotensor.shape
        ), f"Different shape in concat {i+1}"
        assert (
            geo.fill_value_default == first_geotensor.fill_value_default
        ), "Different fill_value_default in concat"

    array_out = np.concatenate([gt.values for gt in geotensors], axis=axis)

    return cls(
        array_out,
        transform=first_geotensor.transform,
        crs=first_geotensor.crs,
        fill_value_default=first_geotensor.fill_value_default,
        attrs=first_geotensor.attrs,
    )

expand_dims(axis)

Expand the dimensions of the GeoTensor values while preserving the spatial dimensions.

This method ensures that no dimensions are added after or in between the spatial dimensions (which are always the last two dimensions).

Parameters:

Name	Type	Description	Default
axis	Union[int, tuple]	Position or positions where new axes should be inserted. Must be less than the number of dimensions minus 2 (to preserve spatial dims). Positions are counted from the first dimension.	required

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with the expanded values.

Raises:

Type	Description
ValueError	If trying to add dimensions at or after the spatial dimensions.

Examples:

>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> # Add a new dimension at axis 0
>>> gt_expanded = gt.expand_dims(0)
>>> assert gt_expanded.shape == (1, 3, 100, 100)

Source code in georeader/geotensor.py

def expand_dims(self, axis: Union[int, tuple]) -> Self:
    """
    Expand the dimensions of the GeoTensor values while preserving the spatial dimensions.

    This method ensures that no dimensions are added after or in between the spatial dimensions
    (which are always the last two dimensions).

    Args:
        axis (Union[int, tuple]): Position or positions where new axes should be inserted.
            Must be less than the number of dimensions minus 2 (to preserve spatial dims).
            Positions are counted from the first dimension.

    Returns:
        GeoTensor: GeoTensor with the expanded values.

    Raises:
        ValueError: If trying to add dimensions at or after the spatial dimensions.

    Examples:
        >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
        >>> # Add a new dimension at axis 0
        >>> gt_expanded = gt.expand_dims(0)
        >>> assert gt_expanded.shape == (1, 3, 100, 100)
    """
    ndim = len(self.shape)

    # Check if axis is valid (not interfering with spatial dimensions)
    if isinstance(axis, int):
        if axis >= ndim - 2 or axis < -ndim:
            raise ValueError(
                f"Cannot add dimension at or after spatial dimensions. "
                f"Axis must be < {ndim - 2} or >= {-ndim}, got {axis}"
            )
        # Convert negative axis to positive
        if axis < 0:
            axis = ndim + axis
    else:  # tuple of axes
        for ax in axis:
            if ax >= ndim - 2 or ax < -ndim:
                raise ValueError(
                    f"Cannot add dimension at or after spatial dimensions. "
                    f"All axes must be < {ndim - 2} or >= {-ndim}, got {ax}"
                )

    # Use numpy expand_dims to add the new dimensions
    expanded_values = np.expand_dims(self.values, axis=axis)

    return GeoTensor(
        expanded_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

footprint(crs=None)

Returns the footprint of the GeoTensor as a Polygon.

Parameters:

Name	Type	Description	Default
crs	Optional[str]	Target coordinate reference system for the footprint. If None, returns footprint in the GeoTensor's native CRS. Common formats: "EPSG:4326" (WGS84), "EPSG:32630" (UTM), CRS object, or WKT string. Use "EPSG:4326" for compatibility with web mapping and GeoJSON. Defaults to None.	None

Returns:

Name	Type	Description
Polygon	Polygon	Shapely Polygon with exactly 5 vertices (4 corners + closing point) representing the rectangular footprint. Coordinates are in the specified CRS (or native CRS if crs=None).

Examples:

>>> import rasterio
>>> from georeader import GeoTensor
>>> import numpy as np
>>>
>>> # Example 1: Get footprint in native CRS
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)  # UTM transform
>>> data = np.random.rand(3, 100, 100)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>> footprint_utm = gt.footprint()
>>> print(f"Bounds (UTM): {footprint_utm.bounds}")
>>> # (xmin, ymin, xmax, ymax) in UTM meters: (500000, 4649000, 501000, 4650000)

>>> # Example 2: Transform footprint to WGS84 for web mapping
>>> footprint_wgs84 = gt.footprint(crs="EPSG:4326")
>>> print(f"Bounds (WGS84): {footprint_wgs84.bounds}")
>>> # (lon_min, lat_min, lon_max, lat_max) in degrees
>>> # Can be used directly in Leaflet, Google Maps, etc.

>>> # Example 3: Check if rasters overlap
>>> gt1 = GeoTensor.load_file('image1.tif')
>>> gt2 = GeoTensor.load_file('image2.tif')
>>> # Get footprints in common CRS
>>> fp1 = gt1.footprint(crs="EPSG:4326")
>>> fp2 = gt2.footprint(crs="EPSG:4326")
>>> if fp1.intersects(fp2):
...     print("Rasters overlap!")
...     overlap_area = fp1.intersection(fp2).area
...     print(f"Overlap area: {overlap_area} square degrees")

>>> # Example 4: Export footprint as GeoJSON
>>> from shapely.geometry import mapping
>>> import json
>>> footprint = gt.footprint(crs="EPSG:4326")
>>> geojson = {
...     "type": "Feature",
...     "geometry": mapping(footprint),
...     "properties": {"name": "Raster extent"}
... }
>>> with open('footprint.geojson', 'w') as f:
...     json.dump(geojson, f)

>>> # Example 5: Check if point is within raster extent
>>> from shapely.geometry import Point
>>> point_of_interest = Point(-3.7038, 40.4168)  # Madrid coordinates
>>> footprint = gt.footprint(crs="EPSG:4326")
>>> if footprint.contains(point_of_interest):
...     print("Point is within raster extent")

>>> # Example 6: Calculate raster area in square kilometers
>>> footprint_utm = gt.footprint()  # In UTM (meters)
>>> area_sqm = footprint_utm.area
>>> area_sqkm = area_sqm / 1_000_000
>>> print(f"Raster covers {area_sqkm:.2f} km²")

Note

• Footprint represents the full raster extent, including nodata regions
• For actual data coverage (excluding nodata), use valid_footprint()
• Polygon always has rectangular shape (4 corners) even for rotated rasters
• CRS transformation is performed if target CRS differs from native CRS
• Footprint is cached-free (computed on-demand from transform and shape)

Source code in georeader/geotensor.py

def footprint(self, crs: Optional[str] = None) -> Polygon:
    """Returns the footprint of the GeoTensor as a Polygon.

    Args:
        crs (Optional[str], optional): Target coordinate reference system for the footprint.
            If None, returns footprint in the GeoTensor's native CRS. Common formats:
            "EPSG:4326" (WGS84), "EPSG:32630" (UTM), CRS object, or WKT string.
            Use "EPSG:4326" for compatibility with web mapping and GeoJSON.
            Defaults to None.

    Returns:
        Polygon: Shapely Polygon with exactly 5 vertices (4 corners + closing point)
            representing the rectangular footprint. Coordinates are in the specified
            CRS (or native CRS if crs=None).

    Examples:
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>> import numpy as np
        >>>
        >>> # Example 1: Get footprint in native CRS
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)  # UTM transform
        >>> data = np.random.rand(3, 100, 100)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
        >>> footprint_utm = gt.footprint()
        >>> print(f"Bounds (UTM): {footprint_utm.bounds}")
        >>> # (xmin, ymin, xmax, ymax) in UTM meters: (500000, 4649000, 501000, 4650000)

        >>> # Example 2: Transform footprint to WGS84 for web mapping
        >>> footprint_wgs84 = gt.footprint(crs="EPSG:4326")
        >>> print(f"Bounds (WGS84): {footprint_wgs84.bounds}")
        >>> # (lon_min, lat_min, lon_max, lat_max) in degrees
        >>> # Can be used directly in Leaflet, Google Maps, etc.

        >>> # Example 3: Check if rasters overlap
        >>> gt1 = GeoTensor.load_file('image1.tif')
        >>> gt2 = GeoTensor.load_file('image2.tif')
        >>> # Get footprints in common CRS
        >>> fp1 = gt1.footprint(crs="EPSG:4326")
        >>> fp2 = gt2.footprint(crs="EPSG:4326")
        >>> if fp1.intersects(fp2):
        ...     print("Rasters overlap!")
        ...     overlap_area = fp1.intersection(fp2).area
        ...     print(f"Overlap area: {overlap_area} square degrees")

        >>> # Example 4: Export footprint as GeoJSON
        >>> from shapely.geometry import mapping
        >>> import json
        >>> footprint = gt.footprint(crs="EPSG:4326")
        >>> geojson = {
        ...     "type": "Feature",
        ...     "geometry": mapping(footprint),
        ...     "properties": {"name": "Raster extent"}
        ... }
        >>> with open('footprint.geojson', 'w') as f:
        ...     json.dump(geojson, f)

        >>> # Example 5: Check if point is within raster extent
        >>> from shapely.geometry import Point
        >>> point_of_interest = Point(-3.7038, 40.4168)  # Madrid coordinates
        >>> footprint = gt.footprint(crs="EPSG:4326")
        >>> if footprint.contains(point_of_interest):
        ...     print("Point is within raster extent")

        >>> # Example 6: Calculate raster area in square kilometers
        >>> footprint_utm = gt.footprint()  # In UTM (meters)
        >>> area_sqm = footprint_utm.area
        >>> area_sqkm = area_sqm / 1_000_000
        >>> print(f"Raster covers {area_sqkm:.2f} km²")

    Note:
        - Footprint represents the full raster extent, including nodata regions
        - For actual data coverage (excluding nodata), use valid_footprint()
        - Polygon always has rectangular shape (4 corners) even for rotated rasters
        - CRS transformation is performed if target CRS differs from native CRS
        - Footprint is cached-free (computed on-demand from transform and shape)
    """
    pol = window_utils.window_polygon(
        rasterio.windows.Window(
            row_off=0, col_off=0, height=self.shape[-2], width=self.shape[-1]
        ),
        self.transform,
    )
    if (crs is None) or window_utils.compare_crs(self.crs, crs):
        return pol

    return window_utils.polygon_to_crs(pol, self.crs, crs)

invalidmask()

Returns a mask of the invalid values of the GeoTensor. The mask is a boolean array with the same shape as the GeoTensor values, where True indicates invalid values and False indicates valid values. The mask is created by comparing the values of the GeoTensor with the self.fill_value_default.

Returns:

Name	Type	Description
Self	Self	GeoTensor with the invalid boolean mask.

Source code in georeader/geotensor.py

def invalidmask(self) -> Self:
    """
    Returns a mask of the invalid values of the GeoTensor. The mask is a boolean array
    with the same shape as the GeoTensor values, where True indicates invalid values and
    False indicates valid values.
    The mask is created by comparing the values of the GeoTensor with the `self.fill_value_default`.

    Returns:
        Self: GeoTensor with the invalid boolean mask.
    """
    if self.fill_value_default is None:
        return GeoTensor(
            np.zeros(self.shape, dtype=bool),
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
        )
    return GeoTensor(
        values=self.values == self.fill_value_default,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

isel(sel)

Select data using dimension-based indexing (xarray-style interface).

Index into the GeoTensor using named dimensions rather than positional indices. Automatically updates the geotransform when slicing spatial dimensions (x, y) to maintain correct georeferencing.

This method provides xarray-compatible syntax while preserving full geospatial metadata. Spatial slices shift the transform origin and optionally rescale resolution if step > 1.

Dimension Mapping::

Standard dimension names for 3D GeoTensor (C, H, W):
┌──────────┬────────────────┬────────────────────┐
│ Dim name │ Array axis      │ Description        │
├──────────┼────────────────┼────────────────────┤
│ "band"   │ axis 0 (C)      │ Spectral bands     │
│ "y"      │ axis 1 (H)      │ Rows (north-south) │
│ "x"      │ axis 2 (W)      │ Cols (east-west)   │
└──────────┴────────────────┴────────────────────┘

Parameters:

Name	Type	Description	Default
sel	Dict[str, Union[slice, list, int]]	Selection dictionary mapping dimension names to slice objects, lists of indices, or integers. For spatial dimensions ("x", "y"), only slice objects are allowed. Common patterns: - {"x": slice(10, 110)}: Columns 10-109 (100 cols) - {"y": slice(20, 120)}: Rows 20-119 (100 rows) - {"band": 0}: Select first band (reduces dimensionality) - {"band": [0, 2, 3]}: Select bands by list - {"x": slice(0, 100, 2)}: Every other column (step=2)	required

Returns:

Name	Type	Description
GeoTensor	Self	Sliced GeoTensor with: - Updated transform: origin shifted to slice start - Updated transform: resolution scaled if step > 1 - Same CRS and fill_value_default

Raises:

Type	Description
NotImplementedError	If dimension name not in self.dims.
NotImplementedError	If spatial dimension uses non-slice indexing.

Examples:

>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>>
>>> # Create 3-band 100x100 GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100).astype(np.float32)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630")
>>> print(f"Original: {gt.shape}, origin: ({gt.transform.c}, {gt.transform.f})")
Original: (3, 100, 100), origin: (500000, 4650000)
>>>
>>> # Slice spatial subset (maintains bands)
>>> subset = gt.isel({"x": slice(20, 80), "y": slice(10, 60)})
>>> print(f"Subset: {subset.shape}")
Subset: (3, 50, 60)
>>> # Transform shifted: new origin at pixel (20, 10)
>>> print(f"New origin: ({subset.transform.c}, {subset.transform.f})")
New origin: (500200.0, 4649900.0)
>>>
>>> # Downsample with step (every 2nd pixel)
>>> downsampled = gt.isel({"x": slice(0, 100, 2), "y": slice(0, 100, 2)})
>>> print(f"Downsampled: {downsampled.shape}")
Downsampled: (3, 50, 50)
>>> # Resolution doubled (10m -> 20m)
>>> print(f"Resolution: {downsampled.res}")
Resolution: (20.0, 20.0)
>>>
>>> # Select single band (reduces to 2D)
>>> red_band = gt.isel({"band": 0})
>>> print(f"Red band: {red_band.shape}")
Red band: (100, 100)
>>>
>>> # Combined: subset + band selection
>>> roi = gt.isel({
...     "band": [1, 2],  # Green and blue
...     "x": slice(25, 75),
...     "y": slice(25, 75)
... })
>>> print(f"ROI: {roi.shape}")
ROI: (2, 50, 50)

See Also

• __getitem__: Positional indexing with tuple syntax: gt[0, 10:60, 20:80]
• read_from_window: Slice using rasterio Window objects
• pad: Add padding to spatial dimensions

Source code in georeader/geotensor.py

def isel(self, sel: Dict[str, Union[slice, list, int]]) -> Self:
    """
    Select data using dimension-based indexing (xarray-style interface).

    Index into the GeoTensor using named dimensions rather than positional
    indices. Automatically updates the geotransform when slicing spatial
    dimensions (x, y) to maintain correct georeferencing.

    This method provides xarray-compatible syntax while preserving full
    geospatial metadata. Spatial slices shift the transform origin and
    optionally rescale resolution if step > 1.

    Dimension Mapping::

        Standard dimension names for 3D GeoTensor (C, H, W):
        ┌──────────┬────────────────┬────────────────────┐
        │ Dim name │ Array axis      │ Description        │
        ├──────────┼────────────────┼────────────────────┤
        │ "band"   │ axis 0 (C)      │ Spectral bands     │
        │ "y"      │ axis 1 (H)      │ Rows (north-south) │
        │ "x"      │ axis 2 (W)      │ Cols (east-west)   │
        └──────────┴────────────────┴────────────────────┘

    Args:
        sel (Dict[str, Union[slice, list, int]]): Selection dictionary mapping
            dimension names to slice objects, lists of indices, or integers.
            For spatial dimensions ("x", "y"), only slice objects are allowed.

            Common patterns:
            - `{"x": slice(10, 110)}`: Columns 10-109 (100 cols)
            - `{"y": slice(20, 120)}`: Rows 20-119 (100 rows)
            - `{"band": 0}`: Select first band (reduces dimensionality)
            - `{"band": [0, 2, 3]}`: Select bands by list
            - `{"x": slice(0, 100, 2)}`: Every other column (step=2)

    Returns:
        GeoTensor: Sliced GeoTensor with:
            - Updated transform: origin shifted to slice start
            - Updated transform: resolution scaled if step > 1
            - Same CRS and fill_value_default

    Raises:
        NotImplementedError: If dimension name not in self.dims.
        NotImplementedError: If spatial dimension uses non-slice indexing.

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>>
        >>> # Create 3-band 100x100 GeoTensor
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> data = np.random.rand(3, 100, 100).astype(np.float32)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630")
        >>> print(f"Original: {gt.shape}, origin: ({gt.transform.c}, {gt.transform.f})")
        Original: (3, 100, 100), origin: (500000, 4650000)
        >>>
        >>> # Slice spatial subset (maintains bands)
        >>> subset = gt.isel({"x": slice(20, 80), "y": slice(10, 60)})
        >>> print(f"Subset: {subset.shape}")
        Subset: (3, 50, 60)
        >>> # Transform shifted: new origin at pixel (20, 10)
        >>> print(f"New origin: ({subset.transform.c}, {subset.transform.f})")
        New origin: (500200.0, 4649900.0)
        >>>
        >>> # Downsample with step (every 2nd pixel)
        >>> downsampled = gt.isel({"x": slice(0, 100, 2), "y": slice(0, 100, 2)})
        >>> print(f"Downsampled: {downsampled.shape}")
        Downsampled: (3, 50, 50)
        >>> # Resolution doubled (10m -> 20m)
        >>> print(f"Resolution: {downsampled.res}")
        Resolution: (20.0, 20.0)
        >>>
        >>> # Select single band (reduces to 2D)
        >>> red_band = gt.isel({"band": 0})
        >>> print(f"Red band: {red_band.shape}")
        Red band: (100, 100)
        >>>
        >>> # Combined: subset + band selection
        >>> roi = gt.isel({
        ...     "band": [1, 2],  # Green and blue
        ...     "x": slice(25, 75),
        ...     "y": slice(25, 75)
        ... })
        >>> print(f"ROI: {roi.shape}")
        ROI: (2, 50, 50)

    See Also:
        - `__getitem__`: Positional indexing with tuple syntax: `gt[0, 10:60, 20:80]`
        - `read_from_window`: Slice using rasterio Window objects
        - `pad`: Add padding to spatial dimensions
    """
    for k in sel:
        if k not in self.dims:
            raise NotImplementedError(f"Axis {k} not in {self.dims}")

    slice_tuple = self._slice_tuple(sel)

    # Compùte the window to shift the transform
    slices_window = []
    for k in ["y", "x"]:
        if k in sel:
            if not isinstance(sel[k], slice):
                raise NotImplementedError(
                    f"Only slice selection supported for x, y dims, found {sel[k]}"
                )
            slices_window.append(sel[k])
        else:
            size = self.width if (k == "x") else self.height
            slices_window.append(slice(0, size))

    window_current = rasterio.windows.Window.from_slices(
        *slices_window, boundless=False, height=self.height, width=self.width
    )

    transform_current = rasterio.windows.transform(
        window_current, transform=self.transform
    )

    # Scale the spatial transform if the step of the slices > 1
    step_rows = slices_window[0].step
    step_cols = slices_window[1].step
    if step_rows is None:
        step_rows = 1

    if step_cols is None:
        step_cols = 1

    if (step_rows != 1) or (step_cols != 1):
        transform_current = transform_current * Affine.scale(step_cols, step_rows)

    return GeoTensor(
        self.values[slice_tuple],
        transform_current,
        self.crs,
        self.fill_value_default,
        attrs=self.attrs,
    )

load_bytes(bytes_read, load_tags=False, rio_env_options=None) classmethod

Load a GeoTensor from a byte stream. It uses rasterio to read the data.

Parameters:

Name	Type	Description	Default
bytes_read	Union[bytes, bytearray, memoryview]	Byte stream to read.	required
load_tags	bool	if True, load the tags of the image. Defaults to False.	False
rio_env_options	Optional[Dict[str, str]]	Rasterio environment options. Defaults to None.	None

Returns:

Name	Type	Description
__class__	Self	GeoTensor object with the loaded data.

Note

The descriptions attribute cannote be loaded from a byte stream. This is a limitation of rasterio. The issue is related to how rasterio.io.MemoryFile handles band descriptions compared to direct file access. This is a known limitation when working with in-memory file representations in GDAL (which rasterio uses under the hood). If you need to load descriptions, you should use georeader.rasterio_reader class.

Source code in georeader/geotensor.py

@classmethod
def load_bytes(
    cls,
    bytes_read: Union[bytes, bytearray, memoryview],
    load_tags: bool = False,
    rio_env_options: Optional[Dict[str, str]] = None,
) -> Self:
    """
    Load a GeoTensor from a byte stream. It uses rasterio to read the data.


    Args:
        bytes_read (Union[bytes, bytearray, memoryview]): Byte stream to read.
        load_tags (bool, optional): if True, load the tags of the image. Defaults to False.
        rio_env_options (Optional[Dict[str, str]], optional): Rasterio environment options. Defaults to None.

    Returns:
        __class__: GeoTensor object with the loaded data.

    Note:
        The `descriptions` attribute cannote be loaded from a byte stream. This is a limitation of `rasterio`.
        The issue is related to how `rasterio.io.MemoryFile` handles band descriptions
        compared to direct file access. This is a known limitation when working
        with in-memory file representations in GDAL (which `rasterio` uses under
        the hood). If you need to load descriptions, you should use `georeader.rasterio_reader`
        class.
    """
    import rasterio.io

    tags = None
    rio_env_options = (
        RIO_ENV_OPTIONS_DEFAULT if rio_env_options is None else rio_env_options
    )
    with rasterio.Env(**rio_env_options):
        with rasterio.io.MemoryFile(bytes_read) as mem:
            with mem.open() as src:
                data = src.read()
                transform = src.transform
                crs = src.crs
                fill_value_default = src.nodata
                if load_tags:
                    tags = src.tags()

    attrs = {}
    if tags is not None:
        attrs["tags"] = tags

    return cls(
        data, transform, crs, fill_value_default=fill_value_default, attrs=attrs
    )

load_file(path, fs=None, load_tags=False, load_descriptions=False, rio_env_options=None) classmethod

Load a GeoTensor from a file. It uses rasterio to read the data. This function loads all the data in memory. For a lazy loading reader use georeader.rasterio_reader.

Parameters:

Name	Type	Description	Default
path	str	Path to the file.	required
fs	Optional[Any]	fsspec.Filesystem object. Defaults to None.	None
load_descriptions	bool	If True, load the description of the image. Defaults to False.	False
load_tags	bool	If True, load the tags of the image. Defaults to False.	False
rio_env_options	Optional[Dict[str, str]]	Rasterio environment options. Defaults to None.	None

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor object with the loaded data.

Source code in georeader/geotensor.py

@classmethod
def load_file(
    cls,
    path: str,
    fs: Optional[Any] = None,
    load_tags: bool = False,
    load_descriptions: bool = False,
    rio_env_options: Optional[Dict[str, str]] = None,
) -> Self:
    """
    Load a GeoTensor from a file. It uses rasterio to read the data. This function
    loads all the data in memory. For a lazy loading reader use `georeader.rasterio_reader`.

    Args:
        path (str): Path to the file.
        fs (Optional[Any], optional): fsspec.Filesystem object. Defaults to None.
        load_descriptions (bool, optional): If True, load the description of the image. Defaults to False.
        load_tags (bool, optional): If True, load the tags of the image. Defaults to False.
        rio_env_options (Optional[Dict[str, str]], optional): Rasterio environment options. Defaults to None.

    Returns:
        GeoTensor: GeoTensor object with the loaded data.
    """

    if fs is not None:
        if load_descriptions:
            raise NotImplementedError(
                """Description loading not supported with `fsspec`. This is because
        the `descriptions` attribute cannote be loaded from a byte stream. This is a limitation of `rasterio`.
        The issue is related to how `rasterio.io.MemoryFile` handles band descriptions
        compared to direct file access. This is a known limitation when working
        with in-memory file representations in GDAL (which `rasterio` uses under
        the hood). If you need to load descriptions, you can use `georeader.rasterio_reader`
        class."""
            )

        with fs.open(path, "rb") as fh:
            return cls.load_bytes(
                fh.read(), load_tags=load_tags, rio_env_options=rio_env_options
            )

    tags = None
    descriptions = None
    rio_env_options = (
        RIO_ENV_OPTIONS_DEFAULT if rio_env_options is None else rio_env_options
    )
    with rasterio.Env(**get_rio_options_path(rio_env_options, path)):
        with rasterio.open(path) as src:
            data = src.read()
            transform = src.transform
            crs = src.crs
            fill_value_default = src.nodata
            if load_tags:
                tags = src.tags()
            if load_descriptions:
                descriptions = tuple(src.descriptions)

    attrs = {}
    if tags is not None:
        attrs["tags"] = tags

    if descriptions is not None:
        attrs["descriptions"] = descriptions

    return cls(
        data, transform, crs, fill_value_default=fill_value_default, attrs=attrs
    )

meshgrid(dst_crs=None)

Create a meshgrid of spatial dimensions of the GeoTensor.

Parameters:

Name	Type	Description	Default
dst_crs	Optional[Any]	output coordinate reference system. Defaults to None.	None

Returns:

Type	Description
Tuple[NDArray, NDArray]	Tuple[NDArray, NDArray]: 2D arrays of xs and ys coordinates.

Source code in georeader/geotensor.py

def meshgrid(self, dst_crs: Optional[Any] = None) -> Tuple[NDArray, NDArray]:
    """
    Create a meshgrid of spatial dimensions of the GeoTensor.

    Args:
        dst_crs (Optional[Any], optional): output coordinate reference system. Defaults to None.

    Returns:
        Tuple[NDArray, NDArray]: 2D arrays of xs and ys coordinates.
    """
    from georeader import griddata

    return griddata.meshgrid(
        self.transform,
        self.width,
        self.height,
        source_crs=self.crs,
        dst_crs=dst_crs,
    )

pad(pad_width, mode='constant', **kwargs)

Add padding around the GeoTensor, updating georeferencing.

Expands the spatial dimensions of the GeoTensor by adding pixels on each side, while automatically updating the geotransform to maintain correct geographic positioning. The new pixels are filled according to the specified mode.

Padding is essential for: - CNN inference requiring fixed input sizes - Avoiding edge artifacts in convolution operations - Creating uniform tile sizes for batch processing

Padding behavior::

pad_width = {"x": (2, 3), "y": (1, 4)}

Original (5×5):           Padded (10×10):
┌─────────────┐           ┌ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ┐
│ █ █ █ █ █ │           │   ← 1 row top           │
│ █ █ █ █ █ │    →      │ ┌─────────────┐          │
│ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
│ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
│ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
└─────────────┘           │ │ █ █ █ █ █ │          │
                          │ │ █ █ █ █ █ │          │
                          │ └─────────────┘          │
                          │   ← 4 rows bottom       │
                          └ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ┘
                            ↑           ↑
                            2 cols      3 cols
                            left        right

Parameters:

Name	Type	Description	Default
pad_width	Union[Dict[str, Tuple[int, int]], List[Tuple[int, int]]]	Padding specification. Can be: Dict: Dimension names to (before, after) padding:: {"x": (left, right), "y": (top, bottom)} List/Tuple: Padding per dimension in order of self.dims:: For 3D GeoTensor with dims (band, y, x): [(band_before, band_after), (y_before, y_after), (x_before, x_after)]	required
mode	str	Padding mode (see numpy.pad). Options: "constant": Fill with constant value (default) "edge": Replicate edge pixels "reflect": Mirror reflection at edges "wrap": Wrap around (toroidal) Defaults to "constant".	'constant'
**kwargs	Any	Additional arguments passed to np.pad: constant_values: Fill value for "constant" mode. Defaults to self.fill_value_default.	{}

Returns:

Name	Type	Description
GeoTensor	GeoTensor	Padded GeoTensor with: - Updated transform: origin shifted by padding amount - Shape: increased by sum of padding in each dimension - Same CRS and fill_value_default

Raises:

Type	Description
ValueError	If list-style pad_width has wrong length.

Examples:

>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>>
>>> # Create 3-band 100x100 GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.random.rand(3, 100, 100).astype(np.float32)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630", fill_value_default=0)
>>> print(f"Original: {gt.shape}")
Original: (3, 100, 100)
>>>
>>> # Pad spatial dimensions (dict style)
>>> padded = gt.pad({"x": (32, 32), "y": (32, 32)})
>>> print(f"Padded: {padded.shape}")
Padded: (3, 164, 164)
>>> # Transform shifted: origin moved by -32 pixels (320m)
>>> print(f"New origin: ({padded.transform.c}, {padded.transform.f})")
New origin: (499680.0, 4650320.0)
>>>
>>> # Pad with list style (must match dims order)
>>> padded = gt.pad([(0, 0), (16, 16), (16, 16)])
>>> print(f"Padded: {padded.shape}")
Padded: (3, 132, 132)
>>>
>>> # Asymmetric padding for edge tiles
>>> edge_padded = gt.pad({"x": (0, 50), "y": (0, 25)})
>>> print(f"Edge padded: {edge_padded.shape}")
Edge padded: (3, 125, 150)
>>>
>>> # Edge replication instead of constant fill
>>> padded_edge = gt.pad({"x": (10, 10), "y": (10, 10)}, mode="edge")

Note

• Transform is updated so padded region has correct georeferencing
• For CNN inference, also consider read_from_window with boundless=True
• Use pad_width={("band": (0, 0), "x": ..., "y": ...} to be explicit

See Also

• pad_array: Lower-level method with dict-only interface
• read_from_window: Read with implicit padding for out-of-bounds
• window_utils.pad_window: Expand window coordinates

Source code in georeader/geotensor.py

def pad(
    self,
    pad_width: Union[Dict[str, Tuple[int, int]], List[Tuple[int, int]]],
    mode: str = "constant",
    **kwargs: Any,
) -> "GeoTensor":
    """
    Add padding around the GeoTensor, updating georeferencing.

    Expands the spatial dimensions of the GeoTensor by adding pixels on each
    side, while automatically updating the geotransform to maintain correct
    geographic positioning. The new pixels are filled according to the
    specified mode.

    Padding is essential for:
    - CNN inference requiring fixed input sizes
    - Avoiding edge artifacts in convolution operations
    - Creating uniform tile sizes for batch processing

    Padding behavior::

        pad_width = {"x": (2, 3), "y": (1, 4)}

        Original (5×5):           Padded (10×10):
        ┌─────────────┐           ┌ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ┐
        │ █ █ █ █ █ │           │   ← 1 row top           │
        │ █ █ █ █ █ │    →      │ ┌─────────────┐          │
        │ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
        │ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
        │ █ █ █ █ █ │           │ │ █ █ █ █ █ │          │
        └─────────────┘           │ │ █ █ █ █ █ │          │
                                  │ │ █ █ █ █ █ │          │
                                  │ └─────────────┘          │
                                  │   ← 4 rows bottom       │
                                  └ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ─ ┘
                                    ↑           ↑
                                    2 cols      3 cols
                                    left        right

    Args:
        pad_width (Union[Dict[str, Tuple[int, int]], List[Tuple[int, int]]]):
            Padding specification. Can be:

            - Dict: Dimension names to (before, after) padding::

                {"x": (left, right), "y": (top, bottom)}
                {"x": (10, 10), "y": (5, 5), "band": (0, 0)}

            - List/Tuple: Padding per dimension in order of self.dims::

                # For 3D GeoTensor with dims (band, y, x):
                [(band_before, band_after), (y_before, y_after), (x_before, x_after)]

        mode (str, optional): Padding mode (see numpy.pad). Options:

            - "constant": Fill with constant value (default)
            - "edge": Replicate edge pixels
            - "reflect": Mirror reflection at edges
            - "wrap": Wrap around (toroidal)

            Defaults to "constant".

        **kwargs: Additional arguments passed to np.pad:

            - constant_values: Fill value for "constant" mode.
              Defaults to self.fill_value_default.

    Returns:
        GeoTensor: Padded GeoTensor with:
            - Updated transform: origin shifted by padding amount
            - Shape: increased by sum of padding in each dimension
            - Same CRS and fill_value_default

    Raises:
        ValueError: If list-style pad_width has wrong length.

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>>
        >>> # Create 3-band 100x100 GeoTensor
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> data = np.random.rand(3, 100, 100).astype(np.float32)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630", fill_value_default=0)
        >>> print(f"Original: {gt.shape}")
        Original: (3, 100, 100)
        >>>
        >>> # Pad spatial dimensions (dict style)
        >>> padded = gt.pad({"x": (32, 32), "y": (32, 32)})
        >>> print(f"Padded: {padded.shape}")
        Padded: (3, 164, 164)
        >>> # Transform shifted: origin moved by -32 pixels (320m)
        >>> print(f"New origin: ({padded.transform.c}, {padded.transform.f})")
        New origin: (499680.0, 4650320.0)
        >>>
        >>> # Pad with list style (must match dims order)
        >>> padded = gt.pad([(0, 0), (16, 16), (16, 16)])
        >>> print(f"Padded: {padded.shape}")
        Padded: (3, 132, 132)
        >>>
        >>> # Asymmetric padding for edge tiles
        >>> edge_padded = gt.pad({"x": (0, 50), "y": (0, 25)})
        >>> print(f"Edge padded: {edge_padded.shape}")
        Edge padded: (3, 125, 150)
        >>>
        >>> # Edge replication instead of constant fill
        >>> padded_edge = gt.pad({"x": (10, 10), "y": (10, 10)}, mode="edge")

    Note:
        - Transform is updated so padded region has correct georeferencing
        - For CNN inference, also consider `read_from_window` with boundless=True
        - Use pad_width={("band": (0, 0), "x": ..., "y": ...} to be explicit

    See Also:
        - `pad_array`: Lower-level method with dict-only interface
        - `read_from_window`: Read with implicit padding for out-of-bounds
        - `window_utils.pad_window`: Expand window coordinates
    """
    if isinstance(pad_width, list) or isinstance(pad_width, tuple):
        if len(pad_width) != len(self.dims):
            raise ValueError(
                f"Expected {len(self.dims)} pad widths found {len(pad_width)}"
            )
        pad_width_dict = {}
        for i, k in enumerate(self.dims):
            pad_width_dict[k] = pad_width[i]
    else:
        pad_width_dict = pad_width
    return self.pad_array(pad_width_dict, mode=mode, **kwargs)

pad_array(pad_width, mode='constant', constant_values=None)

Pad the GeoTensor.

Parameters:

Name	Type	Description	Default
pad_width	_type_	dictionary with Tuple to pad for each dimension {"x": (pad_x_0, pad_x_1), "y": (pad_y_0, pad_y_1)}.	required
mode	str	pad mode (see np.pad or torch.nn.functional.pad). Defaults to "constant".	'constant'
constant_values	Any	description. Defaults to self.fill_value_default.	None

Returns:

Name	Type	Description
GeoTensor	Self	padded GeoTensor.

Examples:

>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> gt.pad_array({"x": (10, 10), "y": (10, 10)})
>>> assert gt.shape == (3, 120, 120)

Source code in georeader/geotensor.py

def pad_array(
    self,
    pad_width: Dict[str, Tuple[int, int]],
    mode: str = "constant",
    constant_values: Optional[Any] = None,
) -> Self:
    """
    Pad the GeoTensor.

    Args:
        pad_width (_type_, optional):  dictionary with Tuple to pad for each dimension
            `{"x": (pad_x_0, pad_x_1), "y": (pad_y_0, pad_y_1)}`.
        mode (str, optional): pad mode (see np.pad or torch.nn.functional.pad). Defaults to "constant".
        constant_values (Any, optional): _description_. Defaults to `self.fill_value_default`.

    Returns:
        GeoTensor: padded GeoTensor.

    Examples:
        >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
        >>> gt.pad_array({"x": (10, 10), "y": (10, 10)})
        >>> assert gt.shape == (3, 120, 120)
    """
    if constant_values is None and mode == "constant":
        if self.fill_value_default is None:
            raise ValueError(
                f"Mode constant either requires constant_values passed or fill_value_default not None in current GeoTensor"
            )
        constant_values = self.fill_value_default

    pad_list_np = []
    for k in self.dims:
        if k in pad_width:
            pad_list_np.append(pad_width[k])
        else:
            pad_list_np.append((0, 0))

    kwargs_extra = {}
    if mode == "constant":
        kwargs_extra["constant_values"] = constant_values
    values_new = np.pad(self.values, tuple(pad_list_np), mode=mode, **kwargs_extra)

    # Compute the new transform
    slices_window = []
    for k in ["y", "x"]:
        size = self.width if (k == "x") else self.height
        if k in pad_width:
            slices_window.append(slice(-pad_width[k][0], size + pad_width[k][1]))
        else:
            slices_window.append(slice(0, size))

    window_current = rasterio.windows.Window.from_slices(
        *slices_window, boundless=True
    )
    transform_current = rasterio.windows.transform(
        window_current, transform=self.transform
    )
    return GeoTensor(
        values_new,
        transform_current,
        self.crs,
        self.fill_value_default,
        attrs=self.attrs,
    )

read_from_window(window, boundless=True)

Extract a spatial subset using a rasterio Window.

Reads data from the GeoTensor corresponding to the specified pixel window, with optional boundless reading that pads out-of-bounds regions with the fill value. This method provides the same interface as rasterio's windowed reading, enabling seamless transitions between lazy file readers and in-memory GeoTensors.

Boundless vs Bounded Reading::

Window extends beyond data:     boundless=True:      boundless=False:
┌───────────────┐                 ┌─────────┐           ┌─────┐
│ ┌─────────┐   │                 │▒▒▒▒▒▒▒▒▒│           │     │
│ │  DATA   │   │                 │▒▒▒███▒▒▒│           │ ███ │
│ │         │   │  ─────────►    │▒▒▒███▒▒▒│           │ ███ │
│ └─────────┘   │                 │▒▒▒▒▒▒▒▒▒│           └─────┘
└───────────────┘                 └─────────┘           Returns only
     Window                   Padded with            intersection
     request                  fill_value

Parameters:

Name	Type	Description	Default
window	Window	Pixel window defining the region to extract. Format: Window(col_off, row_off, width, height). May have negative offsets or extend beyond data bounds if boundless=True.	required
boundless	bool	If True, allows reading windows that extend beyond data boundaries. Out-of-bounds regions are filled with self.fill_value_default. If False, returns only the intersection with the data extent. Defaults to True.	True

Returns:

Name	Type	Description
GeoTensor	Self	Extracted data with: - Shape: (C, window.height, window.width) if boundless=True, or smaller if boundless=False and window exceeds bounds - Transform: Updated to window origin - Same CRS and fill_value_default

Raises:

Type	Description
WindowError	If window does not intersect the data at all (no overlap).

Examples:

>>> import numpy as np
>>> import rasterio
>>> import rasterio.windows
>>> from georeader import GeoTensor
>>>
>>> # Create test GeoTensor
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> data = np.arange(300).reshape(3, 10, 10).astype(np.float32)
>>> gt = GeoTensor(data, transform, crs="EPSG:32630", fill_value_default=-1)
>>> print(f"Data shape: {gt.shape}")
Data shape: (3, 10, 10)
>>>
>>> # Read window fully within data
>>> window = rasterio.windows.Window(2, 3, 5, 4)  # col_off, row_off, width, height
>>> subset = gt.read_from_window(window)
>>> print(f"Subset shape: {subset.shape}")
Subset shape: (3, 4, 5)
>>>
>>> # Boundless read: window extends beyond data
>>> window_edge = rasterio.windows.Window(7, 7, 6, 6)  # Extends 3 pixels past edge
>>> padded = gt.read_from_window(window_edge, boundless=True)
>>> print(f"Padded shape: {padded.shape}")
Padded shape: (3, 6, 6)
>>> # Corners filled with fill_value_default (-1)
>>> print(f"Has padding: {(padded.values == -1).any()}")
Has padding: True
>>>
>>> # Bounded read: only get intersection
>>> clipped = gt.read_from_window(window_edge, boundless=False)
>>> print(f"Clipped shape: {clipped.shape}")
Clipped shape: (3, 3, 3)
>>>
>>> # CNN tiling: process image in overlapping windows
>>> tile_size = 4
>>> stride = 2
>>> for row in range(0, gt.height, stride):
...     for col in range(0, gt.width, stride):
...         tile_window = rasterio.windows.Window(col, row, tile_size, tile_size)
...         tile = gt.read_from_window(tile_window, boundless=True)
...         # Process tile with CNN...

Note

• Boundless reading is essential for processing edge tiles in CNNs
• Transform is always updated to match the output data origin
• Use get_slice_pad directly for manual slice/pad control
• Compatible with RasterioReader interface for polymorphic code

See Also

• isel: Dictionary-based dimension indexing
• window_utils.get_slice_pad: Compute slice and padding parameters
• RasterioReader.read_from_window: Lazy file-based equivalent

Source code in georeader/geotensor.py

def read_from_window(
    self, window: rasterio.windows.Window, boundless: bool = True
) -> Self:
    """
    Extract a spatial subset using a rasterio Window.

    Reads data from the GeoTensor corresponding to the specified pixel window,
    with optional boundless reading that pads out-of-bounds regions with the
    fill value. This method provides the same interface as rasterio's windowed
    reading, enabling seamless transitions between lazy file readers and
    in-memory GeoTensors.

    Boundless vs Bounded Reading::

        Window extends beyond data:     boundless=True:      boundless=False:
        ┌───────────────┐                 ┌─────────┐           ┌─────┐
        │ ┌─────────┐   │                 │▒▒▒▒▒▒▒▒▒│           │     │
        │ │  DATA   │   │                 │▒▒▒███▒▒▒│           │ ███ │
        │ │         │   │  ─────────►    │▒▒▒███▒▒▒│           │ ███ │
        │ └─────────┘   │                 │▒▒▒▒▒▒▒▒▒│           └─────┘
        └───────────────┘                 └─────────┘           Returns only
             Window                   Padded with            intersection
             request                  fill_value

    Args:
        window (rasterio.windows.Window): Pixel window defining the region
            to extract. Format: Window(col_off, row_off, width, height).
            May have negative offsets or extend beyond data bounds if
            boundless=True.
        boundless (bool, optional): If True, allows reading windows that
            extend beyond data boundaries. Out-of-bounds regions are filled
            with `self.fill_value_default`. If False, returns only the
            intersection with the data extent. Defaults to True.

    Returns:
        GeoTensor: Extracted data with:
            - Shape: (C, window.height, window.width) if boundless=True,
              or smaller if boundless=False and window exceeds bounds
            - Transform: Updated to window origin
            - Same CRS and fill_value_default

    Raises:
        rasterio.windows.WindowError: If window does not intersect the data
            at all (no overlap).

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> import rasterio.windows
        >>> from georeader import GeoTensor
        >>>
        >>> # Create test GeoTensor
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> data = np.arange(300).reshape(3, 10, 10).astype(np.float32)
        >>> gt = GeoTensor(data, transform, crs="EPSG:32630", fill_value_default=-1)
        >>> print(f"Data shape: {gt.shape}")
        Data shape: (3, 10, 10)
        >>>
        >>> # Read window fully within data
        >>> window = rasterio.windows.Window(2, 3, 5, 4)  # col_off, row_off, width, height
        >>> subset = gt.read_from_window(window)
        >>> print(f"Subset shape: {subset.shape}")
        Subset shape: (3, 4, 5)
        >>>
        >>> # Boundless read: window extends beyond data
        >>> window_edge = rasterio.windows.Window(7, 7, 6, 6)  # Extends 3 pixels past edge
        >>> padded = gt.read_from_window(window_edge, boundless=True)
        >>> print(f"Padded shape: {padded.shape}")
        Padded shape: (3, 6, 6)
        >>> # Corners filled with fill_value_default (-1)
        >>> print(f"Has padding: {(padded.values == -1).any()}")
        Has padding: True
        >>>
        >>> # Bounded read: only get intersection
        >>> clipped = gt.read_from_window(window_edge, boundless=False)
        >>> print(f"Clipped shape: {clipped.shape}")
        Clipped shape: (3, 3, 3)
        >>>
        >>> # CNN tiling: process image in overlapping windows
        >>> tile_size = 4
        >>> stride = 2
        >>> for row in range(0, gt.height, stride):
        ...     for col in range(0, gt.width, stride):
        ...         tile_window = rasterio.windows.Window(col, row, tile_size, tile_size)
        ...         tile = gt.read_from_window(tile_window, boundless=True)
        ...         # Process tile with CNN...

    Note:
        - Boundless reading is essential for processing edge tiles in CNNs
        - Transform is always updated to match the output data origin
        - Use `get_slice_pad` directly for manual slice/pad control
        - Compatible with RasterioReader interface for polymorphic code

    See Also:
        - `isel`: Dictionary-based dimension indexing
        - `window_utils.get_slice_pad`: Compute slice and padding parameters
        - `RasterioReader.read_from_window`: Lazy file-based equivalent
    """

    window_data = rasterio.windows.Window(
        col_off=0, row_off=0, width=self.width, height=self.height
    )
    if boundless:
        slice_dict, pad_width = window_utils.get_slice_pad(window_data, window)
        need_pad = any(p != 0 for p in pad_width["x"] + pad_width["y"])
        X_sliced = self.isel(slice_dict)
        if need_pad:
            X_sliced = X_sliced.pad(
                pad_width=pad_width,
                mode="constant",
                constant_values=self.fill_value_default,
            )
        return X_sliced
    else:
        window_read = rasterio.windows.intersection(window, window_data)
        slice_y, slice_x = window_read.toslices()
        slice_dict = {"x": slice_x, "y": slice_y}
        slices_ = self._slice_tuple(slice_dict)
        transform_current = rasterio.windows.transform(
            window_read, transform=self.transform
        )
        return GeoTensor(
            self.values[slices_],
            transform_current,
            self.crs,
            self.fill_value_default,
            attrs=self.attrs,
        )

resize(output_shape=None, resolution_dst=None, anti_aliasing=True, anti_aliasing_sigma=None, interpolation='bilinear', mode_pad='constant')

Resize the geotensor to match a certain size output_shape. This function works with GeoTensors of 2D, 3D and 4D. The geoinformation of the output tensor is changed accordingly.

Parameters:

Name	Type	Description	Default
output_shape	Optional[Tuple[int, int]]	output spatial shape if None resolution_dst must be provided. If not provided, the output shape is computed from the resolution_dst rounding to the closest integer.	None
resolution_dst	Optional[Tuple[float, float]]	output resolution if None output_shape must be provided.	None
anti_aliasing	bool	Whether to apply a Gaussian filter to smooth the image prior to downsampling	True
anti_aliasing_sigma	Optional[Union[float, ndarray]]	anti_aliasing_sigma : {float}, optional Standard deviation for Gaussian filtering used when anti-aliasing. By default, this value is chosen as (s - 1) / 2 where s is the downsampling factor, where s > 1	None
interpolation	Optional[str]	Algorithm used for resizing: 'nearest' | 'bilinear' | 'bicubic'	'bilinear'
mode_pad	str	mode pad for resize function	'constant'

Returns:

Type	Description
Self	resized GeoTensor

Examples:

>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> resized = gt.resize((50, 50))
>>> assert resized.shape == (3, 50, 50)
>>> assert resized.res == (2*gt.res[0], 2*gt.res[1])

Source code in georeader/geotensor.py

def resize(
    self,
    output_shape: Optional[Tuple[int, int]] = None,
    resolution_dst: Optional[Tuple[float, float]] = None,
    anti_aliasing: bool = True,
    anti_aliasing_sigma: Optional[Union[float, np.ndarray]] = None,
    interpolation: Optional[str] = "bilinear",
    mode_pad: str = "constant",
) -> Self:
    """
    Resize the geotensor to match a certain size output_shape. This function works with GeoTensors of 2D, 3D and 4D.
    The geoinformation of the output tensor is changed accordingly.

    Args:
        output_shape: output spatial shape if None resolution_dst must be provided. If not provided,
            the output shape is computed from the resolution_dst rounding to the closest integer.
        resolution_dst: output resolution if None output_shape must be provided.
        anti_aliasing: Whether to apply a Gaussian filter to smooth the image prior to downsampling
        anti_aliasing_sigma:  anti_aliasing_sigma : {float}, optional
            Standard deviation for Gaussian filtering used when anti-aliasing.
            By default, this value is chosen as (s - 1) / 2 where s is the
            downsampling factor, where s > 1
        interpolation: Algorithm used for resizing: 'nearest' | 'bilinear' | 'bicubic'
        mode_pad: mode pad for resize function

    Returns:
         resized GeoTensor

    Examples:
        >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
        >>> resized = gt.resize((50, 50))
        >>> assert resized.shape == (3, 50, 50)
        >>> assert resized.res == (2*gt.res[0], 2*gt.res[1])
    """
    input_shape = self.shape
    spatial_shape = input_shape[-2:]
    resolution_or = self.res

    if output_shape is None:
        assert (
            resolution_dst is not None
        ), f"Can't have output_shape and resolution_dst as None"
        output_shape = int(
            round(spatial_shape[0] * resolution_or[0] / resolution_dst[0])
        ), int(round(spatial_shape[1] * resolution_or[1] / resolution_dst[1]))
    else:
        assert (
            resolution_dst is None
        ), f"Both output_shape and resolution_dst can't be provided"
        assert (
            len(output_shape) == 2
        ), f"Expected output shape to be the spatial dimensions found: {output_shape}"
        resolution_dst = (
            spatial_shape[0] * resolution_or[0] / output_shape[0],
            spatial_shape[1] * resolution_or[1] / output_shape[1],
        )

    # Compute output transform
    transform_scale = rasterio.Affine.scale(
        resolution_dst[0] / resolution_or[0], resolution_dst[1] / resolution_or[1]
    )
    transform = self.transform * transform_scale

    from skimage.transform import resize

    # https://scikit-image.org/docs/stable/api/skimage.transform.html#skimage.transform.resize
    output_tensor = np.ndarray(input_shape[:-2] + output_shape, dtype=self.dtype)
    if len(input_shape) == 4:
        for i, j in product(range(0, input_shape[0]), range(0, input_shape[1])):
            if (
                (not anti_aliasing)
                or (anti_aliasing_sigma is None)
                or isinstance(anti_aliasing_sigma, numbers.Number)
            ):
                anti_aliasing_sigma_iter = anti_aliasing_sigma
            else:
                anti_aliasing_sigma_iter = anti_aliasing_sigma[i, j]
            output_tensor[i, j] = resize(
                self.values[i, j],
                output_shape,
                order=ORDERS[interpolation],
                anti_aliasing=anti_aliasing,
                preserve_range=False,
                cval=self.fill_value_default,
                mode=mode_pad,
                anti_aliasing_sigma=anti_aliasing_sigma_iter,
            )
    elif len(input_shape) == 3:
        for i in range(0, input_shape[0]):
            if (
                (not anti_aliasing)
                or (anti_aliasing_sigma is None)
                or isinstance(anti_aliasing_sigma, numbers.Number)
            ):
                anti_aliasing_sigma_iter = anti_aliasing_sigma
            else:
                anti_aliasing_sigma_iter = anti_aliasing_sigma[i]
            output_tensor[i] = resize(
                self.values[i],
                output_shape,
                order=ORDERS[interpolation],
                anti_aliasing=anti_aliasing,
                preserve_range=False,
                cval=self.fill_value_default,
                mode=mode_pad,
                anti_aliasing_sigma=anti_aliasing_sigma_iter,
            )
    else:
        output_tensor[...] = resize(
            self.values,
            output_shape,
            order=ORDERS[interpolation],
            anti_aliasing=anti_aliasing,
            preserve_range=False,
            cval=self.fill_value_default,
            mode=mode_pad,
            anti_aliasing_sigma=anti_aliasing_sigma,
        )

    return GeoTensor(
        output_tensor,
        transform=transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

same_extent(other, precision=0.001)

Check if two GeoTensors have the same georeferencing (crs, transform and spatial dimensions).

Parameters:

Name	Type	Description	Default
other	GeoTensor | GeoData	GeoTensor to compare with. Other GeoData object can be passed (it requires crs, transform and shape attributes)	required
precision	float	precision to compare the transform. Defaults to 1e-3.	0.001

Returns:

Name	Type	Description
bool	bool	True if both GeoTensors have the same georeferencing.

Source code in georeader/geotensor.py

def same_extent(self, other: Self, precision: float = 1e-3) -> bool:
    """
    Check if two GeoTensors have the same georeferencing (crs, transform and spatial dimensions).

    Args:
        other (GeoTensor | GeoData): GeoTensor to compare with. Other GeoData object can be passed (it requires crs, transform and shape attributes)
        precision (float, optional): precision to compare the transform. Defaults to 1e-3.

    Returns:
        bool: True if both GeoTensors have the same georeferencing.
    """
    return (
        self.transform.almost_equals(other.transform, precision=precision)
        and window_utils.compare_crs(self.crs, other.crs)
        and (self.shape[-2:] == other.shape[-2:])
    )

set_dtype(dtype)

Set the dtype of the GeoTensor in-place.

.. deprecated:: Use astype() instead for a cleaner numpy-compatible interface. This method modifies the array in place which can be confusing.

Parameters:

Name	Type	Description	Default
dtype	dtype	Target data type for the array.	required

Raises:

Type	Description
RuntimeError	If the dtype cannot be changed in-place (due to numpy subclass limitations).

Source code in georeader/geotensor.py

def set_dtype(self, dtype: np.dtype) -> None:
    """Set the dtype of the GeoTensor in-place.

    .. deprecated::
        Use astype() instead for a cleaner numpy-compatible interface.
        This method modifies the array in place which can be confusing.

    Args:
        dtype (np.dtype): Target data type for the array.

    Raises:
        RuntimeError: If the dtype cannot be changed in-place (due to numpy subclass limitations).
    """
    warnings.warn(
        "set_dtype() is deprecated. Use astype() for a cleaner interface.",
        DeprecationWarning,
        stacklevel=2
    )
    original_dtype = self.dtype
    requested_dtype = np.dtype(dtype)

    # If requesting the same dtype, nothing to do
    if original_dtype == requested_dtype:
        return

    # Convert using numpy's astype
    result = np.asarray(self).astype(dtype)
    # Since we can't truly change dtype in-place when sizes differ,
    # we update the internal data by using the ndarray resize mechanism
    # This modifies the underlying buffer
    np.ndarray.resize(self, result.shape, refcheck=False)
    self[:] = result

    # Check if dtype actually changed when a change was requested
    if self.dtype != requested_dtype:
        raise RuntimeError(
            f"set_dtype() failed to change dtype from {original_dtype} to {requested_dtype}. "
            f"This is due to numpy subclass limitations. Use astype() instead to create "
            f"a new GeoTensor with the desired dtype."
        )

squeeze(axis=None)

Remove single-dimensional entries from the shape of the GeoTensor values. It does not squeeze the spatial dimensions (last two dimensions).

Returns:

Name	Type	Description
GeoTensor	Self	GeoTensor with the squeezed values.

Source code in georeader/geotensor.py

def squeeze(self, axis=None) -> Self:
    """
    Remove single-dimensional entries from the shape of the GeoTensor values.
    It does not squeeze the spatial dimensions (last two dimensions).

    Returns:
        GeoTensor: GeoTensor with the squeezed values.
    """
    if axis is None:
        axis = tuple(range(self.values.ndim - 2))
    else:
        if isinstance(axis, int):
            axis = (axis,)
        # Check if spatial dimesions will be squeezed
        if self.width == 1:
            if any(a in (-1, self.values.ndim - 1) for a in axis):
                raise ValueError(
                    "Cannot squeeze spatial dimensions. "
                    "Use `squeeze(axis=0)` to squeeze the first dimension."
                )
        elif self.height == 1:
            if any(a in (-2, self.values.ndim - 2) for a in axis):
                raise ValueError(
                    "Cannot squeeze spatial dimensions. "
                    "Use `squeeze(axis=0)` to squeeze the first dimension."
                )

    # squeeze all but last two dimensions
    squeezed_values = np.squeeze(self.values, axis=axis)

    return GeoTensor(
        squeezed_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

stack(geotensors) classmethod

Stack multiple GeoTensors along a new leading dimension.

Creates a new GeoTensor with an additional dimension at the front, containing all input GeoTensors. All inputs must have identical georeferencing (transform, CRS) and shape.

This is useful for: - Creating time-series stacks from multiple acquisitions - Batching for CNN inference - Multi-temporal analysis

Stacking behavior::

Input: 3 GeoTensors each (3, H, W)

gt1: (3, 100, 100)  ─┐
gt2: (3, 100, 100)  ─┼──► stacked: (3, 3, 100, 100)
gt3: (3, 100, 100)  ─┘     ↑
                        new time/batch dim

Parameters:

Name	Type	Description	Default
geotensors	List[GeoTensor]	List of GeoTensors to stack. Requirements: Non-empty list All must have identical shape All must have same transform (same_extent) All must have same CRS All must have same fill_value_default	required

Returns:

Name	Type	Description
GeoTensor	Self	Stacked GeoTensor with: Shape: (len(geotensors),) + original_shape Transform: Same as input GeoTensors CRS: Same as input GeoTensors fill_value_default: Same as input GeoTensors attrs: Copied from first GeoTensor

Raises:

Type	Description
AssertionError	If list is empty.
AssertionError	If GeoTensors have different extents.
AssertionError	If GeoTensors have different shapes.
AssertionError	If GeoTensors have different fill_value_default.

Examples:

>>> import numpy as np
>>> import rasterio
>>> from georeader import GeoTensor
>>>
>>> # Create multiple GeoTensors (e.g., different dates)
>>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
>>> gt_jan = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
>>> gt_feb = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
>>> gt_mar = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
>>>
>>> # Stack into time series
>>> time_series = GeoTensor.stack([gt_jan, gt_feb, gt_mar])
>>> print(f"Time series shape: {time_series.shape}")
Time series shape: (3, 3, 100, 100)
>>> # dims: (time, band, y, x)
>>>
>>> # Access individual time steps
>>> jan_data = time_series[0]  # Returns GeoTensor with shape (3, 100, 100)
>>>
>>> # Compute temporal median (across time dimension)
>>> median = np.nanmedian(time_series.values, axis=0)
>>> # median shape: (3, 100, 100)
>>>
>>> # Stack single GeoTensor (adds dimension)
>>> single_stacked = GeoTensor.stack([gt_jan])
>>> print(f"Single stacked: {single_stacked.shape}")
Single stacked: (1, 3, 100, 100)

Note

• All GeoTensors must have identical georeferencing
• Use concatenate to join along existing dimension instead
• For memory efficiency with large stacks, consider lazy loading
• The first GeoTensor's attrs are preserved

See Also

• concatenate: Join along existing axis (no new dimension)
• same_extent: Check if GeoTensors have matching georeferencing
• numpy.stack: Underlying operation for values

Source code in georeader/geotensor.py

@classmethod
def stack(cls, geotensors: List[Self]) -> Self:
    """
    Stack multiple GeoTensors along a new leading dimension.

    Creates a new GeoTensor with an additional dimension at the front,
    containing all input GeoTensors. All inputs must have identical
    georeferencing (transform, CRS) and shape.

    This is useful for:
    - Creating time-series stacks from multiple acquisitions
    - Batching for CNN inference
    - Multi-temporal analysis

    Stacking behavior::

        Input: 3 GeoTensors each (3, H, W)

        gt1: (3, 100, 100)  ─┐
        gt2: (3, 100, 100)  ─┼──► stacked: (3, 3, 100, 100)
        gt3: (3, 100, 100)  ─┘     ↑
                                new time/batch dim

    Args:
        geotensors (List[GeoTensor]): List of GeoTensors to stack. Requirements:

            - Non-empty list
            - All must have identical shape
            - All must have same transform (same_extent)
            - All must have same CRS
            - All must have same fill_value_default

    Returns:
        GeoTensor: Stacked GeoTensor with:

            - Shape: (len(geotensors),) + original_shape
            - Transform: Same as input GeoTensors
            - CRS: Same as input GeoTensors
            - fill_value_default: Same as input GeoTensors
            - attrs: Copied from first GeoTensor

    Raises:
        AssertionError: If list is empty.
        AssertionError: If GeoTensors have different extents.
        AssertionError: If GeoTensors have different shapes.
        AssertionError: If GeoTensors have different fill_value_default.

    Examples:
        >>> import numpy as np
        >>> import rasterio
        >>> from georeader import GeoTensor
        >>>
        >>> # Create multiple GeoTensors (e.g., different dates)
        >>> transform = rasterio.Affine(10, 0, 500000, 0, -10, 4650000)
        >>> gt_jan = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
        >>> gt_feb = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
        >>> gt_mar = GeoTensor(np.random.rand(3, 100, 100), transform, "EPSG:32630")
        >>>
        >>> # Stack into time series
        >>> time_series = GeoTensor.stack([gt_jan, gt_feb, gt_mar])
        >>> print(f"Time series shape: {time_series.shape}")
        Time series shape: (3, 3, 100, 100)
        >>> # dims: (time, band, y, x)
        >>>
        >>> # Access individual time steps
        >>> jan_data = time_series[0]  # Returns GeoTensor with shape (3, 100, 100)
        >>>
        >>> # Compute temporal median (across time dimension)
        >>> median = np.nanmedian(time_series.values, axis=0)
        >>> # median shape: (3, 100, 100)
        >>>
        >>> # Stack single GeoTensor (adds dimension)
        >>> single_stacked = GeoTensor.stack([gt_jan])
        >>> print(f"Single stacked: {single_stacked.shape}")
        Single stacked: (1, 3, 100, 100)

    Note:
        - All GeoTensors must have identical georeferencing
        - Use `concatenate` to join along existing dimension instead
        - For memory efficiency with large stacks, consider lazy loading
        - The first GeoTensor's attrs are preserved

    See Also:
        - `concatenate`: Join along existing axis (no new dimension)
        - `same_extent`: Check if GeoTensors have matching georeferencing
        - `numpy.stack`: Underlying operation for values
    """
    assert len(geotensors) > 0, "Empty list provided can't concat"

    if len(geotensors) == 1:
        gt = geotensors[0]
        return GeoTensor(
            gt.values[np.newaxis],
            transform=gt.transform,
            crs=gt.crs,
            fill_value_default=gt.fill_value_default,
            attrs=gt.attrs,
        )

    first_geotensor = geotensors[0]
    array_out = np.zeros(
        (len(geotensors),) + first_geotensor.shape, dtype=first_geotensor.dtype
    )
    array_out[0] = first_geotensor.values

    for i, geo in enumerate(geotensors[1:]):
        assert geo.same_extent(first_geotensor), f"Different size in concat {i+1}"
        assert (
            geo.shape == first_geotensor.shape
        ), f"Different shape in concat {i+1}"
        assert (
            geo.fill_value_default == first_geotensor.fill_value_default
        ), "Different fill_value_default in concat"
        array_out[i + 1] = geo.values

    return cls(
        array_out,
        transform=first_geotensor.transform,
        crs=first_geotensor.crs,
        fill_value_default=first_geotensor.fill_value_default,
        attrs=first_geotensor.attrs,
    )

transpose(axes=None)

Permute the dimensions of the GeoTensor while keeping the spatial dimensions at the end.

Parameters:

Name	Type	Description	Default
axes	tuple	If specified, it must be a tuple or list of axes. The last two values must be the original spatial dimensions indices (ndim-2, ndim-1). If None, the non-spatial dimensions are reversed while spatial dimensions remain at the end.	None

Returns:

Name	Type	Description
GeoTensor	Self	A view of the array with dimensions transposed.

Raises:

Type	Description
ValueError	If the spatial dimensions are moved from their last positions.

Examples:

>>> gt = GeoTensor(np.random.rand(3, 4, 100, 200), transform, crs)
>>> # Shape is (3, 4, 100, 200)
>>> gt_t = gt.transpose()
>>> # Shape is now (4, 3, 100, 200)
>>>
>>> # You can also specify axes explicitly:
>>> gt_t = gt.transpose((1, 0, 2, 3))  # Valid: spatial dims remain at end
>>> # But this would raise an error:
>>> # gt.transpose((0, 2, 1, 3))  # Invalid: spatial dims must stay at end

Source code in georeader/geotensor.py

def transpose(self, axes=None) -> Self:
    """
    Permute the dimensions of the GeoTensor while keeping the spatial dimensions at the end.

    Args:
        axes (tuple, optional): If specified, it must be a tuple or list of axes. The last two
            values must be the original spatial dimensions indices (ndim-2, ndim-1).
            If None, the non-spatial dimensions are reversed while spatial dimensions remain at the end.

    Returns:
        GeoTensor: A view of the array with dimensions transposed.

    Raises:
        ValueError: If the spatial dimensions are moved from their last positions.

    Examples:
        >>> gt = GeoTensor(np.random.rand(3, 4, 100, 200), transform, crs)
        >>> # Shape is (3, 4, 100, 200)
        >>> gt_t = gt.transpose()
        >>> # Shape is now (4, 3, 100, 200)
        >>>
        >>> # You can also specify axes explicitly:
        >>> gt_t = gt.transpose((1, 0, 2, 3))  # Valid: spatial dims remain at end
        >>> # But this would raise an error:
        >>> # gt.transpose((0, 2, 1, 3))  # Invalid: spatial dims must stay at end
    """
    ndim = len(self.shape)

    if ndim <= 2:
        # Nothing meaningful to transpose for arrays with only spatial dimensions
        return self.copy()

    # Original spatial dimensions indices
    y_dim = ndim - 2
    x_dim = ndim - 1

    if axes is None:
        # Reverse all dimensions except the spatial ones which stay at the end
        non_spatial_axes = list(range(ndim - 2))
        non_spatial_axes.reverse()
        axes = tuple(non_spatial_axes + [y_dim, x_dim])
    else:
        # Convert to tuple if necessary
        axes = tuple(axes)

        # Check if axes has the right length
        if len(axes) != ndim:
            raise ValueError(
                f"axes should contain {ndim} dimensions, got {len(axes)}"
            )

        # Check if the last two values in axes are the spatial dimensions
        if axes[-2:] != (y_dim, x_dim):
            raise ValueError(
                "Cannot change the position of spatial dimensions. "
                f"The last two axes must be {y_dim} and {x_dim}."
            )

    # Perform the transpose
    transposed_values = np.transpose(self.values, axes)

    return GeoTensor(
        transposed_values,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=self.fill_value_default,
        attrs=self.attrs,
    )

valid_footprint(crs=None, method='all')

vectorizes the valid values of the GeoTensor and returns the footprint as a Polygon.

Parameters:

Name	Type	Description	Default
crs	Optional[str]	Coordinate reference system. Defaults to None.	None
method	str	"all" or "any" to aggregate the channels of the image. Defaults to "all".	'all'

Returns:

Type	Description
Union[MultiPolygon, Polygon]	Polygon or MultiPolygon: footprint of the GeoTensor.

Source code in georeader/geotensor.py

def valid_footprint(
    self, crs: Optional[str] = None, method: str = "all"
) -> Union[MultiPolygon, Polygon]:
    """
    vectorizes the valid values of the GeoTensor and returns the footprint as a Polygon.

    Args:
        crs (Optional[str], optional): Coordinate reference system. Defaults to None.
        method (str, optional): "all" or "any" to aggregate the channels of the image. Defaults to "all".

    Returns:
        Polygon or MultiPolygon: footprint of the GeoTensor.
    """
    valid_values = self.values != self.fill_value_default
    if len(valid_values.shape) > 2:
        if method == "all":
            valid_values = np.all(
                valid_values,
                axis=tuple(np.arange(0, len(valid_values.shape) - 2).tolist()),
            )
        elif method == "any":
            valid_values = np.any(
                valid_values,
                axis=tuple(np.arange(0, len(valid_values.shape) - 2).tolist()),
            )
        else:
            raise NotImplementedError(
                f"Method {method} to aggregate channels not implemented"
            )

    from georeader import vectorize

    polygons = vectorize.get_polygons(valid_values, transform=self.transform)
    if len(polygons) == 0:
        raise ValueError("GeoTensor has no valid values")
    elif len(polygons) == 1:
        pol = polygons[0]
    else:
        pol = MultiPolygon(polygons)
    if crs is None:
        return pol

    return window_utils.polygon_to_crs(pol, self.crs, crs)

validmask()

Returns a mask of the valid values of the GeoTensor. The mask is a boolean array with the same shape as the GeoTensor values, where True indicates valid values and False indicates invalid values. The mask is created by comparing the values of the GeoTensor with the self.fill_value_default.

Returns:

Name	Type	Description
Self	Self	GeoTensor with the valid boolean mask.

Source code in georeader/geotensor.py

def validmask(self) -> Self:
    """
    Returns a mask of the valid values of the GeoTensor. The mask is a boolean array
    with the same shape as the GeoTensor values, where True indicates valid values and
    False indicates invalid values.
    The mask is created by comparing the values of the GeoTensor with the `self.fill_value_default`.

    Returns:
        Self: GeoTensor with the valid boolean mask.
    """
    if self.fill_value_default is None:
        return GeoTensor(
            np.ones(self.shape, dtype=bool),
            transform=self.transform,
            crs=self.crs,
            fill_value_default=self.fill_value_default,
            attrs=self.attrs,
        )
    return GeoTensor(
        values=self.values != self.fill_value_default,
        transform=self.transform,
        crs=self.crs,
        fill_value_default=False,
        attrs=self.attrs,
    )

write_from_window(data, window)

Writes array to GeoTensor values object at the given window position. If window surpasses the bounds of this object it crops the data to fit the object.

Parameters:

Name	Type	Description	Default
data	ndarray	Tensor to write. Expected: spatial dimensions window.width, window.height. Rest: same as self	required
window	Window	Window object that specifies the spatial location to write the data	required

Examples:

>>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
>>> data = np.random.rand(3, 50, 50)
>>> window = rasterio.windows.Window(col_off=7, row_off=9, width=50, height=50)
>>> gt.write_from_window(data, window)

Source code in georeader/geotensor.py

def write_from_window(self, data: np.ndarray, window: rasterio.windows.Window):
    """
    Writes array to GeoTensor values object at the given window position. If window surpasses the bounds of this
    object it crops the data to fit the object.

    Args:
        data: Tensor to write. Expected: spatial dimensions `window.width`, `window.height`. Rest: same as `self`
        window: Window object that specifies the spatial location to write the data

    Examples:
        >>> gt = GeoTensor(np.random.rand(3, 100, 100), transform, crs)
        >>> data = np.random.rand(3, 50, 50)
        >>> window = rasterio.windows.Window(col_off=7, row_off=9, width=50, height=50)
        >>> gt.write_from_window(data, window)

    """
    window_data = rasterio.windows.Window(
        col_off=0, row_off=0, width=self.width, height=self.height
    )
    if not rasterio.windows.intersect(window, window_data):
        return

    assert data.shape[-2:] == (
        window.height,
        window.width,
    ), f"window {window} has different shape than data {data.shape}"
    assert (
        data.shape[:-2] == self.shape[:-2]
    ), f"Dimension of data in non-spatial channels found {data.shape} expected: {self.shape}"

    slice_dict, pad_width = window_utils.get_slice_pad(window_data, window)
    slice_list = self._slice_tuple(slice_dict)
    # need_pad = any(p != 0 for p in pad_width["x"] + pad_width["y"])

    slice_data_spatial_x = slice(
        pad_width["x"][0], None if pad_width["x"][1] == 0 else -pad_width["x"][1]
    )
    slice_data_spatial_y = slice(
        pad_width["y"][0], None if pad_width["y"][1] == 0 else -pad_width["y"][1]
    )
    slice_data = self._slice_tuple(
        {"x": slice_data_spatial_x, "y": slice_data_spatial_y}
    )
    self.values[slice_list] = data[slice_data]