Raster

Raster#

Raster is the GeoUtils object for rasters. Its API is mirrored in the rst, which exposes the same methods and attributes on xarray.DataArray objects.

For the main raster API reference page, see Raster API. Below is the full class description specific to Raster:

class geoutils.Raster(filename_or_dataset, bands=None, is_mask=False, load_data=False, parse_sensor_metadata=False, silent=True, downsample=1, force_nodata=None)[source]

Bases: RasterBase

The georeferenced raster.

Main attributes:
data: np.ndarray

Data array of the raster, with dimensions corresponding to (count, height, width).

transform: affine.Affine

Geotransform of the raster.

crs: pyproj.crs.CRS

Coordinate reference system of the raster.

nodata: int or float

Nodata value of the raster.

All other attributes are derivatives of those attributes, or read from the file on disk. See the API for more details.

property area_or_point: Literal['Area', 'Point'] | None

Pixel interpretation of the raster.

Based on the “AREA_OR_POINT” raster metadata:

  • If pixel interpretation is “Area”, the value of the pixel is associated with its upper left corner.

  • If pixel interpretation is “Point”, the value of the pixel is associated with its center.

When setting with self.area_or_point = new_area_or_point, uses the default arguments of self.set_area_or_point().

astype(dtype, convert_nodata=True, inplace=False)[source]

Convert data type of the raster.

Overloads:

By default, converts the nodata value to the default of the new data type.

Parameters:

  • dtype (DTypeLike) – Any numpy dtype or string accepted by numpy.astype.

  • convert_nodata (bool) – Whether to convert the nodata value to the default of the new dtype.

  • inplace (bool) – Whether to modify the raster in-place.

Returns:

Raster with updated dtype (or None if inplace).

property bands: tuple[int, ...]

Band indexes loaded in memory if they are, otherwise on disk.

property bbox: BoundingBox

Bounding box of the raster.

property bounds: BoundingBox

Bounding box of the raster, provided as an alias of bbox.

clip(mask, all_touched=False, mp_config=None)

Mask raster cells outside an exact clipping geometry.

The raster grid and extent stay unchanged. Cells inside the geometry keep their values, while cells outside become masked values for Raster objects and NaN values for Xarray objects.

GCP/RPC referencing requires reproject() to an affine grid first.

Parameters:

  • mask (Any) – Clipping geometry, vector, point cloud, raster or bounding box. Georeferenced masks are reprojected to this raster’s CRS.

  • all_touched (bool) – Whether to keep every cell touched by the geometry. By default, keep cells whose center is inside the geometry, following Rasterio rasterization behavior.

  • mp_config (MultiprocConfig | None) – Multiprocessing configuration.

Return type:

RasterType

Returns:

Raster with cells outside the clipping geometry as NaN or masked.

coords(grid=True, shift_area_or_point=None, force_offset=None, crs=None)

Get coordinates (x,y) of all pixels in the raster.

Parameters:

  • grid (bool) – Whether to return coordinate matrix, otherwise only X/Y axis coordinates. Required for rotated rasters where X/Y axis coordinates do not suffice to represent grid coordinates.

  • shift_area_or_point (bool | None) – Whether to shift with pixel interpretation, which shifts to center of pixel coordinates if self.area_or_point is “Point” and maintains corner pixel coordinate if it is “Area” or None. Defaults to True. Can be configured with the global setting geoutils.config[“shift_area_or_point”].

  • force_offset (str | None) – Ignore pixel interpretation and force coordinate to a certain offset: “center” of pixel, or any corner (upper-left “ul”, “ur”, “ll”, lr”). Default coordinate of a raster is upper-left.

  • crs (CRS | int | None) – Coordinate reference system in which the coordinates are provided. Defaults to CRS of the input raster.

Returns x,y:

Arrays of the (x,y) coordinates.

Return type:

tuple[NDArray[floating[Any] | integer[Any]], NDArray[floating[Any] | integer[Any]]]

copy(new_array=None, cast_nodata=True, deep=True)[source]

Copy the raster in-memory.

Parameters:

  • new_array (NDArray[floating[Any] | integer[Any]] | None) – New array to use in the copied raster.

  • cast_nodata (bool) – Automatically cast nodata value to the default nodata for the new array type if not compatible. If False, will raise an error when incompatible.

  • deep (bool) – If True, will return a deep copy of the raster.

Return type:

RasterType

Returns:

Copy of the raster.

cosample(other, *, band=1, other_band=1, auxiliary=None, auxiliary_at=None, at=None, mask=None, mask_mode='inside', subsample=1, random_state=None, strategy='topk', raster_point_mode=None, grid_method='linear', resample_method='linear', grid_kwargs=None, resample_kwargs=None, align='raise', mp_config=None)

Sample this raster and another dataset at common finite locations.

By default, raster values are resampled at point locations when either primary input is a point cloud. Use raster_point_mode="grid_points" to grid point values onto the raster instead. An explicit at chooses the exact output locations and must agree with any explicit mode. Raw auxiliary arrays must identify the primary input whose grid or point ordering they follow.

Spatial inputs, explicit output support and raster or point masks must use one family: Raster/PointCloud objects, or DataArray/GeoDataFrame objects. The latter may mix eager and Dask storage. Plain arrays and vector outlines are accepted with either family.

Parameters:

  • other (RasterBase | DataArray | PointCloudBase | GeoDataFrame | ArrayLike) – Dataset to sample alongside this raster. A plain array follows this raster’s grid.

  • band (int) – Band selected from this raster, counting from one.

  • other_band (int) – Band selected from other if it is a raster, counting from one.

  • auxiliary (Mapping[str, Any] | None) – Additional values by output name (e.g. {“slope”: slope_raster}). Select a raster band with {“slope”: (slope_raster, 2)} or a point column with {“intensity”: (points, “intensity”)}. Spatial inputs without a selector use the first raster band or active point values.

  • auxiliary_at (Literal['self', 'other'] | Mapping[str, Literal['self', 'other']] | None) – Input locations followed by plain auxiliary arrays: “self”, “other”, or a choice per name (e.g. {“slope”: “other”}). Spatial auxiliaries use their own coordinates.

  • at (Literal['self', 'other'] | RasterBase | DataArray | PointCloudBase | GeoDataFrame | None) – Output locations: “self”, “other”, or a reference raster/point cloud. Defaults to the other point cloud’s locations when present, or this raster’s grid otherwise.

  • mask (RasterBase | DataArray | geoutils.vector.base.VectorBase | GeoDataFrame | ArrayLike | None) – Locations eligible for sampling, defined by a boolean array, spatial mask, or vector outlines.

  • mask_mode (Literal['inside', 'outside']) – Whether a vector mask keeps locations “inside” or “outside” its geometries.

  • subsample (int | float) – Fraction of common finite locations (e.g. 0.1), or maximum count (e.g. 1000); 1 keeps all.

  • random_state (int | Generator | None) – Seed or random generator for reproducible sampling (e.g. 42).

  • strategy (Literal['sequential', 'topk']) – Sampling with “topk” or “sequential”; “topk” keeps the same seeded sample across chunk sizes.

  • raster_point_mode (Literal['grid_points', 'resample_raster'] | None) – Conversion direction: “grid_points” places points on a raster, “resample_raster” reads rasters at points. Defaults to at’s locations, or point locations when available. Must agree with at.

  • grid_method (Literal['nearest', 'linear', 'cubic', 'idw'] | Literal['mean', 'average', 'minimum', 'min', 'maximum', 'max', 'range', 'count', 'stdev', 'average_distance', 'average_distance_pts'] | Interpolator | Reducer) – Point gridding by SciPy interpolation (“nearest”, “linear”, “cubic”), or circular “idw”, “mean”, “minimum”, “maximum”, “range”, “count”, “stdev”, “average_distance”, “average_distance_pts”. The aliases “average”, “min” and “max” select “mean”, “minimum” and “maximum”.

  • resample_method (Literal['nearest', 'linear', 'cubic', 'quintic', 'slinear', 'pchip', 'splinef2d'] | Literal['reduce']) – Raster interpolation using the SciPy methods “nearest”, “linear”, “cubic”, “quintic”, “slinear”, “pchip” or “splinef2d”. Window reduction (“reduce”) is not implemented.

  • grid_kwargs (Mapping[str, Any] | None) – Options for PointCloud.grid(), e.g. {“radius”: 2, “min_points”: 3} sets a two-pixel radius and minimum of three finite points for circular methods. Other options include “distance_power” for IDW and “engine” (“scipy” or “numba”). Set locations and method with at and grid_method.

  • resample_kwargs (Mapping[str, Any] | None) – Options for Raster.interp_points(), e.g. {“nodata_handling”: “ignore”}. The choices are “nearest”, “ignore”, “propagate”, a half-order rule, or a non-negative distance in pixels. Set locations, band and method with the corresponding cosample() arguments.

  • align (Literal['raise', 'reproject']) – Handling of mismatched grids or coordinate systems: “raise” an error, or “reproject” to match at. Point inputs must still share the same ordered coordinates when sampled at points.

  • mp_config (MultiprocConfig | None) – Worker and tile settings for multiprocessing. Raster output uses its outfile; cannot be combined with Dask inputs.

Return type:

RasterBase | DataArray | PointCloudBase | GeoDataFrame

Returns:

Raster or point cloud on the selected support; Xarray DataArray or eager/lazy GeoDataFrame for accessor calls. Bands or columns contain “self”, “other”, then auxiliaries in mapping order. Raster band names are stored in tags["long_name"] (Xarray attrs["long_name"]). Point outputs use “self” as their active data attribute and preserve the support index and order.

property count: int

Count of bands loaded in memory if they are, otherwise the one on disk.

crop(bbox, inplace=False)

Crop the raster to a given extent.

Match-reference: a reference raster or vector can be passed to match bounds during cropping.

Cropping preserves the original pixel resolution, cropping to the extent that most closely aligns with the current coordinates. To match the extent of another dataset exactly, use reproject().

GCP/RPC referencing requires reproject() to an affine grid first. To select image pixels directly, use icrop().

Parameters:

  • bbox (RasterType | VectorType | list[float] | tuple[float, ...]) – Geometry to crop raster to. Can use either a raster or vector as match-reference, or a list of coordinates. If bbox is a raster or vector, will crop to the bounds. If bbox is a list of coordinates, the order is assumed to be [xmin, ymin, xmax, ymax].

  • inplace (bool) – (DEPRECATED. Use rast = rast.crop() instead) Whether to crop in-place or not.

Return type:

RasterType | None

Returns:

A new cropped raster.

property crs: CRS

Coordinate reference system of the raster.

Returns:

Coordinate reference system.

property data: MaskedArray[Any, dtype[floating[Any] | integer[Any]]]

Array of the raster.

Returns:

Raster array.

Setter accepts:

NDArray[floating[Any] | integer[Any]] | MaskedArray[Any, dtype[floating[Any] | integer[Any]]]

property driver: str | None

Driver used to read a file on disk.

property dtype: DTypeLike

Data type of the raster.

edit(*, crs=<object object>, transform=<object object>, nodata=<object object>, tags=<object object>, area_or_point=<object object>)

Return a copy with selected raster metadata edited together.

This groups the core metadata changes provided by gdal_edit without changing pixel values or the source raster. Omitting an argument retains its existing value, while explicitly passing None clears CRS, nodata, tags or pixel interpretation.

Parameters:

  • crs (CRS | int | str | None | object) – New coordinate reference system.

  • transform (Affine | tuple[float, ...] | object) – New affine geotransform.

  • nodata (int | float | None | object) – New nodata metadata value without changing pixel values or their mask.

  • tags (dict[str, Any] | None | object) – Metadata tags to update, or None to clear all tags.

  • area_or_point (Literal['Area', 'Point'] | None | object) – New pixel interpretation without shifting the geotransform.

Return type:

RasterType

Returns:

Raster copy with updated metadata.

estimate_error_structure(other, *, other_precision='same', predictors=None, components=None, mask=None, mp_config=None, **estimate_kwargs)

Estimate this raster’s error structure from its difference with another dataset.

The two inputs are aligned by cosample(). With other_precision="same", their errors must be independent and have the same magnitude and correlation; the difference is divided by the square root of two. Use "negligible" when the other dataset’s error can be ignored.

Parameters:

  • other (RasterBase | DataArray | PointCloudBase | GeoDataFrame | ArrayLike) – Raster, point cloud, or array of comparable measurements.

  • other_precision (Literal['same', 'negligible']) – Whether the other input has the same precision or negligible error.

  • predictors (Mapping[str, Any] | None) – Named variables controlling a magnitude that varies with location.

  • components (Mapping[str, Mapping[str, Any]] | None) – Named magnitude and correlation models to fit.

  • mask (RasterBase | DataArray | geoutils.vector.base.VectorBase | GeoDataFrame | ArrayLike | None) – Locations used to estimate the error structure.

  • mp_config (MultiprocConfig | None) – Worker and tile settings for multiprocessing.

  • estimate_kwargs (Any) – Additional options accepted by ErrorStructure.estimate().

Return type:

ErrorStructure

Returns:

Fitted error structure for this raster.

fill_nodata(max_search_distance=100.0, smoothing_iterations=0, interpolation=None, mask=None, *, resampling=None, backend='auto')

Fill nodata cells from nearby finite values.

This mirrors gdal_fillnodata for each band. Cells beyond the search distance remain nodata, and the complete array is processed eagerly so interpolation can reach across chunk boundaries.

Parameters:

  • max_search_distance (float) – Maximum search distance in pixels (defaults to 100).

  • smoothing_iterations (int) – Number of 3 by 3 smoothing passes after filling (defaults to 0).

  • interpolation (Literal['inv_dist', 'nearest'] | None) – Legacy alias using “inv_dist” or “nearest”.

  • resampling (Literal['idw', 'inv_dist', 'nearest'] | None) – “idw” (default) for inverse distance weighting, or “nearest”.

  • backend (Literal['auto', 'gdal', 'scipy', 'numba', 'gstools', 'gpytorch']) – “auto” or “gdal” for GDAL gap filling.

  • mask (RasterBase | DataArray | NDArray[bool] | None) – Optional Boolean raster or array where True values can provide interpolation values.

Return type:

RasterType

Returns:

Raster copy with reachable nodata cells filled.

filter(method, window=None, sigma=1, engine=None, outlier_threshold=2.0, mp_config=None, coverage=None, window_shape=None, backend='auto', method_kwargs=None, **kwargs)

Apply a filter to the raster array.

Choose a built-in filter with method and set its neighborhood with window, window_shape and coverage. Advanced callers can supply a Reducer or custom callable. A Reducer can select cells with GridNeighbours, for example filter(Mean(neighborhood=GridNeighbours(size=5, shape="circular", coverage="fractional"))). Its configured window is used unless window, window_shape or coverage is supplied. size and kernel_shape remain deprecated aliases for window and window_shape.

Parameters:

  • method (Literal['nearest', 'linear', 'bilinear', 'cubic', 'cubic_convolution', 'cubic_spline', 'lanczos', 'quintic', 'slinear', 'pchip', 'splinef2d', 'idw', 'inv_dist', 'kriging', 'mean', 'average', 'sum', 'median', 'med', 'quantile', 'q1', 'q3', 'mode', 'minimum', 'min', 'maximum', 'max', 'range', 'count', 'stdev', 'std', 'standard_deviation', 'rms', 'rmse', 'average_distance', 'average_distance_pts'] | Literal['gaussian', 'distance'] | Reducer | Callable[..., NDArray[floating[Any] | integer[Any]]]) – Built-in filter or statistic name (e.g. “gaussian”, “median”, “mean”, “max”, “min”, “distance”), a Reducer, or a custom callable that takes a 2D ndarray and returns one.

  • backend (Literal['auto', 'gdal', 'scipy', 'numba', 'gstools', 'gpytorch']) – “auto” uses SciPy; “scipy” and “numba” select the calculation explicitly.

  • method_kwargs (dict[str, Any] | None) – Settings for a named statistic, such as a quantile q.

  • window (int | None) – Window size for filter, in source pixels.

  • mp_config (MultiprocConfig | None) – Multiprocessing configuration.

  • coverage (Literal['center', 'all_touched', 'fractional'] | None) – Select cells by center, include all touched cells, or weight by covered area.

  • window_shape (Literal['square', 'circular'] | None) – Use a square or circular window for a Reducer.

  • sigma (int) – Optional standard deviation for Gaussian filtering. Only used when method=”gaussian”.

  • engine (Literal['scipy', 'numba'] | None) – Legacy alias for backend (“scipy” or “numba”).

  • outlier_threshold (float) – The minimum difference abs(array - mean) for a pixel to be considered an outlier. Only used when method=”distance”.

  • kwargs (Any) – Additional keyword arguments passed to the underlying filter implementation.

Return type:

RasterType

Returns:

A new Raster instance with the filtered data (or None if inplace).

Raises:

  • ValueError – If the filter name is not one of the predefined options.

  • TypeError – If method is neither a string nor a callable.

property footprint: Vector

Footprint of the raster.

classmethod from_array(data, transform, crs, nodata=None, area_or_point=None, tags=None, cast_nodata=True, *, gcps=None, rpcs=None)[source]

Create a raster from a numpy array and the georeferencing information.

Expects a 2D (single band) or 3D (multi-band) array of the raster. The first axis corresponds to bands.

Parameters:

  • data (NDArray[floating[Any] | integer[Any]] | MaskedArray[Any, dtype[floating[Any] | integer[Any]]] | NDArray[bool]) – Input array, 2D for single band or 3D for multi-band (bands should be first axis).

  • transform (tuple[float, ...] | Affine) – Affine 2D transform. Either a tuple(x_res, 0.0, top_left_x, 0.0, y_res, top_left_y) or an affine.Affine object.

  • crs (CRS | int | None) – Coordinate reference system. Any CRS supported by Pyproj (e.g., CRS object, EPSG integer).

  • nodata (int | float | None) – Nodata value.

  • area_or_point (Literal['Area', 'Point'] | None) – Pixel interpretation of the raster, will be stored in AREA_OR_POINT metadata.

  • tags (dict[str, Any]) – Metadata stored in a dictionary.

  • cast_nodata (bool) – Automatically cast nodata value to the default nodata for the new array type if not compatible. If False, will raise an error when incompatible.

  • gcps (tuple[list[GroundControlPoint], CRS | None] | None) – Ground control points and their CRS as a (points, CRS) tuple, following Rasterio.

  • rpcs (RPC | None) – Rational polynomial coefficients referenced to WGS84.

Return type:

RasterType

Returns:

Raster created from the provided array and georeferencing.

Example:

You have a data array in EPSG:32645. It has a spatial resolution of 30 m in x and y, and its top left corner is X=478000, Y=3108140.

>>> data = np.ones((500, 500), dtype="uint8")
>>> transform = (30.0, 0.0, 478000.0, 0.0, -30.0, 3108140.0)
>>> myim = Raster.from_array(data, transform, 32645)
classmethod from_pointcloud_regular(pointcloud, grid_coords=None, transform=None, shape=None, nodata=None, data_name='b1', area_or_point='Point')

Create a raster from a point cloud with coordinates on a regular grid.

To inform on what grid to create the raster, either pass a tuple of X/Y grid coordinates, or the expected transform and shape. All point cloud coordinates must fall exactly at one of the coordinates of this grid.

Parameters:

  • pointcloud (PointCloudBase | GeoDataFrame | DataArray) – Point cloud.

  • grid_coords (tuple[NDArray[floating[Any] | integer[Any]], NDArray[floating[Any] | integer[Any]]]) – Regular coordinate vectors for the raster, from which the geotransform and shape are deduced.

  • transform (Affine) – Geotransform of the raster.

  • shape (tuple[int, int]) – Shape of the raster.

  • nodata (int | float | None) – Nodata value of the raster.

  • data_name (str) – Name to use for point cloud data attribute, defaults to “bX” where X is the data band number.

  • area_or_point (Literal['Area', 'Point']) – Whether to set the pixel interpretation of the raster to “Area” or “Point”.

Return type:

RasterType

classmethod from_xarray(ds, dtype=None)[source]

Create raster from a xarray.DataArray.

This conversion loads the xarray.DataArray in memory. Use functions of the Xarray accessor directly to avoid this behaviour.

Parameters:
Return type:

RasterType

Returns:

Raster.

property gcps: tuple[list[GroundControlPoint], CRS | None]

Ground control points and their CRS, following Rasterio’s (points, CRS) convention.

georeferenced_grid_equal(other, warn_3d_crs=True)

Check that raster shape, geotransform and CRS are equal.

Parameters:

  • other (RasterType) – Another raster.

  • warn_3d_crs (bool) – Whether to warn if 3D CRS differs.

Return type:

bool

Returns:

Whether the two objects have the same georeferenced grid.

get_bbox_projected(out_crs, densify_points=5000)

Get the raster bounding box projected in a specified CRS.

Parameters:

  • out_crs (CRS) – Output CRS.

  • densify_points (int) – Maximum points to be added between image corners to account for non linear edges. Reduce if time computation is really critical (ms) or increase if extent is not accurate enough.

Return type:

BoundingBox

get_footprint_projected(out_crs, densify_points=5000)

Get raster footprint projected in a specified CRS.

The polygon points of the vector are densified during reprojection to warp the rectangular square footprint of the original projection into the new one.

Parameters:

  • out_crs (CRS) – Output CRS.

  • densify_points (int) – Maximum points to be added between image corners to account for non linear edges. Reduce if time computation is really critical (ms) or increase if extent is not accurate enough.

Return type:

Vector

get_mask()[source]

Get mask of invalid values from the raster.

If the raster is not loaded, reads only the mask from disk to optimize memory usage.

The mask is always returned as a boolean array, even if there is no mask associated to .data (nomask property of masked arrays).

Return type:

NDArray[bool]

Returns:

The mask of invalid values in the raster.

get_metric_crs(local_crs_type='universal', method='centroid', *, center=None)

Get a universal or custom local metric coordinate reference system for the raster.

Parameters:

  • local_crs_type (Literal['universal'] | Literal['custom']) – Whether to get a “universal” local CRS (UTM or UPS) or a “custom” local CRS (stereographic centered on the footprint centroid).

  • method (Literal['centroid'] | Literal['geopandas']) – Method to choose the zone of the CRS, either based on the centroid of the footprint or the extent as implemented in geopandas.GeoDataFrame.estimate_utm_crs(). Ignored for a custom CRS.

  • center (tuple[float, float] | None) – Optional (longitude, latitude) in WGS84 degrees for a custom CRS. By default, use the planar footprint centroid after unwrapping longitudes across the antimeridian. Supply a center for footprints that span a hemisphere or surround a pole.

Return type:

CRS

Returns:

A projected CRS with metre units. A custom CRS has true scale at its center and suits compact regions.

get_nanarray(floating_dtype='float32', *, return_mask=False)

Call to_nanarray() through its deprecated name.

Overloads:
get_stats(stats_name=None, inlier_mask=None, band=None)

Call stats() with the legacy argument names; deprecated in favor of stats().

Return type:

Any

property height: int

Height of the raster in pixels.

icrop(bbox, inplace=False)

Crop raster based on pixel indices (bbox), converting them into georeferenced coordinates.

Parameters:

  • bbox (list[int] | tuple[int, ...]) – Bounding box based on indices of the raster array (colmin, rowmin, colmax, rowmax).

  • inplace (bool) – (DEPRECATED. Use rast = rast.icrop() instead) Whether to crop in-place or not.

Return type:

RasterType | None

Returns:

Cropped raster.

ij2xy(i, j, shift_area_or_point=None, force_offset=None)

Get coordinates (x,y) of indexes (row,column).

By default, the indexes are shifted with the interpretation of pixel coordinates “AREA_OR_POINT” of the raster, to ensure that the indexes of points represent the right location. See parameter description of shift_area_or_point for more details.

This function is reversible with xy2ij for any pixel interpretation.

Parameters:

  • i (ArrayLike) – Row (i) index of pixel.

  • j (ArrayLike) – Column (j) index of pixel.

  • shift_area_or_point (bool | None) – Whether to shift with pixel interpretation, which shifts to center of pixel coordinates if self.area_or_point is “Point” and maintains corner pixel coordinate if it is “Area” or None. Defaults to True. Can be configured with the global setting geoutils.config[“shift_area_or_point”].

  • force_offset (str | None) – Ignore pixel interpretation and force coordinate to a certain offset: “center” of pixel, or any corner (upper-left “ul”, “ur”, “ll”, lr”). Default coordinate of a raster is upper-left.

Returns x, y:

x,y coordinates of i,j in reference system.

Return type:

tuple[NDArray[floating[Any] | integer[Any]], NDArray[floating[Any] | integer[Any]]]

property indexes: tuple[int, ...]

Band indexes (duplicate of .bands attribute, mirroring Rasterio naming “indexes”). Loaded in memory if they are, otherwise on disk.

info(stats=False, verbose=True)

Print summary information about the raster.

Overloads:
Parameters:

  • stats (bool) – Add statistics for each band of the dataset (max, min, median, mean, std. dev.). Default is to not calculate statistics.

  • verbose (bool) – If set to True (default) will directly print to screen and return None

Returns:

Summary string or None.

interp_at_points(points, method=None, nodata_handling=None, band=1, input_latlon=False, as_array=False, shift_area_or_point=None, mp_config=None, error_structure=None, *, resampling=None, backend='auto', resampling_kwargs=None, **kwargs)

Interpolate raster values at a set of points.

Overloads:

Returns a point cloud containing the interpolated values at the point coordinates, or optionally an array.

Uses scipy.ndimage.map_coordinates if the Raster is on an equal grid using “nearest” or “linear” (for speed), otherwise uses scipy.interpn on a regular grid.

Optionally, user can enforce the interpretation of pixel coordinates in self.tags[‘AREA_OR_POINT’] to ensure that the interpolation of points is done at the right location. See parameter description of shift_area_or_point for more details.

This method forwards its options to resample_at_points().

Parameters:

  • points (tuple[NDArray[floating[Any] | integer[Any]], NDArray[floating[Any] | integer[Any]]] | tuple[int | float | integer[Any] | floating[Any], int | float | integer[Any] | floating[Any]] | PointCloudBase | GeoDataFrame | DataArray) – Point(s) at which to interpolate raster value. Can be either a tuple of array-like of X/Y coordinates (same CRS as raster or latitude/longitude, see “input_latlon”) or a pointcloud in any CRS. If points fall outside of image, value returned is nan.

  • method (Literal['nearest', 'linear', 'cubic', 'quintic', 'slinear', 'pchip', 'splinef2d'] | Interpolator) – Interpolator or one of ‘nearest’, ‘linear’, ‘cubic’, ‘quintic’, ‘slinear’, ‘pchip’ or ‘splinef2d’. String methods use the existing SciPy array functions. A custom Interpolator defines the calculation for source cells at fixed GridNeighbours offsets. Defaults to the configured interpolation method.

  • resampling (Literal['nearest', 'linear', 'bilinear', 'cubic', 'cubic_convolution', 'cubic_spline', 'lanczos', 'quintic', 'slinear', 'pchip', 'splinef2d', 'idw', 'inv_dist', 'kriging', 'mean', 'average', 'sum', 'median', 'med', 'quantile', 'q1', 'q3', 'mode', 'minimum', 'min', 'maximum', 'max', 'range', 'count', 'stdev', 'std', 'standard_deviation', 'rms', 'rmse', 'average_distance', 'average_distance_pts'] | Interpolator | None) – Built-in interpolation name, such as “nearest” or “linear”, or an advanced Interpolator. method remains an alias; defaults to the configured interpolation method.

  • backend (Literal['auto', 'gdal', 'scipy', 'numba', 'gstools', 'gpytorch']) – “auto” uses SciPy; “gstools” and “gpytorch” select a kriging solver.

  • resampling_kwargs (dict[str, Any] | None) – Settings for the named interpolation, including a fitted model for kriging.

  • nodata_handling (Literal['nearest', 'ignore', 'propagate'] | Literal['half_order_up', 'half_order_down'] | int | None) – "nearest" masks values whose nearest source cell is missing; "ignore" uses available finite cells; "propagate" masks values using a missing source cell. A non-negative integer or half-order choice instead masks cells within that distance of missing source values. Defaults to geoutils.config[“interpolation_nodata_handling”], which is "nearest" unless changed.

  • band (int) – Band to use (from 1 to self.count).

  • input_latlon (bool) – (Only for tuple point input) Whether to convert input coordinates from latlon to raster CRS.

  • as_array (bool) – Whether to return an array of interpolated values instead of a point cloud.

  • shift_area_or_point (bool | None) – Whether to shift with pixel interpretation, which shifts to center of pixel coordinates if self.area_or_point is “Point” and maintains corner pixel coordinate if it is “Area” or None. Defaults to True. Can be configured with the global setting geoutils.config[“shift_area_or_point”].

  • error_structure (ErrorStructure | None) – Optional observation error model used for weighting and fitting. Use geoutils.uncertainty.propagate() to calculate output uncertainty.

:returns Point cloud of interpolated points, or 1D array of interpolated values.

interp_points(points, method=None, nodata_handling=None, band=1, input_latlon=False, as_array=False, shift_area_or_point=None, mp_config=None, error_structure=None, *, resampling=None, backend='auto', resampling_kwargs=None, **kwargs)

Interpolate raster values at points using the deprecated name for interp_at_points().

Overloads:

Deprecated since version 1.1: Use interp_at_points() instead.

intersection(other, match_ref=True)

Returns the bounding box of intersection between this raster and another.

If the rasters have different projections, the intersection extent is given in self’s projection system.

Parameters:

  • other (RasterType) – Path to second raster.

  • match_ref (bool) – If set to True, returns the smallest intersection that aligns with that of self, i.e. same resolution and offset with self’s origin is a multiple of the resolution

Return type:

tuple[float, float, float, float]

Returns:

Extent of the intersection between the 2 raster (xmin, ymin, xmax, ymax) in self’s coordinate system.

property is_loaded: bool

Whether the raster array is loaded.

property is_mask: bool

Whether the raster array is a boolean data type (a mask).

krige(variogram, ref=None, *, backend='gstools', max_overlap=None, exact=True, pseudo_inverse=True, mp_config=None, error_structure=None)

Interpolate raster values by ordinary kriging on this grid or a reference grid.

The longest fitted variogram range defines how far GeoUtils searches for source cells. max_overlap can limit that distance and the extra source data read around chunk edges. Repeated raster neighborhoods reuse previously solved weights, while edges and missing-data patterns get weights for their exact source cells.

Parameters:

  • variogram (Variogram | VariogramModel) – Fitted GeoUtils variogram used to calculate covariance.

  • ref (RasterBase | DataArray | None) – Optional destination raster grid. By default, interpolate on this raster’s own grid.

  • backend (Literal['gstools', 'gpytorch']) – Library used to solve the kriging system, "gstools" or "gpytorch".

  • max_overlap (float | None) – Optional physical-distance cap on the modeled decorrelation range.

  • exact (bool) – Return a source value exactly when a target has the same coordinates.

  • pseudo_inverse (bool) – Let GSTools use a pseudo inverse when source coordinates repeat or are redundant.

  • mp_config (MultiprocConfig | None) – Optional multiprocessing output configuration.

  • error_structure (ErrorStructure | None) – Optional observation errors added to the covariance used for the kriging fit.

Return type:

RasterType

Returns:

Kriged raster. Use geoutils.uncertainty.propagate() to calculate output uncertainty.

load(bands=None, **kwargs)[source]

Load the raster array from disk.

Parameters:

  • bands (int | list[int] | None) – Band(s) to load. Note that rasterio begins counting at 1, not 0.

  • kwargs (Any) – Optional keyword arguments sent to ‘_load_rio()’.

Raises:
Return type:

None

property name: str | None

Name of the file on disk, if it exists.

property nodata: int | float | None

Nodata value of the raster.

When setting with self.nodata = new_nodata, uses the default arguments of set_nodata().

Returns:

Nodata value.

outside_image(xi, yj, index=True)

Check whether a given point falls outside the bounds of the raster.

Parameters:

  • xi (ArrayLike) – Indices (or coordinates) of x direction to check.

  • yj (ArrayLike) – Indices (or coordinates) of y direction to check.

  • index (bool) – Interpret ij as raster indices (default is True). If False, assumes ij is coordinates.

Returns is_outside:

True if ij is outside the bounds.

Return type:

bool

pairsample(*, band=1, n_pairs=1000000, sampling='loglag', min_distance=None, max_distance=None, random_state=None, mask=None, strategy='chunk_anchors', deduplicate='per_anchor', batch_pairs=2000000, max_rounds=50, max_oversample=8.0, chunks_per_round=8, anchors_per_round=20000, distances_per_anchor=8, angles_per_distance=8, hybrid_local_fraction=0.0, max_local_distance=None, index_dtype=<class 'numpy.int32'>, distance_dtype=<class 'numpy.float64'>, mp_config=None)

Sample cell pairs in the raster.

This function provides different strategies for sampling short and long pairwise distances in large rasters. It supports chunked Dask and Multiprocessing out-of-memory reads, with explicit pair subsampling to also limit the returned data held in memory.

Sampling methods#

With the default sampling="loglag", distances are drawn across a logarithmic scale and directions are drawn uniformly around the first cell, so that short and long distances are both represented. The resulting offsets are rounded to the raster grid. The "independent" strategy draws a new first cell for every pair, while "anchors" reuses a limited set of first cells. The default "chunk_anchors" also draws those first cells from a limited number of source chunks to reduce reads. The "anchor_batched" strategy draws several distances and directions together for each first cell.

With sampling="random_xy", both cells are drawn independently and pairs outside the requested distance range are discarded. Pairs near min_distance or max_distance are usually rare, even though these distance extremes are often important for spatial analysis. Log-lag strategy options do not apply; batch_pairs and max_rounds control how candidates are collected.

Memory and chunked inputs#

n_pairs limits the returned data held in memory, while batch_pairs limits temporary candidate arrays. For Dask and Multiprocessing inputs, finite cells are counted separately by chunk or tile. Pair geometry is generated from the raster grid, then only chunks or tiles containing candidate endpoints are read. The complete raster band is not collected in memory and the Raster stays unloaded. Pair selection and the returned Xarray Dataset are eager for both backends.

type band:

int

param band:

Band to sample, counting from one.

type n_pairs:

int

param n_pairs:

Requested number of pairs with two finite values; fewer may be returned if sampling stops early.

type sampling:

Literal['loglag', 'random_xy']

param sampling:

"loglag" balances short and long distances on a log scale; "random_xy" draws endpoints uniformly.

type min_distance:

float | None

param min_distance:

Smallest distance in CRS units (e.g. meters). Defaults to the smaller pixel spacing.

type max_distance:

float | None

param max_distance:

Largest distance in CRS units. Defaults to the diagonal between outermost cell centers.

type random_state:

int | Generator | None

param random_state:

Seed for reproducible sampling (e.g. 42).

type mask:

RasterBase | DataArray | geoutils.vector.base.VectorBase | GeoDataFrame | ArrayLike | None

param mask:

Eligible cells: True in a mask array or aligned mask raster, or inside vector geometries.

type strategy:

Literal['independent', 'anchors', 'chunk_anchors', 'anchor_batched']

param strategy:

GeoUtils log-lag strategy: "independent" draws each pair separately, "anchors" reuses first endpoints, "chunk_anchors" also limits source chunks, and "anchor_batched" draws several distances and directions from each first endpoint.

type deduplicate:

Literal['none', 'per_anchor', 'global']

param deduplicate:

"none" keeps repeats, "per_anchor" removes repeated targets within each anchor batch, and "global" removes repeated pairs across all batches. "random_xy" always removes repeats.

type batch_pairs:

int

param batch_pairs:

Maximum candidate pairs per batch; smaller batches use less temporary memory.

type max_rounds:

int

param max_rounds:

Maximum attempts to fill the sample after rejecting missing values or out-of-range pairs.

type max_oversample:

float

param max_oversample:

Maximum candidate count as a multiple of the target pair count (e.g. 8).

type chunks_per_round:

int

param chunks_per_round:

Maximum source chunks used when drawing anchors from selected chunks.

type anchors_per_round:

int

param anchors_per_round:

Maximum first endpoints reused per round by "anchors", "chunk_anchors", or local "independent" sampling.

type distances_per_anchor:

int

param distances_per_anchor:

Distances drawn per first endpoint with "anchor_batched".

type angles_per_distance:

int

param angles_per_distance:

Directions drawn per distance with "anchor_batched".

type hybrid_local_fraction:

float

param hybrid_local_fraction:

Fraction of candidate pairs kept within the first endpoint’s chunk (e.g. 0.5). Zero samples across the full raster; one keeps all pairs local.

type max_local_distance:

float | None

param max_local_distance:

Largest proposed local distance in CRS units. Defaults to the largest chunk diagonal.

type index_dtype:

DTypeLike

param index_dtype:

Integer NumPy dtype for returned cell indexes (e.g. "int64" for very large rasters).

type distance_dtype:

DTypeLike

param distance_dtype:

Floating NumPy dtype for returned distances (e.g. "float32" to reduce memory).

type mp_config:

MultiprocConfig | None

param mp_config:

Worker and tile settings for multiprocessing reads from an unloaded raster. Cannot be combined with Dask inputs. The returned Xarray Dataset is eager.

rtype:

Dataset

returns:

Xarray Dataset with pair and endpoint dimensions, containing cell indexes, values, coordinates, and distances.

plot(bands=None, ref=None, cmap=None, max_pixels='auto', vmin=None, vmax=None, alpha=None, title=None, cbar_title=None, add_cbar=True, ax=None, *, resampling=None, return_axes=False, savefig_fname=None, **kwargs)

Plot the raster.

Overloads:

This method performs automatic subsampling to facilitate the plotting of large datasets out-of-memory, then wraps Matplotlib imshow to which keyword arguments are passed.

Parameters:

  • bands (int | tuple[int, ...] | None) – Bands to plot, counting from 1 to self.count. Defaults to all bands.

  • ref (RasterBase | DataArray | geoutils.vector.base.VectorBase | GeoDataFrame | CRS | str | int | None) – Reference geospatial object or CRS to match. A reference object also sets the plotted axis limits to its bounds.

  • cmap (Colormap | str | None) – Colormap to use. Defaults to matplotlib.rcParams['image.cmap'].

  • max_pixels (Literal['auto'] | int | None) – The default "auto" limits output to the Matplotlib axes width and height in display pixels, as set by the figure size and DPI. An integer limits the total number of plotted pixels, and None keeps the native grid size. This creates a temporary grid through reprojection and does not select a specific overview stored in the source file. See Rasterio’s overview documentation.

  • vmin (float | int | None) – Minimum value for the colorbar. Defaults to the plotted data minimum.

  • vmax (float | int | None) – Maximum value for the colorbar. Defaults to the plotted data maximum.

  • alpha (float | int | None) – Raster and colorbar transparency.

  • title (str | None) – Plot title.

  • cbar_title (str | None) – Colorbar label.

  • add_cbar (bool) – Whether to display a colorbar. Multi-band RGB(A) plots never add one.

  • ax (Axes | Literal['new'] | None) – Matplotlib axes, "new" to create axes, or None to use the current axes.

  • resampling (Resampling | str | None) – Rasterio resampling method used when the display grid changes. Defaults to the configured reprojection method, except for boolean rasters whose reprojection uses nearest-neighbor resampling.

  • return_axes (bool) – Whether to return the image and colorbar axes.

  • savefig_fname (str | None) – Optional path at which to save the current figure.

Returns:

None, or the image axes and optional colorbar axes when return_axes=True.

polygonize(target_values='all', connectivity=4, band=1, data_name='id', strategy='label_stitch', mp_config=None)

Polygonize the raster into a vector of polygon geometries delineating target values.

Parameters:
Return type:

Vector

Returns:

Vector containing the polygonized geometries associated to target values.

property profile: dict[str, Any] | None

Basic metadata and creation options of this dataset. May be passed as keyword arguments to rasterio.open() to create a clone of this dataset.

proximity(vector=None, target_values=None, distance_unit='georeferenced', max_distance=None, mp_config=None)

Compute proximity distances to the raster target pixels, or to a vector geometry on the raster grid.

Match-reference: a raster can be passed to match its resolution, bounds and CRS for computing proximity distances.

When passing a vector, its current geometry is used. Apply a vector geometry operation first to calculate proximity to a derived geometry, for example vector.boundary.proximity(raster).

Parameters:

  • vector (VectorType | None) – Vector for which to compute the proximity to geometry, if not provided computed on this raster target pixels.

  • target_values (list[float] | None) – (Only with raster) List of target values to use for the proximity, defaults to all non-zero values.

  • distance_unit (Literal['pixel'] | Literal['georeferenced']) – Distance unit, either ‘georeferenced’ or ‘pixel’.

  • max_distance (float | None) – Largest distance to return, with farther cells set to nodata. This value is required for Dask and multiprocessing execution because it defines the overlap between chunks.

  • mp_config (MultiprocConfig | None) – Multiprocessing parameters. Cannot be combined with Dask input.

Return type:

RasterBase

Returns:

Proximity distances raster.

random_field(error_structure, *, predictors=None, n_fields=1, random_state=None, chunksizes=None, mp_config=None, backend='gpytorch', gpytorch_inducing_points='auto')

Generate one or more error fields across this raster.

Parameters:

  • error_structure (ErrorStructure) – Error model defined by magnitude and correlation components.

  • predictors (Mapping[str, Any] | None) – Named magnitude predictors aligned with this raster.

  • n_fields (int) – Number of independent fields.

  • random_state (int | Generator | None) – Seed or generator used for reproducible fields.

  • chunksizes (tuple[int, int] | None) – Lazy output chunk size for an Xarray raster. Every chunk uses the same seed for a component.

  • mp_config (MultiprocConfig | None) – Worker, tile and output file settings for a Raster.

  • backend (Literal['gstools', 'gpytorch']) – Library used to draw correlated errors. Chunked GPyTorch fields use an inducing grid.

  • gpytorch_inducing_points (int | Literal['auto'] | None) – Target grid size for approximate GPyTorch fields. “auto” uses 256 points for eager and chunked rasters; None requests an exact eager draw.

Return type:

RasterType | list[RasterType]

Returns:

One result of the same raster type, or a list when n_fields is greater than one.

raster_allclose(other, strict_masked=True, warn_failure_reason=False, rtol=1e-05, atol=1e-08)

Check if two rasters are all close.

This means that are all close: - The raster data.

And equal: - The raster’s mask, fill_value and dtype, - The raster’s transform, crs and nodata values.

The relative and absolute tolerance arguments are passed to np.allclose.

Parameters:

  • other (RasterType) – Other raster.

  • strict_masked (bool) – Whether to check if masked cells (in .data.mask) have the same value (in .data.data). Only when comparing two Rasters objects with masked arrays.

  • warn_failure_reason (bool) – Whether to warn for the reason of failure if the check does not pass.

  • rtol (float) – Relative tolerance.

  • atol (float) – Absolute tolerance.

Return type:

bool

raster_equal(other, strict_masked=True, warn_failure_reason=False)

Check if two rasters are equal.

This means that are equal: - The raster’s masked array’s data (including masked values), mask, fill_value and dtype, - The raster’s transform, crs and nodata values.

Parameters:

  • other (RasterType) – Other raster.

  • strict_masked (bool) – Whether to check if masked cells (in .data.mask) have the same value (in .data.data). Only when comparing two Rasters objects with masked arrays.

  • warn_failure_reason (bool) – Whether to warn for the reason of failure if the check does not pass.

Return type:

bool

reduce_at_points(points, reducer_function=None, window=None, input_latlon=False, band=None, masked=False, as_array=False, boundless=True, *, coverage=None, window_shape=None, error_structure=None, mp_config=None, resampling=None, backend='auto', resampling_kwargs=None)

Reduce raster values around point coordinates.

Samples the first band by default, or the requested band. Window reductions share the eager, Dask and multiprocessing paths of resample_at_points(), with overlap covering the selected neighborhood.

Parameters:

  • points (tuple[int | float | integer[Any] | floating[Any], int | float | integer[Any] | floating[Any]] | tuple[NDArray[floating[Any] | integer[Any]], NDArray[floating[Any] | integer[Any]]] | PointCloudBase | GeoDataFrame | DataArray) – Point(s) at which to interpolate raster value. Can be either a tuple of array-like of X/Y coordinates (same CRS as raster or latitude/longitude, see “input_latlon”) or a pointcloud in any CRS. If points fall outside of image, value returned is nan.

  • reducer_function (Callable[[NDArray[floating[Any] | integer[Any]]], float] | Reducer | None) – Legacy method argument, accepting a Reducer or reduction callable.

  • resampling (Literal['nearest', 'linear', 'bilinear', 'cubic', 'cubic_convolution', 'cubic_spline', 'lanczos', 'quintic', 'slinear', 'pchip', 'splinef2d', 'idw', 'inv_dist', 'kriging', 'mean', 'average', 'sum', 'median', 'med', 'quantile', 'q1', 'q3', 'mode', 'minimum', 'min', 'maximum', 'max', 'range', 'count', 'stdev', 'std', 'standard_deviation', 'rms', 'rmse', 'average_distance', 'average_distance_pts'] | Reducer | None) – Statistic name or an advanced Reducer. Defaults to the mean.

  • backend (Literal['auto', 'gdal', 'scipy', 'numba', 'gstools', 'gpytorch']) – “auto” or “scipy” for window reductions.

  • resampling_kwargs (dict[str, Any] | None) – Settings for a named statistic, such as a quantile q.

  • window (int | None) – Window size to read around coordinates. Must be odd.

  • input_latlon (bool) – (Only for tuple point input) Whether to convert input coordinates from latlon to raster CRS.

  • band (int | None) – Band number to extract from (from 1 to self.count).

  • masked (bool) – Whether to return a masked array, or classic array.

  • as_array (bool) – Whether to return an array of reduced values (defaults to a point cloud containing input coordinates).

  • boundless (bool) – Allow partial windows at raster edges. If False, incomplete windows return NaN.

  • coverage (Literal['center', 'all_touched', 'fractional'] | None) – Select cells by center, include all touched cells, or weight by covered area.

  • window_shape (Literal['square', 'circular'] | None) – Use a square or circular window.

  • error_structure (ErrorStructure | None) – Optional observation error model used for weighting and fitting. Use geoutils.uncertainty.propagate() to calculate output uncertainty.

  • mp_config (MultiprocConfig | None) – Multiprocessing settings for reading and reducing raster tiles.

Return type:

Any

Returns:

Point cloud of reduced values, or a scalar/array when as_array=True.

reduce_points(points, reducer_function=None, window=None, input_latlon=False, band=None, masked=False, as_array=False, boundless=True, *, coverage=None, window_shape=None, error_structure=None, mp_config=None, resampling=None, backend='auto', resampling_kwargs=None)

Reduce raster values around points using the deprecated name for reduce_at_points().

Deprecated since version 1.1: Use reduce_at_points() instead.

Return type:

Any

reproject(ref=None, crs=None, res=None, grid_size=None, bounds=None, nodata=None, dtype=None, resampling=None, tolerance=0, apply_vertical=False, transformer_options=None, force_source_nodata=None, silent=False, inplace=False, n_threads=0, memory_limit=64, mp_config=None, nodata_propagation=None, overlap_backend='auto', error_structure=None, window=None, window_shape=None, coverage=None, area_weighting='diagonal_bounds', backend='auto', resampling_kwargs=None)

Reproject raster to a different geotransform (resolution, bounds) and/or coordinate reference system (CRS).

Match-reference: a reference raster can be passed to match resolution, bounds and CRS during reprojection.

Alternatively, the destination resolution, bounds and CRS can be passed individually.

Any resampling algorithm implemented in Rasterio can be passed as a string.

Rasters referenced by ground control points (GCPs) or rational polynomial coefficients (RPCs) use their stored model instead of the affine transform. Eager, lazy Dask and multiprocessing execution all support these models with Rasterio resampling methods. GCPs use their stored CRS; RPCs use WGS84 (EPSG:4326). The result has a regular affine grid.

The reprojection can be computed out-of-memory in multiprocessing by passing a MultiprocConfig object. The reprojected raster is written to disk under the path specified in the configuration

Parameters:

  • ref (RasterType | str | None) – Reference raster to match resolution, bounds and CRS.

  • crs (CRS | str | int | None) – Destination coordinate reference system as a string or EPSG. If ref not set, defaults to this raster’s CRS.

  • res (float | Iterable[float] | None) – Destination resolution (pixel size) in units of destination CRS. Single value or (xres, yres). Do not use with grid_size.

  • grid_size (tuple[int, int] | None) – Destination grid size as (x, y). Do not use with res.

  • bounds (dict[str, float] | BoundingBox | None) – Destination bounds as a Rasterio bounding box, or a dictionary containing left, bottom, right, top bounds in the destination CRS.

  • nodata (int | float | None) – Destination nodata value. If set to None, will use the same as source. If source does not exist, will use GDAL’s default.

  • dtype (DTypeLike | None) – Destination data type of array.

  • resampling (Resampling | Literal['nearest', 'linear', 'bilinear', 'cubic', 'cubic_convolution', 'cubic_spline', 'lanczos', 'quintic', 'slinear', 'pchip', 'splinef2d', 'idw', 'inv_dist', 'kriging', 'mean', 'average', 'sum', 'median', 'med', 'quantile', 'q1', 'q3', 'mode', 'minimum', 'min', 'maximum', 'max', 'range', 'count', 'stdev', 'std', 'standard_deviation', 'rms', 'rmse', 'average_distance', 'average_distance_pts'] | Interpolator | Reducer | None) – Built-in method name, such as “linear”, “mean” or “median”. Advanced callers can supply an Interpolator or Reducer. By default, a Reducer weights source cells by their covered area. With GridNeighbours, it uses a fixed window around the transformed destination center. Defaults to bilinear. Can be configured with the global setting geoutils.config[“reprojection_method”]. See https://rasterio.readthedocs.io/en/stable/api/rasterio.enums.html#rasterio.enums.Resampling for the full list.

  • backend (Literal['auto', 'gdal', 'scipy', 'numba', 'gstools', 'gpytorch']) – “auto” preserves GDAL for supported string methods and uses SciPy for operator settings or other methods. “gdal” and “scipy” select one explicitly; “gstools” and “gpytorch” select kriging.

  • resampling_kwargs (dict[str, Any] | None) – Settings for a named method, such as a quantile q or a fitted kriging model.

  • tolerance (float) – Maximum error in source pixels when approximating a Rasterio warp transformation. Defaults to 0 for an exact transformation and consistent results across chunks. Positive values allow faster approximation but can produce different results across chunks. Requires Rasterio >= 1.5.

  • apply_vertical (bool) – Whether to apply the coordinate transformation’s vertical shift to raster values with Rasterio resampling. Use True for elevation values that need a vertical datum transformation. Applies even when the horizontal grid is unchanged. Defaults to False, so image values are only resampled. This does not control the ground height used to locate RPC pixels; use RPC_HEIGHT or RPC_DEM in transformer_options for that.

  • transformer_options (dict[str, Any] | None) – GDAL transformer options, such as RPC_DEM or RPC_HEIGHT, passed to grid calculation and reprojection. When both GCPs and RPCs are stored, select one with {"SRC_METHOD": "RPC"}, {"SRC_METHOD": "GCP_POLYNOMIAL"} or {"SRC_METHOD": "GCP_TPS"}. Explicit reproject() arguments take precedence over options with the same name.

  • force_source_nodata (int | float | None) – Force a source nodata value (read from the metadata by default).

  • inplace (bool) – (DEPRECATED. Use rast = rast.reproject() instead) Whether to reproject in-place or not.

  • silent (bool) – Whether to print warning statements.

  • n_threads (int) – Number of threads. Defaults to (os.cpu_count() - 1).

  • memory_limit (int) – Memory limit in MB for warp operations. Larger values may perform better.

  • mp_config (MultiprocConfig | None) – Configuration object containing chunk size, output file path, and an optional cluster.

  • nodata_propagation (Literal['nearest', 'ignore', 'propagate'] | None) – How invalid source values affect a custom Interpolator or Reducer. By default, an Interpolator excludes invalid neighbors and masks a destination whose nearest source cell is invalid; a Reducer excludes invalid source cells and returns nodata if none are valid. ignore excludes invalid values without the nearest-cell mask; propagate returns nodata if any contributing value is invalid. nearest explicitly requests the default Interpolator rule.

  • overlap_backend (Literal['auto', 'numba', 'exactextract', 'shapely']) – Method used to calculate a Reducer’s cell overlap. auto uses Numba for small unaligned grids when installed, then ExactExtract when installed, or Shapely. numba, exactextract and shapely select one explicitly. Aligned grids use direct overlap lengths.

  • error_structure (ErrorStructure | None) – Optional observation error model used for weighting and fitting. Use geoutils.uncertainty.propagate() to calculate output uncertainty.

  • window (int | None) – Odd source window size for a Reducer, overriding its configured neighborhood.

  • window_shape (Literal['square', 'circular'] | None) – Square or circular source window for a Reducer.

  • coverage (Literal['center', 'all_touched', 'fractional'] | None) – Select cells by center, include all touched cells, or weight by covered area. Without a fixed neighborhood, the default is fractional.

  • area_weighting (Literal['diagonal_bounds', 'intersection']) – For fractional reductions, diagonal_bounds weights the rectangle between two transformed diagonal corners. At the source boundary, its outer row or column also receives weight for the area extending beyond the raster. intersection weights only the overlap of each source cell with the transformed four-corner destination polygon. Sums use the polygon in either mode.

Return type:

RasterType | None

Returns:

Reprojected raster.

property res: tuple[float | int, float | int]

Resolution (X, Y) of the raster in georeferenced units.

resample_at_points(points, method=None, *, resampling=None, backend='auto', resampling_kwargs=None, coverage=None, band=1, input_latlon=False, as_array=False, nodata_handling=None, error_structure=None, window=None, window_shape=None, masked=False, boundless=True, **kwargs)

Resample raster values at points with a named interpolation or statistic.

Interpolation shares the SciPy and chunked paths of interp_at_points(). Statistics summarize a 3 x 3 window by default; pass window=5 to use five cells along each axis. Advanced callers can supply an Interpolator or Reducer for custom calculations or neighborhoods.

Parameters:

  • points (tuple[NDArray[floating[Any] | integer[Any]], NDArray[floating[Any] | integer[Any]]] | tuple[int | float | integer[Any] | floating[Any], int | float | integer[Any] | floating[Any]] | PointCloudBase | GeoDataFrame | DataArray) – X/Y coordinates or a point cloud defining the requested positions.

  • method (Literal['nearest', 'linear', 'bilinear', 'cubic', 'cubic_convolution', 'cubic_spline', 'lanczos', 'quintic', 'slinear', 'pchip', 'splinef2d', 'idw', 'inv_dist', 'kriging', 'mean', 'average', 'sum', 'median', 'med', 'quantile', 'q1', 'q3', 'mode', 'minimum', 'min', 'maximum', 'max', 'range', 'count', 'stdev', 'std', 'standard_deviation', 'rms', 'rmse', 'average_distance', 'average_distance_pts'] | Interpolator | Reducer | Callable[[NDArray[floating[Any] | integer[Any]]], float] | None) – Legacy alias for resampling, also accepting a reduction callable.

  • resampling (Literal['nearest', 'linear', 'bilinear', 'cubic', 'cubic_convolution', 'cubic_spline', 'lanczos', 'quintic', 'slinear', 'pchip', 'splinef2d', 'idw', 'inv_dist', 'kriging', 'mean', 'average', 'sum', 'median', 'med', 'quantile', 'q1', 'q3', 'mode', 'minimum', 'min', 'maximum', 'max', 'range', 'count', 'stdev', 'std', 'standard_deviation', 'rms', 'rmse', 'average_distance', 'average_distance_pts'] | Interpolator | Reducer | None) – Built-in method name, such as “linear”, “mean” or “median”, or an advanced operator. Defaults to the configured interpolation method.

  • backend (Literal['auto', 'gdal', 'scipy', 'numba', 'gstools', 'gpytorch']) – “auto” uses SciPy; “gstools” and “gpytorch” select a kriging solver.

  • resampling_kwargs (dict[str, Any] | None) – Method settings, such as a quantile q or a fitted kriging model.

  • coverage (Literal['center', 'all_touched', 'fractional'] | None) – Select cells by center, include all touched cells, or weight by covered area.

  • band (int) – Source band number, starting at one.

  • input_latlon (bool) – Convert tuple coordinates from longitude/latitude to the raster CRS.

  • as_array (bool) – Return values rather than a point cloud at the requested coordinates.

  • nodata_handling (Literal['nearest', 'ignore', 'propagate'] | Literal['half_order_up', 'half_order_down'] | int | None) – How missing source values affect the result. Reducers omit them by default. A distance choice also masks results near missing source cells.

  • error_structure (ErrorStructure | None) – Optional observation error model used for weighting and fitting. Use geoutils.uncertainty.propagate() to calculate output uncertainty.

  • window (int | None) – Optional odd window size for a reducer, overriding its configured neighborhood for this call.

  • window_shape (Literal['square', 'circular'] | None) – Use a square or circular window.

  • masked (bool) – Return a masked array when as_array=True.

  • boundless (bool) – Allow partial reduction windows at raster edges. If False, incomplete windows return NaN.

  • kwargs (Any) – Additional interp_at_points() options when method is an Interpolator or method name.

Return type:

Any

Returns:

Resampled values.

property rpcs: RPC | None

Rational polynomial coefficients relating image pixels to WGS84 ground coordinates.

save(filename, driver='GTiff', dtype=None, nodata=None, blank_value=None, co_opts=None, metadata=None, gcps=None, gcps_crs=None)[source]
Return type:

None

set_area_or_point(new_area_or_point, shift_area_or_point=None)

Set new pixel interpretation of the raster.

Overwrites the area_or_point attribute and updates “AREA_OR_POINT” in raster metadata tags.

Optionally, shifts the raster to correct value coordinates in relation to interpretation:

  • By half a pixel (right and downwards) if old interpretation was “Area” and new is “Point”,

  • By half a pixel (left and upwards) if old interpretration was “Point” and new is “Area”,

  • No shift for all other cases.

Parameters:

  • new_area_or_point (Literal['Area', 'Point'] | None) – New pixel interpretation “Area”, “Point” or None.

  • shift_area_or_point (bool | None) – Whether to shift with pixel interpretation, which shifts to center of pixel indexes if self.area_or_point is “Point” and maintains corner pixel indexes if it is “Area” or None. Defaults to True. Can be configured with the global setting geoutils.config[“shift_area_or_point”].

Return type:

None

Returns:

None.

set_crs(new_crs)

Set the coordinate reference system of the raster.

Return type:

None

set_mask(mask)[source]

Set a mask on the raster array.

All pixels where mask is set to True or > 0 will be masked (in addition to previously masked pixels).

Masking is performed in place. The mask must have the same shape as loaded data, unless the first dimension is 1, then it is ignored.

Parameters:

mask (NDArray[bool] | Raster) – The raster array mask.

Return type:

None

set_nodata(new_nodata, update_array=True, update_mask=True)

Set a new nodata value for all bands. This updates the old nodata into a new nodata value in the metadata, replaces the nodata values in the data of the masked array, and updates the mask of the masked array.

Careful! If the new nodata value already exists in the array, the related grid cells will be masked by default.

If the nodata value was not defined in the raster, run this function with a new nodata value corresponding to the value of nodata that exists in the data array and is not yet accounted for. All those values will be masked.

If a nodata value was correctly defined in the raster, and you wish to change it to a new value, run this function with that new value. All values having either the old or new nodata value will be masked.

If the nodata value was wrongly defined in the raster, and you wish to change it to a new value without affecting data that might have the value of the old nodata, run this function with the update_array argument as False. Only the values of the new nodata will be masked.

If you wish to set nodata value without updating the mask, run this function with the update_mask argument as False.

If None is passed as nodata, only the metadata is updated and the mask of old nodata unset.

Parameters:

  • new_nodata (int | float | None) – New nodata value.

  • update_array (bool) – Update the old nodata values into new nodata values in the data array.

  • update_mask (bool) – Update the old mask by unmasking old nodata and masking new nodata (if array is updated, old nodata are changed to new nodata and thus stay masked).

Return type:

None

set_transform(new_transform)

Set the geotransform of the raster.

Return type:

None

property shape: tuple[int, int]

Shape (i.e., height, width) of the raster in pixels.

sieve(size, connectivity=4, mask=None)

Remove connected integer regions smaller than a number of pixels.

This mirrors gdal_sieve for each band. The complete array is processed eagerly because connected regions can cross arbitrary chunk boundaries.

Parameters:

  • size (int) – Minimum connected-region size in pixels.

  • connectivity (Literal[4, 8]) – Pixel connectivity, either 4 or 8 (defaults to 4).

  • mask (RasterBase | DataArray | NDArray[bool] | None) – Optional Boolean raster or array where True values can participate in regions.

Return type:

RasterType

Returns:

Raster copy with small regions replaced by their largest neighboring value.

split_bands(bands=None, deep=True)[source]

Split the bands into separate rasters.

Parameters:

  • bands (list[int] | int | None) – Optional. A list of band indices to extract (from 1 to self.count). Defaults to all.

  • deep (bool) – Whether to return a deep or shallow copy.

Return type:

list[RasterType]

Returns:

A list of Rasters for each band.

stack(rasters, reference=0, resampling_method=None, use_ref_bounds=False)[source]

Stack this raster with one or more rasters into a multi-band raster.

All input rasters are reprojected and resampled to a common grid defined by the reference raster. The reference can be either this raster (reference=0) or another raster (reference>0, defined by its index)

The output multi-band extent is the union of all raster extents, except if use_ref_bounds is used, in which case the reference raster bounds are used. Its number of bands equals the sum of the bands from this raster and all additional rasters.

Note that all rasters will be loaded once in memory. The data is only loaded for reprojection then deleted to optimize memory usage.

Parameters:

  • rasters (Raster | list[Raster]) – Raster or list of rasters to be stacked.

  • reference (int | Raster) – Index of reference raster in the list or separate reference raster. Defaults to this raster.

  • resampling_method (str | Resampling) – Resampling method for reprojection.

  • use_ref_bounds (bool) – If True, will use reference bounds, otherwise will use maximum bounds of all rasters.

Return type:

Raster

Returns:

The merged raster with same CRS and resolution (and optionally bounds) as the reference.

stats(statistics=None, *, by=None, values=None, bins=None, categories=None, at=None, mask=None, mask_mode='inside', subsample=1, subsample_per_group=False, random_state=None, strategy='auto', backend='geoutils', subsampling_strategy='topk', interpolation='linear', align='raise', observed=True, return_masks=False, mp_config=None, fractional=False, overlap_backend='auto', error_structure=None, uncertainty_kwargs=None)

Calculate summary statistics or statistics grouped by categories, bins, or vector zones.

Omit by to summarize raster values. With by, provide categories for discrete groups, bins for continuous groups, or a vector feature column for zonal statistics:

# Global statistics
raster.stats()
raster.stats("mean")
raster.stats(["mean", "std", "nmad"])

# Categorical stats
raster.stats(["mean", "std"], by={"landcover": landcover}, categories={"landcover": classes})
# Binned stats
raster.stats(["mean", "std"], by={"elevation": elevation}, bins={"elevation": elevation_bins})
# Zonal stats
raster.stats(["mean", "std"], by={"glacier": (glacier_outlines, "id")})

# Multiple grouping
raster.stats(["mean", "std"], by={"elevation": elevation, "glacier": (glacier_outlines, "id")},
bins={"elevation": elevation_bins})
Parameters:

  • statistics (str | Callable[[Any], Any] | Reducer | Iterable[str | Callable[[Any], Any] | Reducer] | None) – Statistics to calculate (e.g. “mean”, [“mean”, “nmad”], or np.nanmedian). None returns “min”, “max”, “mean”, “median”, “std”, “nmad”, “validcount”, “totalcount” and “percentagevalidpoints”. “all” also includes “sum”, “sumofsquares”, “90thpercentile”, “iqr”, “le90” and “rmse”, plus inlier counts for masked global statistics. Grouped defaults replace “validcount” with “count”; every grouped result includes “count”.

  • by (Mapping[str, Any] | None) – Named variables to group by (e.g. {“elevation”: elevation}); use {“glacier”: (outlines, “id”)} for vector zones. Arrays must match the selected locations. Omit for global statistics.

  • values (int | Iterable[int] | Mapping[str, Any] | None) – Bands to summarize, counting from one (e.g. [1, 3]); defaults to all bands. Use a mapping for named inputs (e.g. {“elevation”: (dem, 1)}).

  • bins (Mapping[str, Any] | None) – Continuous bins keyed by grouping name (e.g. {“elevation”: 10}). Each definition is a count of equal-width bins, increasing edges (e.g. [0, 2, 5]), or a Pandas IntervalIndex to choose open/closed sides.

  • categories (Mapping[str, Iterable[Any]] | None) – Ordered categories keyed by grouping name (e.g. {“landcover”: [100, 110, 120]}). Values outside these categories are excluded.

  • at (Literal['self'] | RasterBase | DataArray | PointCloudBase | GeoDataFrame | None) – Grid or ordered point locations on which to calculate statistics (e.g. at=reference or at=”self”). Defaults to the first point input, if present, otherwise the source locations.

  • mask (RasterBase | DataArray | geoutils.vector.base.VectorBase | GeoDataFrame | ArrayLike | None) – Locations to include (True in a boolean mask, e.g. mask=dem > 1000, or features in a vector mask). Global counts describe values before this mask; “all” adds counts for values kept by the mask.

  • mask_mode (Literal['inside', 'outside']) – Keep locations “inside” or “outside” vector features; ignored for boolean masks.

  • subsample (int | float) – Fraction (e.g. 0.1 for 10%) or maximum count (e.g. 10000) of eligible locations to use. A value of 1 keeps all locations. Counts describe the sampled locations.

  • subsample_per_group (bool) – Apply subsample within each combined group (True, stratified sampling) or once across all groups (False). Without by, both use one global sample.

  • random_state (int | Generator | None) – Seed to reproduce subsampling (e.g. 42), or an existing random generator.

  • strategy (Literal['auto', 'dense', 'sparse', 'groupwise']) – Combine chunk statistics for all groups (“dense”), only groups present in each chunk (“sparse”), or gather each complete group (“groupwise”). “auto” chooses from the statistics and group count; exact quantiles and custom functions require “auto” or “groupwise” for chunked data.

  • backend (Literal['geoutils', 'flox']) – Use the GeoUtils reducer (“geoutils”) or optional Flox reducer (“flox”) for grouped statistics; see stats() for the Flox restrictions.

  • subsampling_strategy (Literal['sequential', 'topk']) – “topk” keeps the same sampled locations across chunk layouts for a fixed seed; “sequential” draws random locations using traversal order and can depend on the chunks.

  • interpolation (Literal['nearest', 'linear', 'cubic', 'quintic', 'slinear', 'pchip', 'splinef2d']) – Raster values at point locations use interp_points() with SciPy methods “nearest”, “linear”, “slinear”, “cubic”, “quintic”, “pchip” or “splinef2d”. Raster groupers listed in categories use “nearest”.

  • align (Literal['raise', 'reproject']) – “raise” rejects different grids or coordinate systems; “reproject” aligns them to the output locations. Point inputs must still share the same ordered coordinates.

  • observed (bool) – Omit declared group combinations with no eligible locations (True), or include them (False).

  • return_masks (bool) – Also return masks keyed by group labels (e.g. table, masks = raster.stats(…)). Masks cover complete groups before subsampling. Requires by.

  • mp_config (MultiprocConfig | None) – Worker and tile settings for multiprocessing, e.g. MultiprocConfig(chunks=512). Cannot be combined with Dask inputs.

  • fractional (bool) – Use the exact area covered by each polygon for one vector grouping. Overlapping groups are calculated independently, and Reducers such as Mean() and Sum() receive the covered cell fractions.

  • overlap_backend (Literal['auto', 'exactextract', 'shapely']) – Library used to calculate polygon coverage of raster cells. Auto uses ExactExtract when it is installed and the grid is compatible, and otherwise uses Shapely.

  • error_structure (ErrorStructure | None) – Optional source error model for one supplied Reducer.

  • uncertainty_kwargs (Mapping[str, Any] | None) – Options passed to geoutils.uncertainty.propagate().

Return type:

Any

Returns:

A statistic, summary dictionary, grouped dataframe, or grouped dataframe and mask mapping.

subsample(subsample, *, bands=None, mask=None, skip_nodata=True, random_state=None, as_array=False, return_indices=False, strategy='topk', force_pixel_offset='ul', force_output_to_memory=False, mp_config=None, pointcloud_return_type=None)

Randomly sample raster cells without replacement.

Overloads:

The first selected band defines eligible cells when skip_nodata=True. Its selected row and column positions are applied to every requested band, so values and point coordinates remain aligned. Auxiliary band nodata alone does not remove a point.

Raster inputs return a PointCloud by default. Xarray accessor outputs follow config[“pointcloud_return_type”], which defaults to a DataArray. Dask outputs stay lazy. With as_array=True, return values from one band or their source row and column positions.

Parameters:

  • subsample (float | int) – Fraction of eligible cells to extract when at most 1, otherwise the maximum cell count.

  • bands (int | Iterable[int] | Mapping[str, int] | None) – Bands to return. By default, point output includes every band and array output includes band one. An integer selects one band. An iterable selects bands with default b<band> column names, while a mapping gives column names as keys and band numbers as values. The first selected band controls eligible cells. Array output accepts only one integer band.

  • mask (RasterBase | DataArray | geoutils.vector.base.VectorBase | GeoDataFrame | ArrayLike | None) – With as_array=True, restrict eligible raster cells with a boolean or spatial mask.

  • skip_nodata (bool) – Whether to exclude cells where the first selected band is nodata.

  • random_state (int | Generator | None) – Random generator or seed used to make sampling reproducible.

  • as_array (bool) – Return sampled values or positions instead of a point cloud or dataframe.

  • return_indices (bool) – With as_array=True, return row and column positions instead of values.

  • strategy (Literal['sequential', 'topk']) – Use chunk-independent top-k sampling or sequential sampling.

  • force_pixel_offset (Literal['center', 'ul', 'ur', 'll', 'lr']) – Pixel position used to calculate each point coordinate.

  • force_output_to_memory (bool) – Keep the complete output in memory instead of using the automatic chunked path.

  • mp_config (MultiprocConfig | None) – Worker, tile, and output settings for multiprocessing. Point output uses GeoPackage, LAS or LAZ; LAS/LAZ stores the first selected band as elevation. Large array output uses NumPy format.

  • pointcloud_return_type (Literal['dataarray', 'dataset', 'geodataframe'] | None) – Returned representation: dataarray, dataset, or geodataframe. None uses config[“pointcloud_return_type”] for accessors, or returns a PointCloud for GeoUtils objects.

Returns:

Point output with one row per sampled cell, one-dimensional sampled values, or row/column positions.

property tags: dict[str, Any]

Metadata tags of the raster.

Returns:

Dictionary of raster metadata, potentially including sensor information.

Setter accepts:

dict[str, Any] | None

to_file(filename, driver='GTiff', dtype=None, nodata=None, blank_value=None, co_opts=None, metadata=None, gcps=None, gcps_crs=None, rpcs=None)[source]

Write the raster to file.

If blank_value is set to an integer or float, then instead of writing the contents of self.data to disk, write this provided value to every pixel instead.

Compression default value is set to ‘deflate’ (equal to GDALs: COMPRESS=DEFLATE in co_opts). Tiled default value is set to ‘NO’ as the GDAL default value. Raster is saved as a BigTIFF if the output file might exceed 4GB and as classical TIFF otherwise.

Example: dem.to_file(to_file, co_opts={‘TILED’:’YES’, ‘COMPRESS’:’LZW’})

Parameters:

  • filename (str | Path | IO[bytes]) – Filename to write the file to.

  • driver (str) – Driver to write file with.

  • dtype (DTypeLike | None) – Data type to write the image as (defaults to dtype of image data).

  • nodata (int | float | integer[Any] | floating[Any] | None) – Force a nodata value to be used (default to that of raster).

  • blank_value (int | float | None) – Use to write an image out with every pixel’s value. corresponding to this value, instead of writing the image data to disk.

  • co_opts (dict[str, str] | None) – GDAL creation options provided as a dictionary, e.g. {‘TILED’:’YES’, ‘COMPRESS’:’LZW’}.

  • metadata (dict[str, Any] | None) – Pairs of metadata to save to disk, in addition to existing metadata in self.tags.

  • gcps (list[tuple[float, ...]] | list[GroundControlPoint] | None) – GroundControlPoint objects or [row, col, x, y, (z)] tuples. Defaults to stored points; an empty list omits them.

  • gcps_crs (CRS | None) – CRS of the GCPs. Defaults to the stored GCP CRS, then the raster CRS.

  • rpcs (RPC | None) – Rational polynomial coefficients. Defaults to stored coefficients.

Return type:

None

Returns:

None.

to_nanarray(floating_dtype='float32', *, return_mask=False)

Get NaN array from the raster.

Overloads:

Optionally, return the mask from the masked array.

Parameters:

  • floating_dtype (DTypeLike) – Floating dtype to convert to, if masked array is not of floating type.

  • return_mask (bool) – Whether to return the mask of valid data.

:returns Array with masked data as NaNs, (Optional) Mask of invalid data.

to_pointcloud(data_name='b1', data_band=1, auxiliary_data_bands=None, auxiliary_column_names=None, subsample=1, skip_nodata=True, as_array=False, random_state=None, force_pixel_offset='ul', mp_config=None, force_output_to_memory=False, pointcloud_return_type=None)

Convert raster cells to a point cloud, optionally selecting a random subsample.

Overloads:

Cell selection and output construction use the same implementation as subsample(). With the default subsample=1, every eligible cell is returned, while this method keeps conversion arguments first.

Parameters:

  • data_name (str) – Name to use for point cloud data attribute, defaults to “bX” where X is the data band number.

  • data_band (int) – (Only for multi-band rasters) Band to use for data attribute, defaults to first. Band counting starts at 1.

  • auxiliary_data_bands (Iterable[int] | None) – (Only for multi-band rasters) Other bands to save as auxiliary data attributes, defaulting to every band other than data_band. Pass an empty iterable to keep only the main band.

  • auxiliary_column_names (Iterable[str] | None) – (Only for multi-band rasters) Names to use for auxiliary data bands, only if auxiliary data bands is not none, defaults to “b1”, “b2”, etc.

  • subsample (float | int) – Fraction or maximum number of eligible cells to return, with 1 selecting every cell.

  • skip_nodata (bool) – Whether to skip nodata values.

  • as_array (bool) – Return an array instead of a vector.

  • random_state (int | Generator | None) – Random state or seed number.

  • force_pixel_offset (Literal['center', 'ul', 'ur', 'll', 'lr']) – Force offset to derive point coordinate with. Raster coordinates normally only associate to upper-left corner “ul” (“Area” definition) or center (“Point” definition).

  • mp_config (MultiprocConfig | None) – Worker, tile and output settings for multiprocessing. Point output uses GeoPackage, LAS or LAZ; LAS/LAZ stores the main raster value as elevation. Array output does not write a file. Cannot be combined with a Dask source.

  • force_output_to_memory (bool) – Keep the complete output in memory instead of using the automatic chunked path.

  • pointcloud_return_type (Literal['dataarray', 'dataset', 'geodataframe'] | None) – Returned representation: dataarray, dataset, or geodataframe. None uses config[“pointcloud_return_type”] for accessors, or returns a PointCloud for GeoUtils objects.

Raises:

ValueError – If the sample count or fraction is poorly formatted.

Returns:

A point cloud, or array of the shape (N, 2 + count) where N is the sample count.

to_points(**kwargs)
to_rio_dataset()[source]

Export to a Rasterio in-memory dataset.

Return type:

DatasetReader

to_xarray(name=None)[source]

Convert a xarray.DataArray raster.

Integer-type rasters are cast to float32.

Parameters:

name (str | None) – Name attribute for the data array.

Return type:

DataArray

Returns:

Data array.

property transform: Affine

Geotransform of the raster.

Returns:

Affine matrix geotransform.

translate(xoff, yoff, distance_unit='georeferenced', inplace=False)

Translate a raster by a (x,y) offset.

Overloads:

The translation only updates the geotransform (no resampling is performed).

GCP/RPC referencing requires reproject() to an affine grid first.

Parameters:

  • xoff (float) – Translation x offset.

  • yoff (float) – Translation y offset.

  • distance_unit (Literal['georeferenced', 'pixel']) – Distance unit, either ‘georeferenced’ (default) or ‘pixel’.

  • inplace (bool) – Whether to modify the raster in-place.

Returns:

Translated raster (or None if inplace).

variogram(*, band=1, n_pairs=1000000, sampling='loglag', estimator='dowd', bins='log', n_lags=24, min_lag=None, max_lag=None, n_runs=1, model=None, fit_kwargs=None, random_state=None, mask=None, **pair_sampling_kwargs)

Estimate a lightweight empirical variogram from raster pairs.

Parameters:

  • band (int) – Raster band to sample, counting from 1.

  • n_pairs (int) – Number of finite pairs targeted in each run (e.g. 100_000).

  • sampling (Literal['loglag', 'random_xy']) – How to select pairs: "loglag" balances short and long distances, while "random_xy" selects endpoints independently.

  • estimator (str | Callable[[NDArray[floating[Any] | integer[Any]]], float]) – Semivariance estimator from SciKit-GStat: "dowd", "matheron", "cressie", "genton", "minmax", "entropy" or "percentile". A function can instead map absolute pair differences to one value per distance bin.

  • bins (Literal['log', 'uniform'] | Iterable[float]) – Distance bins: "log" for logarithmic spacing, "uniform" for equal widths, or explicit edges (e.g. [1, 10, 100]).

  • n_lags (int) – Number of distance bins when bins is "log" or "uniform".

  • min_lag (float | None) – Minimum sampled distance in CRS units; defaults to the smaller pixel dimension.

  • max_lag (float | None) – Maximum sampled distance in CRS units; defaults to the diagonal between outermost cell centers.

  • n_runs (int) – Independent samples to average; repeat sampling to estimate each distance bin’s standard error.

  • model (str | Callable[..., Any] | list[str | Callable[..., Any]] | None) – SciKit-GStat model to fit: "spherical", "exponential", "gaussian", "cubic", "stable" or "matern", or the corresponding model function. Sum a list of models ordered from short to long range (e.g. [“spherical”, “exponential”]). None keeps only the empirical variogram.

  • fit_kwargs (dict[str, Any] | None) – Options for Variogram.fit(): use_nugget, bounds, p0 or maxfev (e.g. {“use_nugget”: True}); optimization uses SciPy curve_fit().

  • random_state (int | Generator | None) – Seed or NumPy generator for reproducible sampling across runs (e.g. 42).

  • mask (RasterBase | DataArray | geoutils.vector.base.VectorBase | GeoDataFrame | ArrayLike | None) – Cells to keep: True values in a boolean mask or aligned mask raster, or cells inside vector geometries.

  • pair_sampling_kwargs (Any) – Extra pairsample() options (e.g. strategy or max_rounds).

Return type:

Variogram

Returns:

Variogram with distance bins, pair counts and semivariance, plus sampling errors and a fitted model when requested.

property width: int

Width of the raster in pixels.

xy2ij(x, y, op=<class 'numpy.float32'>, precision=None, shift_area_or_point=None)

Get indexes (row,column) of coordinates (x,y).

By default, the indexes are shifted with the interpretation of pixel coordinates “AREA_OR_POINT” of the raster, to ensure that the indexes of points represent the right location. See parameter description of shift_area_or_point for more details.

This function is reversible with ij2xy for any pixel interpretation.

Parameters:

  • x (ArrayLike) – X coordinates.

  • y (ArrayLike) – Y coordinates.

  • op (type) – Operator to compute index.

  • precision (float | None) – Precision passed to rasterio.transform.rowcol().

  • shift_area_or_point (bool | None) – Whether to shift with pixel interpretation, which shifts to center of pixel indexes if self.area_or_point is “Point” and maintains corner pixel indexes if it is “Area” or None. Defaults to True. Can be configured with the global setting geoutils.config[“shift_area_or_point”].

Returns i, j:

Indices of (x,y) in the image.

Return type:

tuple[NDArray[floating[Any] | integer[Any]], NDArray[floating[Any] | integer[Any]]]