Skip to content

Resample interpolation uses edge-aligned coordinates instead of block-centered #1202

@brendancol

Description

@brendancol

Describe the bug

resample() with interpolation methods (nearest, bilinear, cubic) uses scipy.ndimage.zoom's coordinate formula, which maps output pixel o to input pixel o * (N_in - 1) / (N_out - 1). This pins the first output pixel to the first input pixel and the last to the last.

The output coordinate metadata (_new_coords) places pixels at block centers, which is the normal convention for raster data. So there's a spatial mismatch between where the interpolation actually samples and where the coordinates say the data is.

The aggregation path (average, min, max, median, mode) uses block-centered sampling and doesn't have this problem.

Reproducer

With a linear gradient where value equals x-coordinate, a correct resample should produce output values that match their coordinate labels:

import numpy as np
from xrspatial.resample import resample
from xrspatial.tests.general_checks import create_test_raster

data = np.tile(np.arange(4, dtype=np.float32), (4, 1))
agg = create_test_raster(data, attrs={'res': (1.0, 1.0)}, dims=['y', 'x'])

out = resample(agg, scale_factor=0.5, method='bilinear')
print(out.x.values)    # [0.5, 2.5]  -- correct block centers
print(out.values[0])   # [0.0, 3.0]  -- wrong, should be [0.5, 2.5]

The interpolation samples at input positions 0.0 and 3.0 (edge-aligned), but the coordinates say the pixels are at 0.5 and 2.5 (block-centered). Half-pixel shift at the edges.

Expected behavior

Interpolation should sample at positions matching the output coordinate metadata. The correct formula is input_pixel = (out_pixel + 0.5) * (N_in / N_out) - 0.5.

Affected code paths

  • _run_numpy -- uses scipy.ndimage.zoom directly
  • _run_cupy -- uses cupyx.scipy.ndimage.zoom directly
  • _interp_block_np -- dask+numpy path, uses map_coordinates with the zoom formula
  • _interp_block_cupy -- dask+cupy path, same

Fix

Replace scipy.ndimage.zoom in the non-dask paths with map_coordinates using the block-centered formula. Update the dask block functions to use the same formula.

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't working

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions