Rasters: grids, transforms, bands and nodata

Move the lesson’s 10 m image half a pixel, and 750 of 1,000 points sampled from it read a neighbouring pixel, without a single error message. A raster file is a grid of numbers plus a few lines of metadata: six numbers that put every pixel on the ground, a coordinate reference system that says which ground, the bands, the data type, and the value that means “no data”. One wrong line, and every result built on the grid is wrong in the same quiet way.

Here you take apart the metadata of one raster, the two-band image of 800 by 600 pixels that this module’s lessons share. You will move between pixels and map coordinates, choose a data type and a nodata value that cannot clash with real data, work out which blocks a read of part of a raster fetches, and choose how a GeoTIFF you write is tiled and compressed.

Six numbers turn a row and a column into a place

The lesson raster is 800 pixels wide and 600 tall, each pixel 10 m on a side, in UTM zone 43N near Pune. Its affine transformaffine transform: Six numbers, a to f, that turn a raster pixel's column and row into map coordinates: x = a·col + b·row + c and y = d·col + e·row + f give the pixel's top-left corner. rasterio keeps them in this order; GDAL's geotransform holds the same six as c, a, b, f, d, e. is six numbers, which rasterio, the Python library for rasters built on GDAL, keeps in the order a, b, c, d, e, f, and here they are 10, 0, 370000, 0, −10 and 2046000 (rasterio, Georeferencing). A pixel’s column and row become map coordinates with two sums:

In words: start at the raster’s top-left corner, (c, f) = (370000, 2046000), and move a = 10 m east for every column and e = −10 m, 10 m south, for every row. b and d are 0 because the image is north-up. For row 254, column 432, x = 10 × 432 + 370000 = 374,320 and y = −10 × 254 + 2046000 = 2,043,460, the top-left corner of that pixel. Look at where the corner and the centres of the first pixels land:

The sums give a pixel’s corner; its centre is half a pixel further in. GDAL counts from the top-left corner of the top-left pixel, so that pixel’s centre is at (0.5, 0.5) in pixel units (GDAL, Raster Data Model). rasterio’s xy(row, col) returns the centre unless its offset asks for a corner, and its rowcol(x, y) rounds down with floor (rasterio, transform). The programs of this lesson share one module that builds the lesson raster from formulas, so they run anywhere without a file:

"""The lesson's raster: two bands, red and near-infrared, of 800 by 600 pixels of 10 m in UTM zone 43N.

Each value is a reflectance times 10,000, stored as an unsigned 16-bit integer, and 0 marks a pixel without data.
The scene is made by whole-number formulas, so every program that imports it, here or in a browser, gets the
same numbers: a river, fields of crops and bare soil, a town, and a hole of 2,400 pixels with no data.
"""

WIDTH, HEIGHT = 800, 600  # columns, rows
TRANSFORM = (10.0, 0.0, 370_000.0, 0.0, -10.0, 2_046_000.0)  # a, b, c, d, e, f
CRS = "EPSG:32643"  # WGS 84 / UTM zone 43N
NODATA = 0


def noise(row, col):
    """A fixed whole number from -40 to 40 for each row and column, the same on every computer."""
    h = ((row * 73_856_093) ^ (col * 19_349_663)) & 0xFFFFFFFF
    h = ((h ^ (h >> 13)) * 1_274_126_177) & 0xFFFFFFFF
    return h % 81 - 40


def wobble(row, col):
    """Texture that changes smoothly over 16 pixels, as real ground does, plus a little noise of its own."""
    r, c, fr, fc = row // 16, col // 16, row % 16, col % 16
    top = noise(r, c) * (16 - fc) + noise(r, c + 1) * fc
    bottom = noise(r + 1, c) * (16 - fc) + noise(r + 1, c + 1) * fc
    return (top * (16 - fr) + bottom * fr) // 256 + noise(row, col) // 8


def cover(row, col):
    """What lies under a pixel: water, town, soil, crop, or None in the hole without data."""
    if 420 <= row < 460 and 100 <= col < 160:
        return None
    if abs(row - (260 + col * col // 4000)) < 18:
        return "water"
    if 60 <= row < 200 and 520 <= col < 700:
        return "town"
    field_row, field_col = row // 80, col // 100
    if (field_row * 3 + field_col) % 5 == 0:
        return "soil"
    return "crop"


def pixel(row, col):
    """(red, nir) of one pixel, reflectance times 10,000."""
    kind = cover(row, col)
    if kind is None:
        return NODATA, NODATA
    w = wobble(row, col)
    if kind == "water":
        return 350 + w, 200 + w // 2
    if kind == "town":
        return 1_700 + w, 2_200 + w
    if kind == "soil":
        return 2_200 + w, 2_800 + w
    vigour = ((row // 80) * 7 + (col // 100) * 13) % 10  # 0 to 9, one value per field
    return 800 - vigour * 40 + w, 3_000 + vigour * 150 + (col % 100) * 3 + w


def bands():
    """The two bands as lists of rows: red, then near-infrared."""
    red, nir = [], []
    for row in range(HEIGHT):
        pairs = [pixel(row, col) for col in range(WIDTH)]
        red.append([p[0] for p in pairs])
        nir.append([p[1] for p in pairs])
    return red, nir

The first program does both conversions by hand:

"""The lesson raster's affine transform: from a row and a column to map coordinates, and back."""
import math

from scene import HEIGHT, TRANSFORM, WIDTH


def xy(t, row, col, offset="center"):
    """Map x, y of a pixel's centre, or of its upper-left corner with offset="ul"."""
    a, b, c, d, e, f = t
    s = 0.5 if offset == "center" else 0.0
    return a * (col + s) + b * (row + s) + c, d * (col + s) + e * (row + s) + f


def rowcol(t, x, y):
    """The row and column of the pixel that holds x, y in a north-up raster (b = d = 0)."""
    a, _, c, _, e, f = t
    return math.floor((y - f) / e), math.floor((x - c) / a)


def show(p):
    return f"({p[0]:,.1f}, {p[1]:,.1f})"


a, b, c, d, e, f = TRANSFORM
print(f"rasterio order (a, b, c, d, e, f): {TRANSFORM}")
print(f"GDAL order (c, a, b, f, d, e):     {(c, a, b, f, d, e)}")
for row, col in [(0, 0), (HEIGHT - 1, WIDTH - 1)]:
    print(f"row {row}, col {col}: corner {show(xy(TRANSFORM, row, col, 'ul'))}, centre {show(xy(TRANSFORM, row, col))}")

west, north = xy(TRANSFORM, 0, 0, "ul")
east, south = xy(TRANSFORM, HEIGHT, WIDTH, "ul")  # the far corner of the last pixel
print(f"bounds: west {west:,.0f}, south {south:,.0f}, east {east:,.0f}, north {north:,.0f}")
print(f"on the ground: {(east - west) / 1000:g} km by {(north - south) / 1000:g} km")

clinic = (374_321.4, 2_043_456.8)
row, col = rowcol(TRANSFORM, *clinic)
print(f"the clinic {show(clinic)} is in row {row}, col {col}, centre {show(xy(TRANSFORM, row, col))}")

gdal_order = (c, a, b, f, d, e)  # handed on as if it were a, b, c, d, e, f
for row, col in [(0, 0), (1, 2)]:
    print(f"GDAL's order read as rasterio's: row {row}, col {col} lands at {show(xy(gdal_order, row, col, 'ul'))}")

rasterio order (a, b, c, d, e, f): (10.0, 0.0, 370000.0, 0.0, -10.0, 2046000.0)
GDAL order (c, a, b, f, d, e):     (370000.0, 10.0, 0.0, 2046000.0, 0.0, -10.0)
row 0, col 0: corner (370,000.0, 2,046,000.0), centre (370,005.0, 2,045,995.0)
row 599, col 799: corner (377,990.0, 2,040,010.0), centre (377,995.0, 2,040,005.0)
bounds: west 370,000, south 2,040,000, east 378,000, north 2,046,000
on the ground: 8 km by 6 km
the clinic (374,321.4, 2,043,456.8) is in row 254, col 432, centre (374,325.0, 2,043,455.0)
GDAL's order read as rasterio's: row 0, col 0 lands at (0.0, -10.0)
GDAL's order read as rasterio's: row 1, col 2 lands at (740,010.0, 4,091,990.0)

The far corner of the last pixel is column 800, row 600, so the raster spans 370,000 to 378,000 E and 2,040,000 to 2,046,000 N: 8 km by 6 km. Going back from the clinic at 374,321.4 E, 2,043,456.8 N, subtract the corner, divide by the pixel size and round down: row 254, column 432, the pixel whose centre is 374,325, 2,043,455.

GDAL stores the same six numbers in another order, c, a, b, f, d, e: the geotransformgeotransform: Six numbers that place a raster on the ground: the x of its top-left corner, the cell width, a row rotation, the y of that corner, a column rotation and the cell height. More in “Vector and raster data models” of the earlier lesson. The last lines of the output show what the wrong order does. Read as rasterio’s, GDAL’s tuple puts the first pixel at (0, −10) and the pixel in row 1, column 2 at 740,010, 4,091,990, thousands of kilometres away.

Mixing up the corner and the centre moves a whole raster. If a program takes the first pixel’s centre for the raster’s corner, every pixel shifts 5 m east and 5 m south. This program samples 1,000 random points from the raster and from the shifted copy:

"""A raster placed half a pixel off: how many of 1,000 sample points then read a neighbouring pixel."""
import math
import random

from scene import HEIGHT, TRANSFORM, WIDTH


def rowcol(t, x, y):
    a, _, c, _, e, f = t
    return math.floor((y - f) / e), math.floor((x - c) / a)


a, b, c, d, e, f = TRANSFORM
# The first pixel's centre, (370005, 2045995), taken for its corner: the whole grid moves 5 m east and 5 m south.
shifted = (a, b, c + a / 2, d, e, f + e / 2)

rng = random.Random(7)
points = [(rng.uniform(c + 20, c + a * WIDTH - 20), rng.uniform(f + e * HEIGHT + 20, f - 20)) for _ in range(1_000)]
moved = [p for p in points if rowcol(TRANSFORM, *p) != rowcol(shifted, *p)]
both = [p for p in moved if rowcol(TRANSFORM, *p)[0] != rowcol(shifted, *p)[0] and rowcol(TRANSFORM, *p)[1] != rowcol(shifted, *p)[1]]
print(f"shift of the grid: {math.hypot(a / 2, e / 2):.2f} m (half a pixel east and half a pixel south)")
print(f"points that read another pixel: {len(moved)} of {len(points)}")
print(f"  of them, a diagonal neighbour: {len(both)}")
print(f"points that still read the right pixel: {len(points) - len(moved)}")

shift of the grid: 7.07 m (half a pixel east and half a pixel south)
points that read another pixel: 750 of 1000
  of them, a diagonal neighbour: 236
points that still read the right pixel: 250

750 of the 1,000 points read another pixel, 236 of them a diagonal neighbour, and nothing failed. GDAL records how a file’s values were sampled in its AREA_OR_POINT metadata, “Area” by default or “Point” for a value measured at the pixel’s centre, but its georeferencing stays corner-based either way; GeoTIFF stores the same choice as PixelIsArea or PixelIsPoint (OGC GeoTIFF 1.1).

The metadata says what every number means

The transform places the grid; the rest of the metadata says what the grid holds. rasterio gathers it in a profile, a dictionary of the settings a file was written with, which a new file can be written from (rasterio, Profiles). The program prints the lesson raster’s profile and checks a point against it:

"""What the lesson raster's metadata says, and a check that a point and the raster share one CRS."""
from scene import CRS, HEIGHT, NODATA, TRANSFORM, WIDTH, pixel

profile = {
    "driver": "GTiff",
    "width": WIDTH,
    "height": HEIGHT,
    "count": 2,
    "dtype": "uint16",
    "nodata": NODATA,
    "crs": CRS,
    "transform": TRANSFORM,
}
for key, value in profile.items():
    print(f"{key:>9}: {value}")

a, _, c, _, e, f = TRANSFORM
west, east, north, south = c, c + a * WIDTH, f, f + e * HEIGHT
print(f"resolution: {a:g} m by {-e:g} m; bounds {west:,.0f} to {east:,.0f} E, {south:,.0f} to {north:,.0f} N")
print("band 1 is red, band 2 is near-infrared; row 254, col 432 holds", pixel(254, 432))


def inside(x, y):
    return west <= x < east and south < y <= north


for name, (x, y) in [("clinic in UTM 43N metres", (374_321.4, 2_043_456.8)), ("same clinic in degrees", (73.80964, 18.47775))]:
    print(f"{name}: {x}, {y} -> {'inside the raster' if inside(x, y) else 'outside: transform the point first'}")

   driver: GTiff
    width: 800
   height: 600
    count: 2
    dtype: uint16
   nodata: 0
      crs: EPSG:32643
transform: (10.0, 0.0, 370000.0, 0.0, -10.0, 2046000.0)
resolution: 10 m by 10 m; bounds 370,000 to 378,000 E, 2,040,000 to 2,046,000 N
band 1 is red, band 2 is near-infrared; row 254, col 432 holds (660, 3526)
clinic in UTM 43N metres: 374321.4, 2043456.8 -> inside the raster
same clinic in degrees: 73.80964, 18.47775 -> outside: transform the point first

• driver: GTiff. The file is a GeoTIFFGeoTIFF: A TIFF image file with tags that place it on the Earth: its coordinate reference system and the position and size of its pixels. GDAL reads and writes it with its GTiff driver, in strips by default or in tiles., read and written by GDAL’s GTiff driverGDAL driver: The part of GDAL that reads or writes one format, such as GPKG, GeoJSON or GTiff. ogrinfo --formats and gdalinfo --formats list the drivers a copy of GDAL has. More in “GDAL and OGR on the command line”.
• width, height and count. 800 columns, 600 rows and 2 bandsband: One layer of values in a raster, such as the red, green or blue channel of an image, or the single layer of an elevation model. More in “Vector and raster data models”. A band number is only a position: band 1 is red and band 2 near-infrared because the lesson says so. Sentinel-2’s own files name theirs, B4 for red at about 665 nm and B8 for near-infrared at about 833 nm (SentiWiki, S2 Mission).
• dtype and nodata. Unsigned 16-bit integers, with 0 marking a pixel without data; the next section is about both.
• crs and transform. EPSG:32643, WGS 84 / UTM zone 43N, so x and y are metres east and north in that zone, and the transform’s a and −e give the resolution, 10 m by 10 m.

The CRScoordinate reference system (CRS): The definition that gives coordinates a meaning on the Earth: what the numbers measure, in which units and from which origin, including the datum and, for flat maps, the map projection. More in “What GIS is and the questions it answers” is what turns the numbers into places, and it is also the commonest reason a point misses its pixel. The clinic given as 73.80964 E, 18.47775 N fails the bounds check, because 73.8 is not an easting between 370,000 and 378,000. Transform such points into the raster’s CRS before you look them up. When a whole vector layer and a raster meet, transform the vector layer: its points and lines are moved vertex by vertex and keep their values, while a reprojected raster must be resampled, every pixel’s value estimated again, which the lesson on resampling measures.

A data type sets the range, and nodata must lie outside it

Every pixel of a band has the same data typedata type (dtype): The kind of number every pixel of a band holds, such as uint8 (whole numbers 0 to 255), uint16 (0 to 65,535), int16 (-32,768 to 32,767) or float32. It fixes the range of values, their precision and the bytes each pixel takes., and the type decides three things: which values fit, how exact they are and how many bytes each pixel takes. The program lists the common ones with the cost of one band of the lesson raster in each, and then tries two nodata values on a coastal elevation model:

"""What each raster data type can hold, what a band of the lesson raster costs in it, and a nodata value gone wrong."""
import struct

from scene import HEIGHT, WIDTH

TYPES = [  # name, bytes, smallest and largest value
    ("uint8", 1, 0, 2**8 - 1),
    ("uint16", 2, 0, 2**16 - 1),
    ("int16", 2, -(2**15), 2**15 - 1),
    ("int32", 4, -(2**31), 2**31 - 1),
]
for name, size, low, high in TYPES:
    print(f"{name:>7}: {size} byte(s), {WIDTH * HEIGHT * size:>9,} bytes a band, {low:,} to {high:,}")
print(f"{'float32':>7}: 4 byte(s), {WIDTH * HEIGHT * 4:>9,} bytes a band, about 7 significant digits")


def as_float32(x):
    return struct.unpack("<f", struct.pack("<f", x))[0]


for value in (0.4164, 2_046_000.37):
    print(f"float32 stores {value} as {as_float32(value)!r}")

# A coastal elevation model in whole metres (int16): the beach and the mudflats are at 0 m.
dem = [[max(-6, (col - 30) // 4 + (row % 7 == 0)) for col in range(100)] for row in range(100)]
zeros = sum(row.count(0) for row in dem)
print(f"coastal DEM: {zeros:,} of 10,000 cells are at exactly 0 m")
for nodata in (0, -(2**15)):
    hidden = sum(row.count(nodata) for row in dem)
    print(f"  with nodata = {nodata}: {hidden:,} real cells would read as missing")

  uint8: 1 byte(s),   480,000 bytes a band, 0 to 255
 uint16: 2 byte(s),   960,000 bytes a band, 0 to 65,535
  int16: 2 byte(s),   960,000 bytes a band, -32,768 to 32,767
  int32: 4 byte(s), 1,920,000 bytes a band, -2,147,483,648 to 2,147,483,647
float32: 4 byte(s), 1,920,000 bytes a band, about 7 significant digits
float32 stores 0.4164 as 0.4163999855518341
float32 stores 2046000.37 as 2046000.375
coastal DEM: 400 of 10,000 cells are at exactly 0 m
  with nodata = 0: 400 real cells would read as missing
  with nodata = -32768: 0 real cells would read as missing

• uint8 holds 0 to 255 in one byte: classes, masks, photographs.
• uint16 holds 0 to 65,535 in two bytes. The lesson raster stores reflectance times 10,000 in it, as Sentinel-2’s products store their digital numbers in 16-bit integers with a quantification value of 10,000 (SentiWiki, S2 Processing).
• int16 holds −32,768 to 32,767: elevations in metres, which can be below sea level.
• float32 takes four bytes and keeps about seven significant digits. That suits an index: 0.4164 is stored as 0.41639998. It does not suit large values with decimals: 2,046,000.37 comes back as 2,046,000.375.

Integer types have a second trap: a value that does not fit wraps round instead of failing. NumPy “follows C casting rules”, so 300 converted to an 8-bit integer becomes 44, and arithmetic past a type’s range gives a wrong result, not an error (NumPy, Data types). The next lesson shows what that does to a vegetation index.

The nodata valuenodata value: A value, such as -9999, that marks a raster cell with no measurement. It is a code, not data, and must be left out of sums, means and counts. More in “Vector and raster data models” works only if real data can never take it. In the coastal model, 400 of its 10,000 cells, the beach and the mudflats, are at exactly 0 m: with nodata = 0, all 400 would read as missing, and the shoreline would vanish from every statistic. The int16 minimum, −32,768, is an elevation no place has, so it hides nothing. Sentinel-2 chose the same way: its digital number 0 is reserved for no data, and valid values run from 1 up (SentiWiki, S2 Products).

A windowed read fetches whole blocks

A band of the lesson raster is 960,000 bytes, and one 10 m band of a Sentinel-2 tile, about 11,000 pixels square, is about 242 MB uncompressed. Most work needs only a part: a field, a ward, the square round a clinic. A windowwindow: A rectangle of a raster's pixels, given by its first column, its first row, its width and its height, that a program reads or writes instead of the whole raster. names that part in pixels, by its column offset, its row offset, its width and its height (rasterio, Windowed reading). Reading it does not read only those pixels, because a file is stored in blocksblock (tile or strip): The unit in which a raster file stores its pixels: a tile, such as 256 by 256 pixels, or a strip of whole rows. A reader fetches every block that a window overlaps, whole.. GDAL writes a GeoTIFF in strips of whole rows unless told otherwise, each strip 8 KB or less, and with TILED=YES in tiles of 256 by 256 (GDAL, GTiff). rasterio warns that to fill a window it reads “the entirety of one or more chunks”, and that a file without chunks is read whole, even for one pixel. Here are the lesson raster’s 12 tiles and one window:

The window starts at column 300 and row 200 and is 256 pixels square. It covers no tile completely, but it overlaps four of them, and those four are what a reader must fetch. Step through the read:

The program works the same out with numbers, for the window, for the 1 km square round the clinic, and for strips of two widths:

"""Which blocks a windowed read touches, and how many bytes it reads from strips or from 256 x 256 tiles."""
import math

from scene import HEIGHT, TRANSFORM, WIDTH


def window_for_bounds(west, south, east, north, t, width, height):
    """(col_off, row_off, cols, rows) of the pixels a box touches, rounded outwards and clipped to the raster."""
    a, _, c, _, e, f = t
    col0, col1 = max(0, math.floor((west - c) / a)), min(width, math.ceil((east - c) / a))
    row0, row1 = max(0, math.floor((north - f) / e)), min(height, math.ceil((south - f) / e))
    return col0, row0, col1 - col0, row1 - row0


def tiles(window, size=256):
    """The (block row, block col) of every tile a window overlaps."""
    col0, row0, cols, rows = window
    return [(r, c) for r in range(row0 // size, (row0 + rows - 1) // size + 1) for c in range(col0 // size, (col0 + cols - 1) // size + 1)]


def bytes_read(window, width, size=256, sample=2):
    """Bytes of one band read through strips (GDAL's default height: one strip is 8 KB or less) and through tiles."""
    col0, row0, cols, rows = window
    row_bytes = width * sample
    per_strip = max(1, 8192 // row_bytes)
    strips = (row0 + rows - 1) // per_strip - row0 // per_strip + 1
    return strips * per_strip * row_bytes, len(tiles(window, size)) * size * size * sample


x, y = 374_321.4, 2_043_456.8  # the clinic
box = window_for_bounds(x - 500, y - 500, x + 500, y + 500, TRANSFORM, WIDTH, HEIGHT)
print(f"1 km box round the clinic: window {box}, tiles {tiles(box)}")
square = (300, 200, 256, 256)
print(f"window {square}: tiles {tiles(square)}")
print(f"the raster is cut into {math.ceil(WIDTH / 256)} x {math.ceil(HEIGHT / 256)} = {math.ceil(WIDTH / 256) * math.ceil(HEIGHT / 256)} tiles of 256 x 256")

for name, width in [("lesson band, 800 pixels wide", WIDTH), ("a band 11,000 pixels wide", 11_000)]:
    strips, tiled = bytes_read(square, width)
    print(f"{name}: window {square} reads {strips:,} bytes from strips, {tiled:,} from tiles")
aligned = (256, 256, 256, 256)
print(f"aligned window {aligned}: {len(tiles(aligned))} tile, {bytes_read(aligned, 11_000)[1]:,} bytes")

1 km box round the clinic: window (382, 204, 101, 101), tiles [(0, 1), (1, 1)]
window (300, 200, 256, 256): tiles [(0, 1), (0, 2), (1, 1), (1, 2)]
the raster is cut into 4 x 3 = 12 tiles of 256 x 256
lesson band, 800 pixels wide: window (300, 200, 256, 256) reads 416,000 bytes from strips, 524,288 from tiles
a band 11,000 pixels wide: window (300, 200, 256, 256) reads 5,632,000 bytes from strips, 524,288 from tiles
aligned window (256, 256, 256, 256): 1 tile, 131,072 bytes

The 1 km square becomes a window of 101 by 101 pixels, not 100, because its edges fall inside pixels and the window rounds outwards; it overlaps two tiles. For the 256-pixel window a reader fetches four tiles, 524,288 bytes before compression, four times the window’s own 131,072. Strips can be cheaper while a raster is narrow: the lesson band is only 800 pixels wide, and its strips cost 416,000 bytes for the same window. At 11,000 pixels wide one row is already 22,000 bytes, more than 8 KB, so each strip is a single row and the window needs 256 of them, 5,632,000 bytes, while tiles still cost 524,288. A window aligned with the tiles, with offsets that are multiples of 256, reads one tile, 131,072 bytes; rasterio’s block_windows hands you exactly those windows, because reads are “most efficient when the windows match the dataset’s own block structure”.

Tiles and compression are chosen when the file is written

The block layout of a file is decided once, by whoever writes it. GDAL takes it as creation optionscreation option: A format-specific setting that GDAL passes to a driver when it writes a file or a layer, given with -co, -dsco or -lco, such as RFC7946=YES for GeoJSON or ENCODING=UTF-8 for a shapefile. More in “GDAL and OGR on the command line” of the GTiff driver: TILED=YES with BLOCKXSIZE and BLOCKYSIZE (256 by default, and tile sizes must be multiples of 16), COMPRESS (NONE by default; DEFLATE, LZW, ZSTD and others), ZLEVEL for DEFLATE (6 by default) and PREDICTOR, where 1 is none, 2 horizontal differencing and 3 floating-point prediction (GDAL, GTiff). rasterio passes the same choices as keywords when it opens a file for writing, such as tiled=True, blockxsize=256 and compress="deflate" (rasterio, Profiles).

Compression works block by block, so a tiled file can be both small and quick to read in windows. The program cuts the red band into the file’s 12 tiles and compresses each one, as plain values and after a predictorpredictor: A step before compression that rewrites each value as its difference from a neighbour, which compresses better when neighbouring pixels are alike. In GDAL, PREDICTOR=2 is horizontal differencing for integers and 3 is for floating-point data; reading undoes it exactly.:

"""The red band cut into 256 x 256 tiles and compressed with DEFLATE, without and with horizontal differencing."""
import struct
import zlib

from scene import HEIGHT, NODATA, WIDTH, bands

SIZE = 256
red, _ = bands()


def tile(block_row, block_col):
    """The values of one tile, row by row; an edge tile is filled out with nodata to its full 256 x 256."""
    out = []
    for r in range(block_row * SIZE, (block_row + 1) * SIZE):
        line = red[r][block_col * SIZE:(block_col + 1) * SIZE] if r < HEIGHT else []
        out.append(line + [NODATA] * (SIZE - len(line)))
    return out


def differenced(rows):
    """TIFF predictor 2: each value replaced by its difference from the value to its left (modulo 2**16)."""
    return [[row[0]] + [(row[i] - row[i - 1]) % 65_536 for i in range(1, SIZE)] for row in rows]


def packed(rows):
    return b"".join(struct.pack(f"<{SIZE}H", *row) for row in rows)


blocks = [(r, c) for r in range(-(-HEIGHT // SIZE)) for c in range(-(-WIDTH // SIZE))]
raw = plain = predicted = 0
for r, c in blocks:
    values = tile(r, c)
    raw += len(packed(values))
    plain += len(zlib.compress(packed(values), 6))
    predicted += len(zlib.compress(packed(differenced(values)), 6))

print(f"{len(blocks)} tiles of {SIZE} x {SIZE} uint16 values: {raw:,} bytes")
print(f"DEFLATE, level 6:               {plain:,} bytes ({raw / plain:.2f} times smaller)")
print(f"DEFLATE, level 6, predictor 2:  {predicted:,} bytes ({raw / predicted:.2f} times smaller)")
print(f"one tile, uncompressed: {SIZE * SIZE * 2:,} bytes")

12 tiles of 256 x 256 uint16 values: 1,572,864 bytes
DEFLATE, level 6:               476,745 bytes (3.30 times smaller)
DEFLATE, level 6, predictor 2:  361,117 bytes (4.36 times smaller)
one tile, uncompressed: 131,072 bytes

The 12 tiles hold 1,572,864 bytes. DEFLATE at level 6 brings them to 476,745 bytes, 3.30 times smaller, and with predictor 2 first to 361,117, 4.36 times smaller. The predictor helps because neighbouring pixels of real ground are alike: replacing each value by its difference from the one to its left turns a row such as 812, 815, 813 into 812, 3, −2 (stored modulo 65,536), small numbers that repeat, which DEFLATE stores in few bits. DEFLATE loses nothing, so every value comes back exactly. Use predictor 2 for integer bands and 3 for floating-point ones; a file written for reading in parts over a network is better still as a cloud-optimised GeoTIFF, a cloud-native formatcloud-native format: A file format that a program can use over HTTP by fetching only the byte ranges it needs, guided by metadata or an index inside the file, instead of downloading the whole file first. More in “GeoParquet and cloud-native vector data” that a later lesson of this module builds.

Practise: windows and tiles

Turn the window and block arithmetic of this lesson into code, step by step, with sample tests on the lesson raster: first the window a box on the map touches, then the tiles a window makes a reader fetch and the bytes they cost.

import math


def window_for_bounds(west, south, east, north, transform, width, height):
    """(col_off, row_off, cols, rows) of the pixels a box touches, clipped to the raster, or None outside it."""
    # Replace this line: find the first and the last column and row from the box's edges, then clip them.
    return None

from window import window_for_bounds

T = (10, 0, 370_000, 0, -10, 2_046_000)  # the lesson raster: 800 columns, 600 rows of 10 m


def test_box_round_the_clinic():
    """takes all 101 x 101 pixels that a 1 km box touches"""
    assert window_for_bounds(373_821.4, 2_042_956.8, 374_821.4, 2_043_956.8, T, 800, 600) == (382, 204, 101, 101)


def test_box_on_pixel_edges():
    """takes no extra pixel when the box edges lie on the lines between pixels"""
    assert window_for_bounds(370_000, 2_045_000, 371_000, 2_046_000, T, 800, 600) == (0, 0, 100, 100)
    assert window_for_bounds(373_000, 2_043_440, 373_010, 2_043_450, T, 800, 600) == (300, 255, 1, 1)


def test_box_partly_outside():
    """clips a box that reaches over the raster's east and south edges"""
    assert window_for_bounds(377_500, 2_039_000, 379_000, 2_040_500, T, 800, 600) == (750, 550, 50, 50)


def test_box_outside():
    """returns None for a box that misses the raster"""
    assert window_for_bounds(380_000, 2_041_000, 381_000, 2_042_000, T, 800, 600) is None
    assert window_for_bounds(371_000, 2_046_000, 372_000, 2_047_000, T, 800, 600) is None

def blocks_for_window(col_off, row_off, cols, rows, block_width, block_height):
    """[(block_row, block_col), …] of every block a window overlaps, sorted."""
    # Replace this line: find the first and the last block column and row the window reaches.
    return []


def bytes_for_window(col_off, row_off, cols, rows, block_width, block_height, sample_bytes):
    """Bytes read for one band when every block the window overlaps is read whole."""
    # Replace this line: count the blocks and multiply by the size of one block.
    return 0

from blocks import blocks_for_window, bytes_for_window


def test_window_inside_one_tile():
    """finds the one tile that holds a small window"""
    assert blocks_for_window(10, 10, 100, 100, 256, 256) == [(0, 0)]
    assert blocks_for_window(256, 256, 256, 256, 256, 256) == [(1, 1)]


def test_window_across_four_tiles():
    """finds the four tiles under the lesson's window at column 300, row 200"""
    assert blocks_for_window(300, 200, 256, 256, 256, 256) == [(0, 1), (0, 2), (1, 1), (1, 2)]


def test_window_on_the_last_column_of_a_tile():
    """does not reach into the next tile when the window ends on a tile's last column"""
    assert blocks_for_window(0, 0, 256, 1, 256, 256) == [(0, 0)]
    assert blocks_for_window(255, 0, 2, 1, 256, 256) == [(0, 0), (0, 1)]


def test_strips_and_empty_windows():
    """treats strips as blocks one raster wide, and an empty window as no blocks"""
    assert blocks_for_window(300, 200, 256, 10, 800, 5) == [(40, 0), (41, 0)]
    assert blocks_for_window(300, 200, 0, 10, 256, 256) == []


def test_bytes():
    """counts whole blocks: 4 tiles of 256 x 256 two-byte values are 524,288 bytes"""
    assert bytes_for_window(300, 200, 256, 256, 256, 256, 2) == 524_288
    assert bytes_for_window(300, 200, 256, 256, 11_000, 1, 2) == 5_632_000

Then a short self-check across the five ideas, from the transform to the predictor; it begins with the idea you have answered least well:

Interview questions

Warm-up (fresher): What does a raster’s transform tell you, and how do you find a pixel’s map coordinates?

The transform is six numbers that place the grid on the ground. In rasterio’s order, x = a·col + b·row + c and y = d·col + e·row + f: c and f are the x and y of the raster’s top-left corner, a is the pixel width, e the pixel height (negative, because rows run south) and b and d are 0 for a north-up image. Those sums give a pixel’s top-left corner; add half a pixel to the column and the row for its centre. For the lesson raster, (10, 0, 370000, 0, −10, 2046000), the pixel in row 254, column 432 has its centre at 374,325 E, 2,043,455 N. Going back, subtract the corner, divide by the pixel size and round down. GDAL stores the same numbers in the order c, a, b, f, d, e.

Two longer questions start from this lesson’s mistakes, points that miss their pixels and windows that read too much; the bank gives what each tests and a first step:

Key takeaways

• A raster’s affine transform turns a column and a row into a pixel’s top-left corner, x = a·col + b·row + c and y = d·col + e·row + f; its centre is half a pixel further in, and rasterio’s xy returns the centre by default.
• GDAL keeps the same six numbers in the order c, a, b, f, d, e, and a half-pixel slip sends three samples in four to a neighbouring pixel without any error.
• The profile names the driver, the size, the bands, the data type, the nodata value and the CRS, and a point must be in the raster’s CRS before you look it up.
• The data type fixes the range, the precision and the bytes per pixel, and the nodata value must be one the data can never take.
• A windowed read fetches every block it overlaps, whole: tiles keep that cheap on wide rasters, where strips one row high read the full width of every row.
• Tiling, compression and the predictor are chosen when a GeoTIFF is written, and predictor 2 before DEFLATE made the lesson’s tiles 4.36 times smaller instead of 3.30.