Vector and raster data models

A lake drawn as a polygon covers 69.39 ha. Stored as a grid of 100 m cells, the same lake covers 68.00 ha, and the narrow inlet on its east side has disappeared. Both numbers are honest: they come from two different ways of storing the world, and each way decides what you can measure.

A GIS stores everything either as vector features, shapes that each carry a row of facts, or as rasters, grids of cells that each hold a value. Below you take each model apart, choose between them for the thing you want to store, compute an area from a raster and see what the cell size does to it, and judge whether data made at a given map scale can decide the question in front of you.

A vector feature is a geometry with a row of attributes

A featurefeature: One thing in a vector layer, such as a ward, a road or a clinic: a geometry stored together with a row of attributes. is one thing in a layer: a ward, a road, a clinic. The QGIS introduction puts it in two parts: “A vector feature has its shape represented using geometry”, and vector features also have attributesattribute: A value stored with a feature that is not part of its shape, such as a name or a count. The attributes of a layer's features form its attribute table., “text or numerical information that describe the features” (QGIS, Vector Data). A layer is therefore a table whose rows each hold a shape.

The geometry comes in three basic types. A point is one position. A line is a list of two or more positions. A polygon is a ring of positions whose last one repeats the first, so it closes around an area, and it may hold further rings inside it that cut holes. Each type has a multi-part form, such as a multipolygon for a ward split by a river. Step through one polygon with a hole, Lake Park, from its row in the layer to its text:

The last step showed well-known textwell-known text (WKT): A plain-text form of a geometry from the OGC Simple Features standard, such as POINT (4000 5000). The same name is also used for a text format of coordinate reference systems., the plain-text form of a geometry from the OGC’s Simple Feature Access standard (OGC), which the PostGIS manual uses for its examples, such as POINT(0 0) and LINESTRING(0 0,1 1,1 2) (PostGIS manual). The program below holds four features of the town, measures each and writes it as text:

"""Four vector features: each one a geometry (a point, a line or the rings of a polygon) plus a row of attributes."""
import math

FEATURES = [
    {"type": "Point", "coordinates": (4_000, 5_000), "attributes": {"name": "C2", "kind": "clinic", "beds": 12}},
    {"type": "LineString", "coordinates": [(3_500, 6_000), (3_500, 5_000), (3_500, 2_800), (3_500, 0)],
     "attributes": {"name": "Station Road", "lanes": 2}},
    {"type": "Polygon", "coordinates": [[(2_000, 4_000), (4_000, 4_000), (4_000, 6_000), (2_000, 6_000), (2_000, 4_000)]],
     "attributes": {"name": "W02", "kind": "ward"}},
    {"type": "Polygon", "coordinates": [
        [(4_400, 4_300), (5_200, 4_300), (5_200, 4_900), (4_400, 4_900), (4_400, 4_300)],  # the park's outline
        [(4_700, 4_500), (4_700, 4_700), (4_900, 4_700), (4_900, 4_500), (4_700, 4_500)],  # a pond: a hole
    ], "attributes": {"name": "Lake Park", "kind": "park"}},
]


def signed_area(ring):
    """Shoelace formula: positive for an anticlockwise ring, negative for a clockwise one."""
    return sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(ring, ring[1:])) / 2


def length(points):
    return sum(math.dist(a, b) for a, b in zip(points, points[1:]))


def wkt(feature):
    """Well-known text, the plain-text form of a geometry in the OGC Simple Features standard."""
    pairs = lambda points: ", ".join(f"{x} {y}" for x, y in points)
    kind, coords = feature["type"], feature["coordinates"]
    if kind == "Point":
        return f"POINT ({coords[0]} {coords[1]})"
    if kind == "LineString":
        return f"LINESTRING ({pairs(coords)})"
    return "POLYGON (" + ", ".join(f"({pairs(ring)})" for ring in coords) + ")"


for f in FEATURES:
    kind, coords = f["type"], f["coordinates"]
    if kind == "LineString":
        measure = f"length {length(coords) / 1_000:.2f} km"
    elif kind == "Polygon":
        area = sum(signed_area(ring) for ring in coords)  # a hole is clockwise, so its area is taken away
        measure = f"area {area / 10_000:.0f} ha, {len(coords) - 1} hole(s)"
    else:
        measure = "no length or area"
    print(f"{f['attributes']['name']}: {kind}, {measure}, attributes {f['attributes']}")
    print(f"  {wkt(f)}")

C2: Point, no length or area, attributes {'name': 'C2', 'kind': 'clinic', 'beds': 12}
  POINT (4000 5000)
Station Road: LineString, length 6.00 km, attributes {'name': 'Station Road', 'lanes': 2}
  LINESTRING (3500 6000, 3500 5000, 3500 2800, 3500 0)
W02: Polygon, area 400 ha, 0 hole(s), attributes {'name': 'W02', 'kind': 'ward'}
  POLYGON ((2000 4000, 4000 4000, 4000 6000, 2000 6000, 2000 4000))
Lake Park: Polygon, area 44 ha, 1 hole(s), attributes {'name': 'Lake Park', 'kind': 'park'}
  POLYGON ((4400 4300, 5200 4300, 5200 4900, 4400 4900, 4400 4300), (4700 4500, 4700 4700, 4900 4700, 4900 4500, 4700 4500))

A point has no length or area; Station Road is 6.00 km long; ward W02 covers 400 ha, and Lake Park 44 ha, because its pond’s ring runs clockwise and its signed area is taken away. GeoJSON’s standard requires exactly that orientation, anticlockwise outer rings and clockwise holes (RFC 7946, section 3.1.6), but it also asks readers not to reject polygons wound the other way, so code that reads other people’s files should subtract each hole’s area instead of trusting its sign.

Features that touch, such as two wards with a common edge, each store that edge’s positions. When the two copies do not match exactly, thin gaps or overlaps called slivers appear along the border, one of the errors of digitised data that the QGIS introduction describes. Keeping shared boundaries identical is called keeping the data’s topology clean.

A raster is a grid of cells tied to the ground by six numbers

A rasterraster: A grid of equal cells, each holding one value per band, placed on the ground by a geotransform and a coordinate reference system. It suits imagery and quantities with a value everywhere. is “a grid of regularly sized pixels” whose fixed size sets its spatial resolution (QGIS, Raster Data). The QGIS introduction gives examples from 50 cm aerial photographs to 500 m satellite images. The width of one cell on the ground is the raster’s cell sizecell size (resolution): The width of one raster cell on the ground, such as 10 m or 100 m. Smaller cells show more detail and need more storage.. Each cell holds one value per bandband: 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.: a photograph has red, green and blue bands, an elevation model a single band. A cell with no measurement holds the band’s 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., a code such as -9999 that is not data.

The cells are tied to the ground by the raster’s 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., six numbers that GDAL writes as GT(0) to GT(5). A cell’s column and row turn into map coordinates with x = GT(0) + column × GT(1) + row × GT(2) and y = GT(3) + column × GT(4) + row × GT(5), counted from the top-left corner of the top-left cell (GDAL, Raster Data Model). For a north-up image GT(2) and GT(4) are zero, GT(1) is the cell width and GT(5) the cell height. Here is a small elevation raster of 100 m cells, with one cell missing its value:

Its geotransform is (5000, 100, 0, 3000, 0, −100). The height is negative because rows are written from the north edge down, so y falls as the row number grows. The program below goes from a cell to the map and back, and computes two means of the elevations:

"""A raster: one value per cell of a grid, tied to the ground by six numbers (the geotransform)."""
import math

NODATA = -9999  # the value that marks a cell with no measurement
# Elevation in metres on a grid of 5 rows and 6 columns with 100 m cells; row 0 is the north edge.
ELEVATION = [
    [512, 515, 519, 522, 524, 523],
    [514, 518, 524, 529, 530, 527],
    [515, 521, 529, 537, 535, 529],
    [513, 519, 526, NODATA, 531, 526],
    [511, 515, 520, 523, 524, 522],
]
# x of the west edge, cell width, row rotation, y of the north edge, column rotation, cell height
TRANSFORM = (5_000, 100, 0, 3_000, 0, -100)


def cell_centre(row, col, t=TRANSFORM):
    """Map coordinates of the centre of a cell: the geotransform applied to (col + 0.5, row + 0.5)."""
    x = t[0] + (col + 0.5) * t[1] + (row + 0.5) * t[2]
    y = t[3] + (col + 0.5) * t[4] + (row + 0.5) * t[5]
    return x, y


def cell_at(x, y, t=TRANSFORM):
    """The (row, col) of the cell that holds (x, y), for a grid without rotation."""
    return math.floor((y - t[3]) / t[5]), math.floor((x - t[0]) / t[1])


rows, cols = len(ELEVATION), len(ELEVATION[0])
print(f"{rows} rows x {cols} columns of {TRANSFORM[1]} m cells: x {TRANSFORM[0]} to {TRANSFORM[0] + cols * TRANSFORM[1]},"
      f" y {TRANSFORM[3] + rows * TRANSFORM[5]} to {TRANSFORM[3]}")
print(f"Cell (row 2, col 3) is centred on {cell_centre(2, 3)} and holds {ELEVATION[2][3]} m")
row, col = cell_at(5_480, 2_515)
print(f"The point (5480, 2515) is in cell (row {row}, col {col}), which holds {ELEVATION[row][col]} m")
values = [v for line in ELEVATION for v in line]
valid = [v for v in values if v != NODATA]
print(f"Mean of all 30 cells, nodata counted as a value: {sum(values) / len(values):.1f} m")
print(f"Mean of the {len(valid)} cells with data: {sum(valid) / len(valid):.1f} m")
above = sum(v > 525 for v in valid)
print(f"{above} cells are above 525 m: {above * TRANSFORM[1] ** 2 / 10_000:.0f} ha")

5 rows x 6 columns of 100 m cells: x 5000 to 5600, y 2500 to 3000
Cell (row 2, col 3) is centred on (5350.0, 2750.0) and holds 537 m
The point (5480, 2515) is in cell (row 4, col 4), which holds 524 m
Mean of all 30 cells, nodata counted as a value: 171.8 m
Mean of the 29 cells with data: 522.5 m
10 cells are above 525 m: 10 ha

Add half a cell to reach a centre: the cell in row 2, column 3 is centred on x 5000 + 3.5 × 100 = 5,350 and y 3000 − 2.5 × 100 = 2,750. Going back, divide the distance from the corner by the cell size and round down: the point x 5480, y 2515 falls in row 4, column 4.

Objects or fields: the choice decides what you can ask

Some things in the world are objects, separate and countable, with edges: wells, roads, parcels, clinics. Others are fields, quantities with a value at every place: elevation, rainfall, noise. Objects suit vectors, which keep each one’s position and identity. Fields suit rasters, which hold a value everywhere, and that is where the QGIS introduction recommends them: for information that is continuous across an area and “cannot easily be divided into vector features”.

Rasters come in two kinds, and they need different maths. A continuous raster holds measurements, such as metres of elevation, which can be averaged and compared. A categorical raster holds class codes, such as 1 for water and 3 for forest, which are names: the average of 1 and 3 is 2, the code of a class that neither cell holds. Count categorical cells, never average them.

Choosing a model also chooses what you can still ask. The program below stores the same 120 street lights three ways, as points, as counts in 50 m cells and as four totals, and asks each about one house:

"""The same 120 street lights as points, as a grid of 50 m cells and as four totals: what each model can answer."""
import math
import random

rng = random.Random(11)
# Lights along two streets that cross a 1 km square block: one runs east-west, one north-south.
lights = [(round(rng.uniform(0, 1_000), 1), round(500 + rng.uniform(-6, 6), 1)) for _ in range(60)]
lights += [(round(300 + rng.uniform(-6, 6), 1), round(rng.uniform(0, 1_000), 1)) for _ in range(60)]

CELL = 50
grid = {}  # (column, row) -> number of lights in that cell
for x, y in lights:
    grid[(int(x // CELL), int(y // CELL))] = grid.get((int(x // CELL), int(y // CELL)), 0) + 1
quarters = {}  # four totals, the way a report by area would store them
for x, y in lights:
    name = ("north" if y >= 500 else "south") + ("-east" if x >= 500 else "-west")
    quarters[name] = quarters.get(name, 0) + 1

house = (512.0, 531.0)
print(f"Points: the nearest light is {min(math.dist(house, p) for p in lights):.1f} m from the house,"
      f" and {sum(math.dist(house, p) <= 30 for p in lights)} are within 30 m")
cell = (int(house[0] // CELL), int(house[1] // CELL))
print(f"Grid: the house's 50 m cell {cell} holds {grid.get(cell, 0)} lights, but not where in the cell they stand")
print(f"Totals: north-east {quarters['north-east']}, north-west {quarters['north-west']},"
      f" south-east {quarters['south-east']}, south-west {quarters['south-west']}; nothing about any one street")
print(f"Stored: {len(lights)} points, {len(grid)} cells with a light out of {(1_000 // CELL) ** 2}, 4 totals")

Points: the nearest light is 30.1 m from the house, and 0 are within 30 m
Grid: the house's 50 m cell (10, 10) holds 2 lights, but not where in the cell they stand
Totals: north-east 15, north-west 50, south-east 12, south-west 43; nothing about any one street
Stored: 120 points, 62 cells with a light out of 400, 4 totals

Only the points can say that the nearest light is 30.1 m away, so none is within 30 m. The grid knows that the house’s cell holds 2 lights, but not where in the cell they stand. The totals know only that the north-east quarter has 15. Each step of aggregation is cheaper to store and answers fewer questions. Decide for your own data with the table below:

Most answers sort themselves out: separate objects with exact edges become points, lines or polygons, measurements with a value everywhere become a continuous raster, and classes with a value everywhere a categorical one.

Later modules of this track meet models beyond these two. Point clouds, such as those from laser scanning, hold millions of measured points, each with a height; PDAL is “a C++ library for translating and manipulating point cloud data” (PDAL). Networks store roads as junctions and links, as in the first lesson’s drive. Discrete global grids cover the Earth with cells of fixed shapes; H3, for example, indexes the world in hexagonal cells (H3).

Coarser cells change the area you measure

Turning a polygon into a raster, rasterising it, needs a rule for the cells on its edge. The common one keeps a cell when its centre lies inside the polygon; rasterio’s rasterize burns in “only pixels whose center is within the polygon” unless its all_touched option is set (rasterio). Here is one lake, a polygon with a 40 m inlet, rasterised that way with 30 m and with 100 m cells:

The 30 m cells follow the shore in small steps and keep part of the inlet; the 100 m cells keep none of it. The program behind the figure measures the polygon and its rasters with 10, 30 and 100 m cells, under the cell-centre rule and under the rule that counts every cell the lake touches:

"""One lake, stored as a polygon and as rasters of 10, 30 and 100 m cells: what each measures."""
import math

# The lake's outline in metres, anticlockwise, with a narrow inlet 40 m wide reaching east from (1200, 500).
LAKE = [
    (200, 300), (500, 120), (900, 100), (1200, 260), (1200, 480), (1600, 480), (1600, 520), (1200, 520),
    (1180, 720), (950, 880), (600, 900), (300, 780), (150, 560), (200, 300),
]


def inside(ring, x, y):
    """Ray casting: True if (x, y) is inside the ring."""
    crossings = 0
    for (x1, y1), (x2, y2) in zip(ring, ring[1:]):
        if (y1 > y) != (y2 > y) and x < x1 + (y - y1) * (x2 - x1) / (y2 - y1):
            crossings += 1
    return crossings % 2 == 1


def shoelace(ring):
    """The area enclosed by a ring."""
    return abs(sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(ring, ring[1:]))) / 2


def outline_length(ring):
    return sum(math.dist(a, b) for a, b in zip(ring, ring[1:]))


def crosses(a, b, xmin, ymin, xmax, ymax):
    """True if the segment a-b meets the box (Liang-Barsky clipping)."""
    (x1, y1), (x2, y2) = a, b
    dx, dy = x2 - x1, y2 - y1
    t0, t1 = 0.0, 1.0
    for p, q in ((-dx, x1 - xmin), (dx, xmax - x1), (-dy, y1 - ymin), (dy, ymax - y1)):
        if p == 0 and q < 0:
            return False
        if p < 0:
            t0 = max(t0, q / p)
        elif p > 0:
            t1 = min(t1, q / p)
    return t0 <= t1


def rasterise(ring, cell, all_touched=False):
    """Cells of a grid aligned to (0, 0): those whose centre is inside, or with all_touched every cell the lake meets."""
    xs, ys = [x for x, _ in ring], [y for _, y in ring]
    cells = set()
    for c in range(int(min(xs) // cell), int(max(xs) // cell) + 1):
        for r in range(int(min(ys) // cell), int(max(ys) // cell) + 1):
            box = (c * cell, r * cell, (c + 1) * cell, (r + 1) * cell)
            if inside(ring, (c + 0.5) * cell, (r + 0.5) * cell) or (
                all_touched and any(crosses(a, b, *box) for a, b in zip(ring, ring[1:]))
            ):
                cells.add((c, r))
    return cells


def raster_perimeter(cells, cell):
    """Edges between a water cell and a land cell, times the cell size: the staircase outline."""
    edges = sum((c + dc, r + dr) not in cells for c, r in cells for dc, dr in ((1, 0), (-1, 0), (0, 1), (0, -1)))
    return edges * cell


area, outline = shoelace(LAKE), outline_length(LAKE)
staircase = sum(abs(x2 - x1) + abs(y2 - y1) for (x1, y1), (x2, y2) in zip(LAKE, LAKE[1:]))
print(f"Polygon:      area {area / 10_000:.2f} ha, outline {outline / 1_000:.2f} km, inlet 1.60 ha")
print(f"Its outline as a staircase, |dx| + |dy| along every side: {staircase / 1_000:.2f} km")
for cell in (10, 30, 100):
    cells = rasterise(LAKE, cell)
    a, p = len(cells) * cell**2, raster_perimeter(cells, cell)
    inlet = sum(cell**2 for c, r in cells if (c + 0.5) * cell > 1_200)  # water cells east of x = 1200
    print(f"{cell:>3} m cells: area {a / 10_000:.2f} ha ({100 * (a - area) / area:+.2f} %),"
          f" outline {p / 1_000:.2f} km ({100 * (p - outline) / outline:+.0f} %), inlet {inlet / 10_000:.2f} ha")
    touched = len(rasterise(LAKE, cell, all_touched=True)) * cell**2
    print(f"             every cell the lake touches: {touched / 10_000:.2f} ha ({100 * (touched - area) / area:+.2f} %)")

Polygon:      area 69.39 ha, outline 3.84 km, inlet 1.60 ha
Its outline as a staircase, |dx| + |dy| along every side: 4.50 km
 10 m cells: area 69.42 ha (+0.04 %), outline 4.50 km (+17 %), inlet 1.60 ha
             every cell the lake touches: 72.28 ha (+4.16 %)
 30 m cells: area 69.03 ha (-0.52 %), outline 4.50 km (+17 %), inlet 1.17 ha
             every cell the lake touches: 77.67 ha (+11.93 %)
100 m cells: area 68.00 ha (-2.00 %), outline 3.60 km (-6 %), inlet 0.00 ha
             every cell the lake touches: 95.00 ha (+36.91 %)

Four things show in the output, each measured rather than assumed:

• The area drifts, but slowly. It moves from 69.42 ha at 10 m to 69.03 ha at 30 m and 68.00 ha at 100 m, against the polygon’s 69.39 ha. Cells on the shore are all water or all land, and their errors partly cancel.
• Narrow features vanish. The 40 m inlet keeps all of its 1.60 ha at 10 m, 1.17 ha at 30 m and nothing at 100 m, because the strip, from y 480 to 520, falls between two rows of 100 m cell centres, at 450 and 550.
• The rule matters as much as the size. Counting every cell the lake touches, the idea behind rasterio’s all_touched option, gives 72.28 ha at 10 m and 95.00 ha at 100 m, 37 % too much, because every shore cell now counts as water.
• The outline never comes right. A raster’s outline is a staircase of cell edges, and a staircase measures each side as its east-west plus its north-south extent. For this lake that is 4.50 km, at 10 m and at 30 m alike, 17 % more than the polygon’s 3.84 km. The 100 m outline comes out shorter, 3.60 km, only because the inlet and the western tip of the lake have dropped out of its cells. Smaller cells fix the area, not the perimeter.

Scale and accuracy travel with the data

Data made from a map keeps that map’s limits. A map scalemap scale: The ratio between a length on a map and the same length on the ground. At 1:50,000, 1 mm on the map stands for 50 m. is a ratio of lengths: at 1:50,000, 1 mm on the map is 50 m on the ground. Features too small for that scale were simplified or left out when the map was drawn, and the QGIS introduction warns that vector data digitised from a map inherits the scale problems of that map. The program below turns scales and a stated error into distances on the ground:

"""What a map scale and a stated positional accuracy mean on the ground."""
import math
import random

for scale in (10_000, 50_000, 250_000):
    print(f"1:{scale:,}: 1 mm on the map is {scale / 1_000:,.0f} m on the ground, a 0.25 mm line is {0.25 * scale / 1_000:,.1f} m wide")

RMSE = 12.5  # the boundary layer's stated horizontal error in metres, as a root mean square distance
sigma = RMSE / math.sqrt(2)  # the same error split equally between the east and north directions
rng = random.Random(5)
errors = sorted(math.hypot(rng.gauss(0, sigma), rng.gauss(0, sigma)) for _ in range(100_000))
print(f"With an RMSE of {RMSE} m, 95 % of 100,000 simulated errors are within {errors[94_999]:.1f} m"
      f" (by formula {RMSE * math.sqrt(-math.log(0.05)):.1f} m)")
for d in (5, 10, 20, 30):
    chance = 0.5 * math.erfc(d / (sigma * math.sqrt(2)))  # the error across the line must exceed d
    print(f"A house {d} m inside the drawn boundary is truly on the other side with a chance of {100 * chance:.2f} %")

1:10,000: 1 mm on the map is 10 m on the ground, a 0.25 mm line is 2.5 m wide
1:50,000: 1 mm on the map is 50 m on the ground, a 0.25 mm line is 12.5 m wide
1:250,000: 1 mm on the map is 250 m on the ground, a 0.25 mm line is 62.5 m wide
With an RMSE of 12.5 m, 95 % of 100,000 simulated errors are within 21.7 m (by formula 21.6 m)
A house 5 m inside the drawn boundary is truly on the other side with a chance of 28.58 %
A house 10 m inside the drawn boundary is truly on the other side with a chance of 12.89 %
A house 20 m inside the drawn boundary is truly on the other side with a chance of 1.18 %
A house 30 m inside the drawn boundary is truly on the other side with a chance of 0.03 %

At 1:50,000 a pen line 0.25 mm wide is itself 12.5 m wide on the ground. Accuracy is stated separately, often as an RMSERMSE (root mean square error): The square root of the mean of the squared errors. For positions, the typical size of the distance between where things are recorded and where they really are., the root mean square of the distances between recorded and true positions. The program takes a boundary layer that states 12.5 m RMSE and assumes its errors are random, alike in every direction and independent from point to point. Then 95 % of position errors fall within 21.6 m, which it computes as the RMSE times √(−ln 0.05) and confirms with 100,000 simulated errors (21.7 m). That factor, 1.7308, is the one the US National Standard for Spatial Data Accuracy uses to report horizontal accuracy at 95 % confidence (FGDC).

The last four lines answer the question such data is often asked to decide. A house 30 m inside the drawn boundary is safely inside, with a 0.03 % chance of being on the other side. A house 10 m inside has a 12.89 % chance, and one 5 m inside 28.58 %. A boundary like this is right for district maps and statistics, and wrong for deciding which side of it a particular house stands. That needs data captured at a larger scale, such as a parcel survey.

Practise: compute with rasters

Turn this lesson’s raster arithmetic into code in two steps, both tested on the page. One measures the area of each class in a land-cover raster and leaves the nodata cells out; the other goes from a cell to the map and back with a geotransform, including points on the line between two cells and points just outside the raster:

def area_by_class(grid, cell_size, nodata):
    """{class code: area in hectares} of a categorical raster; cells holding nodata are left out."""
    # Replace this line: count the cells of each class, then turn the counts into hectares.
    return {}

from classes import area_by_class

GRID = [
    [1, 1, 2, 2],
    [1, 3, 2, 2],
    [0, 3, 3, 2],
]


def test_each_class_gets_its_area():
    """counts the cells of each class and turns them into hectares (100 m cells are 1 ha each)"""
    assert area_by_class(GRID, 100, 0) == {1: 3.0, 2: 5.0, 3: 3.0}


def test_nodata_is_left_out():
    """leaves out the nodata code, whatever number it is"""
    got = area_by_class([[255, 1], [1, 255]], 30, 255)
    assert set(got) == {1}
    assert abs(got[1] - 0.18) < 1e-9


def test_only_nodata():
    """returns an empty dictionary when every cell is nodata"""
    assert area_by_class([[-9999, -9999]], 10, -9999) == {}


def test_cell_size_matters():
    """gives 0.01 ha per 10 m cell, 0.09 ha per 30 m cell and 1 ha per 100 m cell"""
    water = [[1] * 10 for _ in range(10)]  # 100 cells of water
    for size, hectares in [(10, 1.0), (30, 9.0), (100, 100.0)]:
        assert abs(area_by_class(water, size, 0)[1] - hectares) < 1e-9

import math


def cell_centre(row, col, transform):
    """(x, y) of the centre of a cell of a north-up raster with the geotransform (x0, width, 0, y0, 0, height)."""
    # Replace this line: start at the corner (x0, y0) and move half a cell further than col and row.
    return 0, 0


def cell_at(x, y, transform, rows, cols):
    """(row, col) of the cell that holds (x, y), or None when the point is outside the raster."""
    # Replace this line: divide the distance from the corner by the cell size and round down.
    return None

from cells import cell_at, cell_centre

DEM = (5_000, 100, 0, 3_000, 0, -100)  # the lesson's elevation raster: 5 rows, 6 columns


def test_centres_of_the_lesson_raster():
    """puts cell (2, 3) at (5350, 2750) and cell (0, 0) half a cell from the corner"""
    assert cell_centre(2, 3, DEM) == (5_350, 2_750)
    assert cell_centre(0, 0, DEM) == (5_050, 2_950)
    assert cell_centre(4, 5, DEM) == (5_550, 2_550)


def test_points_inside():
    """finds the cell of a point, rows counting southwards"""
    assert cell_at(5_480, 2_515, DEM, 5, 6) == (4, 4)
    assert cell_at(5_001, 2_999, DEM, 5, 6) == (0, 0)


def test_points_on_cell_edges():
    """puts a point on a line between cells into the cell east or south of it"""
    assert cell_at(5_100, 2_950, DEM, 5, 6) == (0, 1)
    assert cell_at(5_050, 2_900, DEM, 5, 6) == (1, 0)
    assert cell_at(5_000, 3_000, DEM, 5, 6) == (0, 0)


def test_points_outside():
    """returns None beyond every edge, just west of the raster included"""
    for x, y in [(4_990, 2_950), (5_600, 2_950), (5_050, 3_001), (5_050, 2_500), (9_000, 0)]:
        assert cell_at(x, y, DEM, 5, 6) is None, (x, y)


def test_another_raster():
    """works for any cell size: a raster of 30 m cells with its corner at (0, 600)"""
    small = (0, 30, 0, 600, 0, -30)
    assert cell_centre(0, 0, small) == (15, 585)
    assert cell_at(61, 539, small, 20, 20) == (2, 2)

A short self-check follows, over the five ideas from features to scale, and it opens with the idea you have answered least well:

Interview questions

Warm-up (fresher): When would you store data as a raster rather than as vectors?

Store it as a raster when the thing has a value at every place, such as elevation, temperature or rainfall, or when it arrives as a grid already, as satellite and aerial images do. Rasters make cell-by-cell maths easy: adding layers, comparing them, averaging within an area. Store it as vectors when the things are separate objects with exact edges and identities, such as parcels, roads or clinics, or when lengths and boundaries must be exact: a lake polygon measured 69.39 ha, while 100 m cells gave 68.00 ha and lost a 40 m inlet. Many projects use both, for example a land-cover raster summarised inside ward polygons.

The lake and the 1:50,000 boundary come back in two longer questions; for each, the bank names what it tests and gives a first step:

Key takeaways

• A vector feature is a geometry, a point, a line or a polygon with optional holes, stored with a row of attributes, and well-known text writes the geometry as one line of text.
• A raster is a grid of equal cells with one value per band, tied to the ground by its geotransform, and cells holding the nodata value must be left out of every statistic.
• Objects suit vectors and fields suit rasters; categorical rasters are counted, never averaged, and every step of aggregation answers fewer questions.
• Coarser cells make a raster’s area drift and make narrow features vanish, and a raster’s staircase outline never shrinks to the true shoreline, however small the cells get.
• Data keeps the scale and the accuracy it was made with: a 1:50,000 boundary with 12.5 m RMSE can put a house 10 m inside it on the wrong side about one time in eight.