Windowed & tiled reads¶

Read only the pixels you need instead of the whole raster — the basis of out-of-core and tile-server workflows:

• read_array(window=[xoff, yoff, xsize, ysize]) — a pixel window.
• read_part(bbox, bbox_crs=...) — a geographic bounding-box subset.
• get_tile(size) — iterate fixed-size tiles across the raster.
• get_block_arrangement() — the raster's native on-disk block layout.

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 / 'geotiff' / 'sentinel-classes.tif')
ds.shape, ds.epsg

from pyramids.dataset import Dataset ds = Dataset.read_file(DATA / 'geotiff' / 'sentinel-classes.tif') ds.shape, ds.epsg

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

Out[2]:

((1, 256, 256), 4326)

The full 256x256 raster — each windowed read below pulls a sub-region of this image.

In [3]:

ds.plot(band=0, title="full raster — Sentinel-2 classes (256x256)")

ds.plot(band=0, title="full raster — Sentinel-2 classes (256x256)")

Out[3]:

<cleopatra.array_glyph.ArrayGlyph at 0x7fe27894cad0>

A pixel window — read_array(window=...)¶

[xoff, yoff, xsize, ysize] in pixels — reads only that block.

In [4]:

patch = ds.read_array(window=[50, 50, 100, 100])
np.asarray(patch).shape

patch = ds.read_array(window=[50, 50, 100, 100]) np.asarray(patch).shape

Out[4]:

(100, 100)

A geographic subset — read_part¶

Give a bounding box (with its CRS); the reader clips and reads just that area.

In [5]:

xmin, ymin, xmax, ymax = ds.bbox
dx, dy = (xmax - xmin) / 4, (ymax - ymin) / 4
part = ds.read_part((xmin + dx, ymin + dy, xmax - dx, ymax - dy), bbox_crs=ds.epsg)
np.asarray(part).shape

xmin, ymin, xmax, ymax = ds.bbox dx, dy = (xmax - xmin) / 4, (ymax - ymin) / 4 part = ds.read_part((xmin + dx, ymin + dy, xmax - dx, ymax - dy), bbox_crs=ds.epsg) np.asarray(part).shape

Out[5]:

(129, 129)

read_part returns a bare array; cropping ds to the same inner bbox gives a Dataset we can plot, so you can see exactly which sub-region the window covers against the full raster above.

In [6]:

subset = ds.crop(bbox=(xmin + dx, ymin + dy, xmax - dx, ymax - dy), epsg=ds.epsg)
subset.plot(band=0, title="inner-bbox subset")

subset = ds.crop(bbox=(xmin + dx, ymin + dy, xmax - dx, ymax - dy), epsg=ds.epsg) subset.plot(band=0, title="inner-bbox subset")

Out[6]:

<cleopatra.array_glyph.ArrayGlyph at 0x7fe278a92350>

Iterate tiles — get_tile¶

Yields fixed-size windows, so you can stream a large raster tile by tile.

In [7]:

tiles = list(ds.get_tile(size=128))
len(tiles), tiles[0].shape

tiles = list(ds.get_tile(size=128)) len(tiles), tiles[0].shape

Out[7]:

(4, (128, 128))

Native block layout — get_block_arrangement¶

The blocks GDAL actually stores on disk — read along these for the least I/O.

In [8]:

blocks = ds.get_block_arrangement()
list(blocks.columns), len(blocks)

blocks = ds.get_block_arrangement() list(blocks.columns), len(blocks)

Out[8]:

(['x_offset', 'y_offset', 'window_xsize', 'window_ysize'], 16)

Notes¶

• For lazy, parallel block processing use read_array(chunks=...) — see the Dask quickstart — Dataset.
• See also: Raster algebra (map_blocks applies a function per tile).