Raster algebra & apply¶

Run computations over raster values:

• apply(func) — apply a NumPy function to a band, returning a new Dataset.
• map_blocks(func, tile_size=...) — apply a function tile-by-tile (memory-friendly).
• overlay(classes_map) — group a raster's values by the classes in a second raster.

(For lazy, Dask-backed algebra see the Dask quickstart.)

Setup¶

In [1]:

%matplotlib inline

import tempfile
from pathlib import Path

import numpy as np

DATA = Path('../../../examples/data')
WORK = Path(tempfile.mkdtemp(prefix='pyramids-t3-'))
DATA.is_dir(), WORK.is_dir()

%matplotlib inline import tempfile from pathlib import Path import numpy as np DATA = Path('../../../examples/data') WORK = Path(tempfile.mkdtemp(prefix='pyramids-t3-')) DATA.is_dir(), WORK.is_dir()

Out[1]:

(True, True)

In [2]:

from pyramids.dataset import Dataset

ds = Dataset.read_file(DATA / 'acc4000.tif')
ds.shape

from pyramids.dataset import Dataset ds = Dataset.read_file(DATA / 'acc4000.tif') ds.shape

2026-07-11 14:42:04 | INFO | pyramids.base.config | Logging is configured.

Out[2]:

(1, 13, 14)

The source raster — flow accumulation values that the algebra below operates on.

In [3]:

ds.plot(band=0, title="Source (acc4000)", cmap="viridis")

ds.plot(band=0, title="Source (acc4000)", cmap="viridis")

Out[3]:

<cleopatra.array_glyph.ArrayGlyph at 0x7f2163d58ad0>

/home/runner/work/pyramids/pyramids/.pixi/envs/docs/lib/python3.14/site-packages/matplotlib/colors.py:816: RuntimeWarning: overflow encountered in multiply
  xa *= self.N

Element-wise — apply¶

Apply a function to every cell of a band; returns a new Dataset.

In [4]:

doubled = ds.apply(lambda a: a * 2)
type(doubled).__name__, float(
    np.nanmax(
        np.where(
            doubled.read_array() == doubled.no_data_value[0],
            np.nan,
            doubled.read_array(),
        )
    )
)

doubled = ds.apply(lambda a: a * 2) type(doubled).__name__, float( np.nanmax( np.where( doubled.read_array() == doubled.no_data_value[0], np.nan, doubled.read_array(), ) ) )

Out[4]:

('Dataset', 176.0)

The apply result — every cell doubled. Same spatial pattern, values scaled by 2.

In [5]:

doubled.plot(band=0, title="apply: a * 2", cmap="viridis")

doubled.plot(band=0, title="apply: a * 2", cmap="viridis")

Out[5]:

<cleopatra.array_glyph.ArrayGlyph at 0x7f2163ec6350>

/home/runner/work/pyramids/pyramids/.pixi/envs/docs/lib/python3.14/site-packages/matplotlib/colors.py:816: RuntimeWarning: overflow encountered in multiply
  xa *= self.N

Tile-by-tile — map_blocks¶

Process the raster in tile_size blocks — the same result, a bounded memory footprint.

In [6]:

sen = Dataset.read_file(DATA / 'geotiff' / 'sentinel_crop.tif')
incremented = sen.map_blocks(lambda a: a + 1, tile_size=128)
type(incremented).__name__, incremented.shape

sen = Dataset.read_file(DATA / 'geotiff' / 'sentinel_crop.tif') incremented = sen.map_blocks(lambda a: a + 1, tile_size=128) type(incremented).__name__, incremented.shape

Out[6]:

('Dataset', (1, 256, 256))

The tile-by-tile map_blocks output — the Sentinel crop with every value incremented by 1.

In [7]:

incremented.plot(band=0, title="map_blocks: a + 1", cmap="viridis")

incremented.plot(band=0, title="map_blocks: a + 1", cmap="viridis")

Out[7]:

<cleopatra.array_glyph.ArrayGlyph at 0x7f2163b1cf50>

Group by class — overlay¶

overlay reads a second raster of class labels and returns a dict mapping each class value to the list of data-raster values that fall in it — a raster-on-raster zonal summary.

In [8]:

classes = Dataset.read_file(DATA / 'geotiff' / 'sentinel-classes.tif')
grouped = sen.overlay(classes)
type(grouped).__name__, sorted(grouped), [len(v) for v in grouped.values()]

classes = Dataset.read_file(DATA / 'geotiff' / 'sentinel-classes.tif') grouped = sen.overlay(classes) type(grouped).__name__, sorted(grouped), [len(v) for v in grouped.values()]

Out[8]:

('dict', [np.uint16(2)], [5771])

Notes¶

• apply(..., inplace=True) mutates the dataset; overlay takes exclude_value= to drop a sentinel before grouping.
• See also: Zonal statistics (the vector-polygon equivalent).