GIS & Geospatial Data (QGIS, GeoPandas, PostGIS) Module 6 – Raster data and remote sensing
Resampling and reprojecting rasters
Choose nearest, bilinear, cubic, average or mode by what a raster holds, predict the grid a reprojection makes, and put rasters on one grid.
What you will learn
- Choose nearest, bilinear, cubic, average or mode resampling for the values a raster holds
- Predict the grid of a reprojected raster, and set it explicitly with -tr, -tap and -te
- Match rasters to one grid of CRS, resolution, origin and extent before combining them
- Explain why pixels change size and shape on the ground when a raster is reprojected to degrees
Before you start
What you’ll learn · Before you start
What you will learn
- Choose nearest, bilinear, cubic, average or mode resampling for the values a raster holds
- Predict the grid of a reprojected raster, and set it explicitly with -tr, -tap and -te
- Match rasters to one grid of CRS, resolution, origin and extent before combining them
- Explain why pixels change size and shape on the ground when a raster is reprojected to degrees
Before you start
On this page
Move a land-cover raster of two classes half a cell with bilinear resampling, and 198 of its new cells become class 2, a town that is not there. Every resampling decides values the old grid never held, and reprojecting rasters always resamples them, because the pixels of a new CRS never land on the old ones.
Here you choose a method by what the values mean, see the grid a reprojected raster gets and why its pixels are square in degrees but not in metres, and follow one pixel through a warp (raster reprojection): Putting a raster onto a new grid, usually in another CRS: the new grid is chosen first, then each of its pixels is filled by finding its centre on the old raster and resampling there.. Then you measure each method’s error, put rasters on one grid before combining them, and compute areas on a grid of degrees.
Every resampling method decides values the old grid never held
resampling: Computing the cells of a new grid from an old one, whenever the cell size, the grid's position or the CRS changes. Each method, such as nearest, bilinear, cubic, average or mode, decides the new values in its own way. computes the cells of a new grid from an old one whenever the cells change: a
coarser or finer cell size, a grid moved or turned, a new CRS. rasterio’s guide puts it plainly: “all methods involve
some interpolation” (rasterio, Resampling).
gdalwarp names the common methods (GDAL, gdalwarp):
- Nearest takes the value of the old cell under the new cell’s centre: gdalwarp’s default for most formats, the fastest, and of the “worst interpolation quality”.
- Bilinear weighs the four old cell centres round the new centre by their distance.
- Cubic fits a smooth curve through the sixteen old centres round it.
- Average takes the mean of all the old cells inside the new cell, leaving nodata out.
- Mode takes the value that occurs most often among them.
Which one fits depends on what the numbers mean. A categorical raster holds codes, 1 for water and 3 for soil; a continuous one holds measurements, such as metres or NDVI. Here is a grid moved half a cell over two classes:
A grid moved half a cell over two classes, resampled with nearest and with bilinear
Text description of the diagram
At the top, a grid of 4 by 4 old cells: the two left columns hold code 1, water, shaded light blue, and the two right columns hold code 3, soil, shaded grey. Over them, a dashed blue grid of 3 by 3 new cells is moved half a cell to the right and half a cell down, so each new cell's centre, marked with a blue dot, sits on the corner where four old cells meet.
Below, two small grids of 3 by 3 show the new cells. With nearest, each new cell takes the code of one old cell: the columns read 1, 3 and 3, and only codes that exist appear. With bilinear, each new cell averages its four old cells: the left column reads 1, the right column 3, and the middle column, whose centres sit between water and soil, reads 2, shaded and with a dashed outline. The caption says that bilinear averages 1 and 3 into 2, a class neither input holds, while nearest keeps 1 and 3.
The program does the same to the field of the lesson raster that the river crosses, rows 240 to 319 and columns 100 to 199, where only water and soil occur:
"""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"""A land-cover raster of two classes, resampled four ways: what each method does to the classes."""
from collections import Counter
from scene import cover
CODES = {"water": 1, "town": 2, "soil": 3, "crop": 4}
# Rows 240 to 319 and columns 100 to 199 of the lesson raster: a bare field with the river through it.
grid = [[CODES[cover(r, c)] for c in range(100, 200)] for r in range(240, 320)]
rows, cols = len(grid), len(grid[0])
print("source, 10 m:", dict(sorted(Counter(v for row in grid for v in row).items())))
def nearest_shifted(g):
"""Each new cell, 5 m east and 5 m south of an old one, takes the value of the old cell its centre falls in."""
return [[g[r + 1][c + 1] for c in range(cols - 1)] for r in range(rows - 1)]
def bilinear_shifted(g):
"""Each new cell's centre lies midway between four old centres: the mean of the four, rounded to a whole code."""
return [[int((g[r][c] + g[r][c + 1] + g[r + 1][c] + g[r + 1][c + 1]) / 4 + 0.5) for c in range(cols - 1)] for r in range(rows - 1)]
def blocks(g, size=3):
"""The 3 x 3 blocks of old cells that make each new 30 m cell."""
return [[[g[r + i][c + j] for i in range(size) for j in range(size)] for c in range(0, cols - size + 1, size)] for r in range(0, rows - size + 1, size)]
def count(g):
return dict(sorted(Counter(v for row in g for v in row).items()))
print("shifted 5 m, nearest: ", count(nearest_shifted(grid)))
print("shifted 5 m, bilinear:", count(bilinear_shifted(grid)))
coarse = blocks(grid)
print("30 m cells, average: ", count([[int(sum(b) / len(b) + 0.5) for b in row] for row in coarse]))
print("30 m cells, mode: ", count([[Counter(b).most_common(1)[0][0] for b in row] for row in coarse]))
print("code 2 is 'town': every 2 above is a town that is not there")Output
source, 10 m: {1: 3500, 3: 4500}
shifted 5 m, nearest: {1: 3465, 3: 4356}
shifted 5 m, bilinear: {1: 3359, 2: 198, 3: 4264}
30 m cells, average: {1: 364, 2: 42, 3: 452}
30 m cells, mode: {1: 385, 3: 473}
code 2 is 'town': every 2 above is a town that is not there
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 classes.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
Bilinear made 198 cells of class 2, the lesson’s code for a town, and average at 30 m made 42; nearest and mode kept only 1 and 3. Codes must never be interpolated: for classes, use nearest when the cells keep their size and mode when they grow. For measurements, rasterio’s guide warns that nearest “may not be suitable for continuous data” and points to bilinear and cubic; when the cells grow, average gives each new cell the mean of the old cells inside it.
Checkpoint
Check: Resampling methods and what they keep
Reprojecting rasters builds a new grid, then resamples it
Reprojecting rasters starts with a grid, not with values. GDAL documents the rule for the grid it suggests: walk the old raster’s edges, about 20 points along each, transform them, and take the box round them as the new extent; then pick one square pixel size, so that the diagonal from the top-left corner to the bottom-right one holds as many pixels as before (GDAL, GDALSuggestedWarpOutput). The program applies that rule to the lesson raster, 800 by 600 pixels of 10 m in UTM zone 43N, for EPSG:4326, with a module of its own for the UTM projection:
"""UTM zone 43N on WGS 84, both ways, and the size of a degree on the ground: what the lesson's programs reproject with.
The transverse Mercator projection is computed with Krüger's series to the sixth order of the ellipsoid's third
flattening, as Karney (2011) gives it, which is good to well under a millimetre across the zone.
"""
import math
A, F = 6_378_137.0, 1 / 298.257223563 # WGS 84 semi-major axis (m) and flattening
K0, LON0, FALSE_EASTING = 0.9996, math.radians(75.0), 500_000.0 # zone 43N: central meridian 75° E
E = math.sqrt(F * (2 - F)) # eccentricity
N = F / (2 - F) # third flattening
RECT = A / (1 + N) * (1 + N**2 / 4 + N**4 / 64 + N**6 / 256) # the rectifying radius
ALPHA = [
N / 2 - 2 * N**2 / 3 + 5 * N**3 / 16 + 41 * N**4 / 180 - 127 * N**5 / 288 + 7891 * N**6 / 37800,
13 * N**2 / 48 - 3 * N**3 / 5 + 557 * N**4 / 1440 + 281 * N**5 / 630 - 1983433 * N**6 / 1935360,
61 * N**3 / 240 - 103 * N**4 / 140 + 15061 * N**5 / 26880 + 167603 * N**6 / 181440,
49561 * N**4 / 161280 - 179 * N**5 / 168 + 6601661 * N**6 / 7257600,
34729 * N**5 / 80640 - 3418889 * N**6 / 1995840,
212378941 * N**6 / 319334400,
]
BETA = [
N / 2 - 2 * N**2 / 3 + 37 * N**3 / 96 - N**4 / 360 - 81 * N**5 / 512 + 96199 * N**6 / 604800,
N**2 / 48 + N**3 / 15 - 437 * N**4 / 1440 + 46 * N**5 / 105 - 1118711 * N**6 / 3870720,
17 * N**3 / 480 - 37 * N**4 / 840 - 209 * N**5 / 4480 + 5569 * N**6 / 90720,
4397 * N**4 / 161280 - 11 * N**5 / 504 - 830251 * N**6 / 7257600,
4583 * N**5 / 161280 - 108847 * N**6 / 3991680,
20648693 * N**6 / 638668800,
]
def to_utm(lon, lat):
"""Easting and northing (m) in UTM zone 43N of a longitude and latitude in degrees."""
phi, lam = math.radians(lat), math.radians(lon) - LON0
t = math.sinh(math.atanh(math.sin(phi)) - E * math.atanh(E * math.sin(phi)))
xi0, eta0 = math.atan2(t, math.cos(lam)), math.atanh(math.sin(lam) / math.sqrt(1 + t * t))
xi = xi0 + sum(a * math.sin(2 * j * xi0) * math.cosh(2 * j * eta0) for j, a in enumerate(ALPHA, 1))
eta = eta0 + sum(a * math.cos(2 * j * xi0) * math.sinh(2 * j * eta0) for j, a in enumerate(ALPHA, 1))
return FALSE_EASTING + K0 * RECT * eta, K0 * RECT * xi
def from_utm(x, y):
"""Longitude and latitude in degrees of an easting and northing in UTM zone 43N."""
xi, eta = y / (K0 * RECT), (x - FALSE_EASTING) / (K0 * RECT)
xi0 = xi - sum(b * math.sin(2 * j * xi) * math.cosh(2 * j * eta) for j, b in enumerate(BETA, 1))
eta0 = eta - sum(b * math.cos(2 * j * xi) * math.sinh(2 * j * eta) for j, b in enumerate(BETA, 1))
tau0 = math.sin(xi0) / math.sqrt(math.sinh(eta0) ** 2 + math.cos(xi0) ** 2)
tau = tau0
for _ in range(5): # Newton's method for the conformal latitude
s = math.sinh(E * math.atanh(E * tau / math.sqrt(1 + tau * tau)))
guess = tau * math.sqrt(1 + s * s) - s * math.sqrt(1 + tau * tau)
tau += (tau0 - guess) / math.sqrt(1 + guess * guess) * (1 + (1 - E * E) * tau * tau) / ((1 - E * E) * math.sqrt(1 + tau * tau))
return math.degrees(LON0 + math.atan2(math.sinh(eta0), math.cos(xi0))), math.degrees(math.atan(tau))
def metres_per_degree(lat):
"""(east-west, north-south) length on the ground of one degree at a latitude, on the ellipsoid."""
phi = math.radians(lat)
w = 1 - E * E * math.sin(phi) ** 2
across = A * math.cos(phi) / math.sqrt(w) # radius of the parallel
along = A * (1 - E * E) / w**1.5 # radius of curvature of the meridian
return math.radians(1) * across, math.radians(1) * along"""The grid the lesson raster gets in EPSG:4326, by the rule GDAL documents for its suggested warp output."""
import math
from scene import HEIGHT, TRANSFORM, WIDTH
from utm import from_utm, metres_per_degree
a, _, c, _, e, f = TRANSFORM
west, east, north, south = c, c + a * WIDTH, f, f + e * HEIGHT
# 1. Walk round the edges, 20 steps along each, and transform the points: the box round them is the output extent.
steps = [i / 20 for i in range(20)]
ring = [(west + (east - west) * s, north) for s in steps] + [(east, north + (south - north) * s) for s in steps]
ring += [(east + (west - east) * s, south) for s in steps] + [(west, south + (north - south) * s) for s in steps]
lons, lats = zip(*(from_utm(x, y) for x, y in ring))
print(f"extent: {min(lons):.6f} to {max(lons):.6f} E, {min(lats):.6f} to {max(lats):.6f} N")
# 2. One square pixel size, so that the top-left to bottom-right diagonal keeps its 1,000 pixels.
(lon1, lat1), (lon2, lat2) = from_utm(west, north), from_utm(east, south)
size = math.hypot(lon2 - lon1, lat2 - lat1) / math.hypot(WIDTH, HEIGHT)
cols, rows = round((max(lons) - min(lons)) / size), round((max(lats) - min(lats)) / size)
print(f"pixel size: {size:.4e} degrees; grid of {cols} columns and {rows} rows (was {WIDTH} by {HEIGHT})")
# 3. The same square of degrees, measured on the ground at the raster's middle latitude.
across, along = metres_per_degree((max(lats) + min(lats)) / 2)
print(f"one pixel on the ground: {size * across:.2f} m east-west by {size * along:.2f} m north-south")
# 4. The source is a slightly turned shape inside that box: the box's corners hold no data.
footprint = abs(sum(lons[i] * lats[i - 1] - lons[i - 1] * lats[i] for i in range(len(ring)))) / 2
box = (max(lons) - min(lons)) * (max(lats) - min(lats))
print(f"the source covers {footprint / box:.2%} of the box; about {round((1 - footprint / box) * cols * rows):,} of the {cols * rows:,} new pixels get nodata")Output
extent: 73.768549 to 73.844682 E, 18.446250 to 18.500944 N pixel size: 9.3189e-05 degrees; grid of 817 columns and 587 rows (was 800 by 600) one pixel on the ground: 9.84 m east-west by 10.31 m north-south the source covers 98.64% of the box; about 6,505 of the 479,579 new pixels get nodata
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 grid.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
The new grid has 817 by 587 pixels of 9.3189e-05° each, square in degrees. On the ground at 18.47° N that square is 9.84 m east to west and 10.31 m north to south, because a degree of longitude shrinks with the latitude’s cosine while a degree of latitude barely changes. The source also sits slightly turned inside the new box and covers 98.64 % of it, so about 6,505 new pixels in its corners get nodata.
Each new pixel is then filled by working backwards, from the new grid to the old one. Follow one:
The drawing shows the whole scene; its 7 steps are listed below.
Every step as text
- Start from the new grid. A warp works backwards. For every pixel of the new grid, 817 by 587 pixels in degrees, it asks where on the old raster that pixel's ground lies.
- The centre in degrees. The new grid's corner and pixel size, 9.3189e-05 degrees, put the centre of row 300, column 400 at 73.805871 E, 18.472941 N.
- Back to the old CRS. The inverse transformation takes that point into UTM zone 43N, the old raster's CRS, at 373,920.3 E and 2,042,927.5 N.
- A place in the old grid. The old transform turns the point into column 392.03 and row 307.25, a point inside the old pixel in row 307, column 392.
- Resample. Nearest takes the value of that one pixel, water with red 333. Bilinear would weigh the four old pixels round the point, and cubic sixteen.
- Write the value. The value goes into row 300, column 400 of the new grid. Every new pixel is filled this way, so the new grid has no gaps, which mapping the old pixels forwards would leave.
- A corner of the new grid. Row 0, column 0 comes back to row −4.77 of the old grid, outside it. That pixel gets nodata, like about 6,505 of the 479,579 new pixels, in the corners of the box.
Key idea
A reprojected raster gets a new grid first, square in the new CRS’s units, and every pixel of it is resampled from the old raster by working backwards.
Checkpoint
Check: The grid of a reprojected raster
Each method's error, measured, and how to choose Pro: opens with a Pro, Premium or Ultimate pass
The rule of thumb above has numbers behind it. Here every method is measured on a hill whose true heights are known everywhere and on the lesson’s NDVI, and the results become a table that answers your own case.
In this part:
- Finer cells: nearest, bilinear and cubic against a known hill
- Coarser cells: average keeps the mean, nearest keeps one pixel
- Choosing with the table, and the defaults to watch
1 runnable example (Python) · about 8 minutes
Opening this part with your pass…
This part opens with any MySmartCoPilot pass.
- Pro pass · 1 week₹79$2.99
- 30 days₹179$5.99
- 1 year₹1,699$45
One payment; a pass never renews. The price shown is the total you pay: no fees or charges are added at checkout.
Passes can’t be bought right now, so this part can’t be opened at the moment. The free part above stands on its own. If you have a pass, enter your licence key.
It also opens every other in-depth part in Learn and every Pro tool. Sales are final, except where the law gives you a right to withdraw, for example in the EU and the UK (refund policy).
Putting rasters on one grid before combining them Pro: opens with a Pro, Premium or Ultimate pass
Cell-by-cell algebra pairs cells by row and column, so two rasters must share one grid before they are combined. Here you measure what pairing two grids 15 m apart costs, and build one grid with target-aligned pixels.
In this part:
- What one grid means
- Two grids 15 m apart, measured
- -tr, -tap and -te make the grid explicit
- Put two rasters on one grid with target-aligned pixels
- Results of the sample tests
1 runnable example (Python) · 1 exercise with sample tests · about 5 minutes
Opening this part with your pass…Pro part. It opens with the same pass as the part above.
Areas on a grid of degrees, with numbers Pro: opens with a Pro, Premium or Ultimate pass
A pixel of degrees covers less ground the further it lies from the equator. Here the lesson’s river is measured after reprojection three ways, and one pixel’s area is followed across India.
In this part:
- The river measured three ways
- One pixel of 0.0001 degrees across India
- When to measure on a grid of degrees at all
1 runnable example (Python) · about 3 minutes
Opening this part with your pass…Pro part. It opens with the same pass as the part above.
Practise: the grid of a reprojection
The exercise turns GDAL’s rule for the suggested grid into code and checks it on the lesson raster; then the self-check covers the first two sections, the methods and the new grid:
Exercise · Easy · Python
Work out the grid a reprojected raster gets
Write suggested_grid(extent, top_left, bottom_right, width, height) in suggest.py, following the rule GDAL documents for the output it suggests when a raster is reprojected:
extentis(west, south, east, north), the box round the raster's edges after they were transformed;top_leftandbottom_rightare where the raster's top-left and bottom-right corners land, as(x, y);- the pixel size is square, chosen so that the diagonal from
top_lefttobottom_rightholds as many pixels as the diagonal of the old raster ofwidthbyheightpixels; - the new grid has
round((east - west) / size)columns andround((north - south) / size)rows.
Return (size, columns, rows). The sample tests use the lesson raster, 800 by 600 pixels in UTM zone 43N, reprojected to degrees, and a raster whose corners do not move at all.
Starter code · suggest.py
import math
def suggested_grid(extent, top_left, bottom_right, width, height):
"""(pixel size, columns, rows) of the grid GDAL's rule suggests for a reprojected raster."""
# Replace this line: one square pixel size from the two diagonals, then the size of the grid.
return 0.0, 0, 0The sample tests · test_suggest.py
from suggest import suggested_grid
EXTENT = (73.76854884148034, 18.446249951458224, 73.84468186084324, 18.50094444522425)
TOP_LEFT = (73.76854884148034, 18.500466524524956)
BOTTOM_RIGHT = (73.84468186084324, 18.446726373384184)
def test_lesson_raster_in_degrees():
"""gives the lesson's grid in degrees: 817 by 587 pixels of 9.3189e-05 degrees"""
size, cols, rows = suggested_grid(EXTENT, TOP_LEFT, BOTTOM_RIGHT, 800, 600)
assert round(size, 9) == 9.3189e-05
assert (cols, rows) == (817, 587)
def test_a_raster_that_does_not_move():
"""keeps the grid of a raster whose corners stay where they were"""
size, cols, rows = suggested_grid((0, 0, 8000, 6000), (0, 6000), (8000, 0), 800, 600)
assert (size, cols, rows) == (10.0, 800, 600)
def test_rounding_of_the_grid():
"""rounds the columns and rows to the nearest whole number"""
size, cols, rows = suggested_grid((0, 0, 1004, 1006), (0, 1000), (1000, 0), 100, 100)
assert (size, cols, rows) == (10.0, 100, 101)A hint
The old diagonal has math.hypot(width, height) pixels, 1,000 for 800 by 600. The new diagonal is math.hypot(x2 - x1, y2 - y1) in the new CRS's units, so one pixel is the second divided by the first.
Results of the sample tests
| Test | Result | Details |
|---|
What your code printed
The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.
Practice that adapts to you
Practise this lesson
6 questions about 6 ideas of this lesson. Each answer is explained at once, and what you get right is kept in this browser.
Interview questions
Three longer questions start from this lesson, from a class that appears on every riverbank to areas summed in degrees; the bank gives what each tests and a first step:
-
Fresher to mid
After a land-cover map was reprojected and saved for a web map, every riverbank shows a thin strip of "town", a class that the original map never had there. What happened, and how do you produce the map again correctly?
-
Mid to senior
You combine a 10 m NDVI raster with a 30 m elevation model of the same area, both in one UTM zone, and slope never lines up with field edges. How do you check the two grids, and how do you combine them correctly?
-
Mid to senior
A colleague computes forest area for every district of India from a 10 m land-cover raster reprojected to EPSG:4326, counting forest pixels and multiplying by 100 m². Which errors do you expect, and how would you compute the areas instead?
Flashcards · 6 cards
Recall it, then check
Say the answer to yourself before you look: recalling it is what makes it stick.
-
How does each of the five common resampling methods fill a new cell?
Nearest copies the old cell under the new centre; bilinear weighs the four nearest old centres by distance; cubic fits a smooth curve through sixteen; average takes the mean of the old cells inside, leaving nodata out; mode takes the commonest value among them.
-
Which resampling methods suit a land-cover raster, and why?
Nearest when the cells keep their size and mode when they grow, because class codes are names and must never be interpolated. Bilinear turned 198 cells between water and soil into class 2, a town that was not there.
-
Which resampling methods suit continuous values such as heights or NDVI?
Bilinear or cubic when the cells keep their size or shrink, and average when they grow, so each new cell gets the mean of the old ones inside it; nearest may not suit continuous data.
-
How does GDAL choose the grid of a reprojected raster?
It transforms about 20 points along each edge of the old raster and takes the box round them as the new extent, then chooses one square pixel size such that the diagonal from corner to corner holds as many pixels as before.
-
Is a pixel that is square in degrees also square on the ground?
No. Reprojected to EPSG:4326, the lesson raster got pixels of 9.3189e-05°, which at 18.47° N are 9.84 m from west to east and 10.31 m from south to north: degrees of longitude get shorter with cos(latitude) while degrees of latitude barely change.
-
How is each pixel of a reprojected raster filled?
By working backwards: the new pixel's centre is transformed into the old CRS and the old raster is resampled there. Corners of the new grid that the old raster does not reach get nodata, about 6,505 pixels in the lesson's case.
Key takeaways
- Every resampling decides new values: nearest copies one cell, bilinear and cubic interpolate, average takes the mean and mode the commonest value.
- Codes must never be interpolated, so classes take nearest or mode, while measurements take bilinear or cubic for finer cells and average for coarser ones.
- Reprojecting rasters builds a new grid first, its extent from the transformed edges and one square pixel size from the diagonal, then fills each pixel by working backwards to the old grid.
- A grid that is square in degrees is not square on the ground: the lesson raster’s 9.3189e-05° pixels are 9.84 m by 10.31 m.
Terms in this lesson
- resampling
- Computing the cells of a new grid from an old one, whenever the cell size, the grid's position or the CRS changes. Each method, such as nearest, bilinear, cubic, average or mode, decides the new values in its own way.
- warp (raster reprojection)
- Putting a raster onto a new grid, usually in another CRS: the new grid is chosen first, then each of its pixels is filled by finding its centre on the old raster and resampling there.
References
- gdalwarp: image reprojection and warping utility (Open Source Geospatial Foundation (GDAL))
- GDAL algorithms C API: GDALSuggestedWarpOutput (Open Source Geospatial Foundation (GDAL))
- COG: Cloud Optimized GeoTIFF generator (Open Source Geospatial Foundation (GDAL))
- Resampling (rasterio documentation) (rasterio developers)
- Reprojection (rasterio documentation) (rasterio developers)
- Transverse Mercator with an accuracy of a few nanometers (Karney) (Journal of Geodesy)
- Cubic convolution interpolation for digital image processing (Keys) (IEEE Transactions on Acoustics, Speech, and Signal Processing)
Report a problem with this lesson
Kept only in this browser. Your Learn progress