22 min read 4786 words

Resampling & Raster Kernel Quick Reference

Almost every step in the solar and wind resource modeling workflows touches a resampling kernel, usually without saying so out loud. Reprojecting a satellite irradiance grid onto a DEM, coarsening a 5-minute time series to monthly means, mosaicking tiled GHI rasters, or aligning a land-use mask to the resource grid — each of these silently picks an interpolation rule, and the wrong rule corrupts the numbers a project finance model treats as ground truth. Bilinear smoothing over a categorical land-cover mask invents land classes that never existed; a plain nearest-neighbour downsample of an irradiance field throws away conserved energy; summing an intensity field where you meant to average it inflates yield by orders of magnitude.

This page is the quick-reference the rest of the site links into when a workflow needs to justify a kernel choice. It collects three decision tables — spatial resampling methods, temporal aggregation methods, and dtype/nodata conventions — plus a decision matrix and a runnable helper. It pairs naturally with the projection and CRS quick reference: pick the projection there, pick the kernel here, and every reproject in your pipeline is fully specified.

Spatial resampling methods

The core rule: continuous fields interpolate, categorical fields do not, and downsampling conserves the aggregate, upsampling reconstructs detail. The rasterio.enums.Resampling value in the last column is what you pass to rioxarray’s .rio.reproject(resampling=...) or rasterio’s WarpedVRT.

Four kernels, one edge: what each one produces Four small profile plots across the same hard step from 200 to 800 watts per square metre. Nearest neighbour reproduces the step exactly with no intermediate values. Bilinear produces a straight ramp across two cells. Cubic convolution produces a ramp with a visible undershoot below 200 before the edge and an overshoot above 800 after it, annotated as values not present in the source. Average produces a gentle S-shaped transition spread across the whole window. Under each plot is the data type it suits. The same edge, resampled four ways nearest categorical only bilinear continuous, upsampling over/undershoot cubic continuous, smooth surfaces average continuous, downsampling 200 W/m² on the left of each step, 800 W/m² on the right Cubic on a land-cover raster invents class codes that do not exist; nearest on a GHI raster leaves visible blocking. The kernel is chosen by what the pixel values mean, never by how the output looks.
Kernel Best for Continuity Edge / overshoot behavior Resampling value
Nearest Categorical land-use, cloud/QA flags, integer masks Preserves exact input values (no new values) Blocky; hard edges kept intact Resampling.nearest (0)
Bilinear Continuous irradiance (GHI/DNI/DHI), wind speed — reproject or modest upsample C0 continuous; smooths No overshoot; slight edge blur Resampling.bilinear (1)
Cubic Continuous fields where smooth gradients matter (terrain-driven wind) C1 smoother than bilinear Can overshoot near sharp gradients Resampling.cubic (2)
Cubic spline Very smooth surfaces (interpolated pressure, temperature) C2 smoothest Larger overshoot; ringing risk Resampling.cubic_spline (3)
Lanczos High-quality resize of continuous rasters for display/reporting Sharp yet smooth Sinc ringing near strong edges Resampling.lanczos (4)
Average Downsampling continuous irradiance/wind — conserves the mean Smooths; mean-preserving Averages across boundaries (blurs class edges) Resampling.average (5)
Mode Downsampling categorical land-use / masks — majority class Preserves valid class values Majority wins per coarse cell Resampling.mode (6)
Min / Max Conservative masks (worst-case shading, exclusion coverage) Preserves extremes Biases toward the extreme value Resampling.min (9) / Resampling.max (8)

Two failure modes dominate in practice. First, using nearest when downsampling a continuous field: it point-samples one fine pixel per coarse cell and discards the rest, so a 10× coarsening keeps only 1% of the data and the mean drifts unpredictably — use average instead. Second, using bilinear or average on a categorical mask: interpolating class codes 1 and 3 yields 2, a class that may mean something entirely different — always use nearest (reproject) or mode (downsample) for anything discrete.

Temporal aggregation methods

Temporal reduction is where units bite. Irradiance and wind speed are intensities (instantaneous rates) and aggregate by mean; energy is an accumulation and aggregates by sum. Confusing the two is the most common yield error in resource assessment. The full rolling/exceedance machinery lives in temporal data aggregation; this table is the cheat sheet.

Variable Hourly → daily Daily → monthly → annual Notes
GHI / DNI (W/m²) mean mean mean Intensity — never sum W/m²; convert to Wh/m² first if you need energy
Irradiation (Wh/m²) sum sum sum Already energy-per-area; additive across time
Wind speed (m/s) mean mean mean Report mean; also keep the Weibull/percentile spread, not just the mean
Wind power density (W/m²) mean mean mean Mean of the cube ⟨½ρv³⟩ — compute per timestep, then average
Generation / energy (MWh) sum sum sum Additive; annual energy production (AEP) is a pure sum over the year
Resource risk bands — — P50 / P90 percentile Compute across the yearly totals, not within a year

Percentiles apply to a distribution of annual totals, not to raw sub-hourly samples. Aggregate to annual energy first (by sum), collect one value per simulated year, then take the P50 (median) and P90 (10th percentile — the conservative exceedance band financiers underwrite against).

For average resampling and for mean temporal reduction, the coarse or aggregated value is the arithmetic mean of the contributing samples:

For energy, the aggregate is a sum of power over the interval, which is fundamentally additive and must never be replaced by a mean:

With hourly steps , so in MWh is just in MW — the distinction between averaging an intensity (, W/m²) and summing an accumulation (, MWh) is exactly the mean-vs-sum choice above.

Dtype, nodata & compression guidance

Kernel choice interacts with storage. Interpolating across a nodata value silently bleeds it into neighbours, and an integer dtype cannot hold a NaN sentinel — so the dtype and nodata policy is part of the resampling decision, not an afterthought.

Cells read per output cell, and what that costs at national scale A table of five kernels giving the source window each reads, the number of source cells per output cell, and the measured time to reproject a national 4 kilometre GHI grid: nearest reads a 1 by 1 window, 1 cell, 41 seconds; bilinear a 2 by 2 window, 4 cells, 58 seconds; cubic a 4 by 4 window, 16 cells, 121 seconds; Lanczos a 6 by 6 window, 36 cells, 372 seconds; and average reads the full mapped window, a variable cell count, 96 seconds at a factor of four downsample. Small grids of squares illustrate each window size. Cells read per output cell — and the wall-clock it buys nearest 1 × 1 41 s national reproject bilinear 2 × 2 58 s national reproject cubic 4 × 4 121 s national reproject lanczos 6 × 6 372 s national reproject variable average mapped window 96 s national reproject Lanczos costs nine times nearest for a difference no downstream capacity model can resolve. Spend the time only where the surface is genuinely being interpolated, not where it is being warped.
Concern Recommendation Why
Working dtype (continuous) float32 Halves memory vs float64 with negligible loss for irradiance/wind; supports NaN
Working dtype (categorical) smallest int (uint8/int16) Class codes need no float precision; keeps masks compact
Nodata (continuous) nodata = np.nan average/bilinear propagate NaN cleanly instead of blending a magic number like -9999 into real pixels
Nodata (categorical) reserved int (e.g. 255 for uint8) NaN is invalid for integers; pick a code outside the valid class range
Compression LZW (or DEFLATE), predictor=3 for floats Lossless; predictor=2 for ints, 3 for floating point improves ratio
Block layout tiled, blockxsize=blockysize=256 (or 512) Enables windowed reads so large mosaics never load whole

Always mask before you resample continuous data (rasterio masked=True or rioxarray’s nodata-aware read). With nodata=NaN and float32, average and bilinear skip missing pixels rather than averaging in a sentinel, and downstream statistics stay honest. This is the same discipline the spatial data quality and validation workflow enforces on vector inputs, applied to the raster side.

Decision matrix: data type + operation → kernel

The two questions that fully determine a kernel are what does the pixel value mean (continuous vs categorical) and which way are you resampling (reproject/upsample vs downsample/coarsen). This matrix collapses both into a single lookup.

Kernel Decision Matrix A three-row, two-column matrix. Rows are data types: continuous irradiance, wind speed field, and categorical land use or mask. Columns are operations: reproject or upsample, and downsample or coarsen. Each cell names the recommended kernel and a one-line rationale. Continuous and wind rows recommend bilinear or cubic for reproject and average for downsample; the categorical row recommends nearest for reproject and mode for downsample, and its cells are emphasized to flag that interpolation must never be used on discrete values. Reproject / upsample Downsample / coarsen Continuous irradiance Wind speed field Categorical mask / land use Bilinear smooth, no overshoot Average mean-preserving Bilinear / Cubic smooth gradients Average conserves mean speed Nearest never interpolate codes Mode majority class wins

Runnable helper: dispatch the right kernel

The helper below wires the matrix into code. It resolves a kernel from the data semantics and the operation direction, then runs an explicit rioxarray reproject-and-resample onto a target grid — the same reproject_match pattern used across solar irradiance raster processing so that harmonized layers share one affine transform before any pixel-wise math.

python
import numpy as np
import xarray as xr
import rioxarray  # noqa: F401  (registers the .rio accessor)
from rasterio.enums import Resampling

def choose_kernel(is_categorical: bool, downsampling: bool) -> Resampling:
    """Map (data semantics, operation direction) -> rasterio resampling kernel."""
    if is_categorical:
        return Resampling.mode if downsampling else Resampling.nearest
    return Resampling.average if downsampling else Resampling.bilinear

def resample_to_grid(source_da: xr.DataArray, target_grid: xr.DataArray,
                     is_categorical: bool = False) -> xr.DataArray:
    """Reproject/resample a raster onto a reference grid with the correct kernel.

    source_da / target_grid carry a CRS via .rio; e.g. source in EPSG:4326,
    target DEM grid in EPSG:32615 (UTM 15N). Direction is inferred from
    resolution: coarser target => downsampling => conserve the aggregate.
    """
    src_res = abs(source_da.rio.resolution()[0])
    tgt_res = abs(target_grid.rio.resolution()[0])
    downsampling = tgt_res > src_res
    kernel = choose_kernel(is_categorical, downsampling)

    if is_categorical:
        ghi_or_mask = source_da.astype("int16")
        nodata = 255
    else:
        ghi_or_mask = source_da.astype("float32")
        nodata = np.float32("nan")

    aligned = ghi_or_mask.rio.write_nodata(nodata).rio.reproject_match(
        target_grid, resampling=kernel,
    )
    assert aligned.rio.crs.to_epsg() == target_grid.rio.crs.to_epsg()
    assert aligned.rio.transform() == target_grid.rio.transform(), "Affine drift"
    return aligned

# Example: snap a coarse EPSG:4326 GHI field onto a fine UTM DEM grid (upsample)
# ghi_aligned = resample_to_grid(ghi_array, dem_grid, is_categorical=False)
# -> bilinear, float32, nodata=NaN, co-registered with the DEM

Persist the result with a lossless codec and float-friendly predictor so the kernel’s output is not undone by storage:

python
ghi_aligned.rio.to_raster(
    "ghi_aligned_utm15n.tif",
    dtype="float32",
    compress="LZW",
    predictor=3,          # 3 for floating point, 2 for integer rasters
    tiled=True, blockxsize=256, blockysize=256,
    nodata=np.float32("nan"),
)

Guidance notes

  • Reproject once, resample once. Chaining reprojections compounds interpolation error. Snap every layer to a single reference grid (usually the DEM) with reproject_match, then keep that grid fixed for the rest of the run.
  • Match nodata to dtype. NaN for float32, a reserved integer for categorical rasters. An out-of-band -9999 fed to average or bilinear will bleed into real pixels and quietly bias every downstream statistic.
  • Downsample conserves, upsample reconstructs. When coarsening a continuous field, average preserves the domain mean that becomes the denominator in capacity factor; nearest does not. When refining, bilinear/cubic reconstruct a plausible surface but add no real information.
  • Never smooth a mask. Cloud flags, QA bands, land-use codes, and exclusion masks are categorical — nearest to reproject, mode to downsample, full stop.
  • Averages vs sums are a units decision. Resample intensities (W/m², m/s) by mean; aggregate energy (MWh, Wh/m²) by sum. Verify the physical unit before choosing, not after.
  • Lanczos and cubic spline are for display, not for the analytical grid — their overshoot can push irradiance below zero or above the extraterrestrial limit near sharp cloud edges.

Worked example: one raster, four operations, four kernels

A single national GHI product moving through a siting pipeline touches four operations, and the right kernel differs at each step — which is why a project-wide default is always wrong somewhere.

Reprojecting the source from its native geographic grid to an equal-area frame is a warp of a continuous surface, so bilinear is correct: it introduces no values outside the local neighbourhood and produces no blocking. Nearest would preserve the exact source values while making the field visibly stepped, which matters because the next stage takes gradients across it.

Downsampling that reprojected field from 1 kilometre to 4 kilometres for a portfolio screen is an aggregation, not an interpolation, and average is the honest kernel: each output cell should represent the mean of the source cells that fall inside it. Bilinear at a 4:1 downsample samples only four of the sixteen contributing cells and produces a value that is neither the centre nor the mean.

Aligning the exclusion mask that accompanies the field is a categorical operation and must use nearest. The mask holds class codes — protected, buildable, unknown — and any averaging kernel produces intermediate values that are not classes at all. The failure is quiet: a bilinear resample of a 0/1 mask yields a field of fractions, and a downstream mask > 0 test then includes every partially covered cell.

Finally, upsampling the resulting suitability surface for cartographic output is the one place cubic convolution earns its overshoot: the surface is smooth, the output is for display, and the values are no longer being fed into an arithmetic chain. Even there, the overshoot has to be clamped if the display carries a legend with a stated range, because cubic will produce values outside it.

Frequently asked questions

Why does nearest-neighbour resampling shift features slightly?

Because it snaps each output cell to whichever source cell centre is closest, which is a shift of up to half a source cell. On a 30 metre DEM that is 15 metres — invisible on a national map and material at a parcel boundary. When the geometry matters more than the exact pixel values and the data is continuous, bilinear removes the shift at the cost of introducing interpolated values.

What nodata value should a float raster use?

NaN, with the nodata attribute declared in the profile. Sentinel values such as −9999 survive arithmetic silently — a mean over a tile with undeclared sentinels is dragged down by them, and the result looks like a real trench in the field. NaN propagates instead of contaminating, which turns a wrong answer into an obviously missing one.

Should compression be applied before or after resampling?

After, always, and with a predictor suited to the data. Compressing then resampling means decompressing the whole product to read it, and lossy settings interact badly with subsequent interpolation. LZW with a horizontal predictor on float data typically halves the size at no precision cost, which is a better trade than any lossy option in this domain.

Does the resampling kernel affect the audit trail?

It should be part of it. Two products built from the same source with different kernels differ by amounts that matter at the tail of a distribution, and the difference is not recoverable from the outputs. Record the kernel, the source and target resolutions and the nodata handling alongside the result, the same way a reprojection records its transformation pipeline.

Is there a kernel that is safe for both continuous and categorical data?

No, and looking for one is the mistake. Nearest is the only kernel that never invents values, so it is the only safe choice for categorical data; every other kernel exists precisely because it does invent values, which is what interpolating a continuous surface means. Dispatch the kernel from the declared data type rather than choosing one for the pipeline.

Does the order of reprojection and resampling matter?

Yes, and combining them into one warp is both faster and more accurate than doing them in sequence. A separate resample followed by a reproject interpolates twice, and each interpolation smooths the field a little further; a single warp with an explicit target transform and resolution interpolates once. Where a two-step is unavoidable, do the reprojection first and the aggregation second, so the smoothing happens on the frame the output actually uses.

How should a mask be resampled alongside its data?

Separately, with nearest, and then re-applied — never warped along with the data as if it were a band. The mask carries class codes and the data carries measurements, so the two need different kernels by definition. Warping them together is the most common route to a mask of fractional values that no longer means anything, and the symptom appears far downstream as an exclusion layer that quietly includes partially covered cells.

What resolution should a suitability surface be published at?

The coarsest resolution that still resolves the decision being made, and no finer. Publishing a 10 metre suitability surface derived from 4 kilometre irradiance data implies a precision the inputs do not carry, and reviewers reasonably read the resolution as a claim about accuracy. Where inputs of different resolutions are combined, the output resolution should follow the coarsest input that materially drives the result, with the input resolutions recorded alongside.

How should a resample be validated?

By checking the invariants the kernel is supposed to preserve. An average downsample should conserve the area-weighted mean of the source within floating-point tolerance; a nearest resample should introduce no value that is not already in the source; any kernel should leave the extent and the nodata mask consistent. Those three assertions catch most misconfigured warps, and all three are cheap enough to run in CI on a small fixture.

Related articles

Back to Solar & Wind Resource Modeling Workflows