Skip to content

clip_domain raises ValueError (broadcast mismatch) for extents strictly inside the domain #557

Description

@martinguthrie93

Bug description

pysteps.utils.dimension.clip_domain raises a ValueError when the requested clipping extent lies strictly inside the source domain (i.e. when none of the extent edges coincides with the source domain boundary).

This is a core GIS utility used to clip a precipitation field to a geographic bounding box, so the failure affects any workflow that crops a radar/NWP domain to an interior region of interest.

How to reproduce

import numpy as np
from pysteps.utils import dimension

# 10x10 grid, 10-unit pixels, domain spanning 0..100 in both directions
R = np.arange(100, dtype=float).reshape(10, 10)
metadata = {
    "x1": 0.0, "x2": 100.0, "y1": 0.0, "y2": 100.0,
    "xpixelsize": 10.0, "ypixelsize": 10.0,
    "yorigin": "upper", "zerovalue": 0.0,
}

# Clip to an interior 40x40 box that does not touch the domain edges
dimension.clip_domain(R, metadata, extent=(30, 70, 30, 70))

Result:

ValueError: could not broadcast input array from shape (1,1,4,4) into shape (1,1,5,5)

Root cause

clip_domain computes the source window and the destination window independently:

  • the destination size is derived from the requested extent: dim_x_ = int((right_ - left_) / xpixelsize);
  • the source indices are found with a strict pixel-center test np.where((x_coord < right_) & (x_coord > left_)), evaluated against the requested extent, while the destination indices are found with the same kind of test against the full source domain.

These two windows are only guaranteed to have the same length when a clip edge happens to line up with the source-domain boundary. When the extent is strictly interior, the strict </> comparison drops a pixel on one side, the two windows differ in size by one, and the assignment

R_[..., idx_y_[0]:idx_y_[-1]+1, idx_x_[0]:idx_x_[-1]+1] = R[..., idx_y[0]:idx_y[-1]+1, idx_x[0]:idx_x[-1]+1]

fails to broadcast.

Why the tests miss it

The existing test_clip_domain cases only use extent=(2, 4, 2, 4) on a 0..4 domain, i.e. the clip's right/top edge coincides with the domain boundary. In that special case the two windows happen to match, so the bug is never triggered.

Proposed fix

Because clip_domain preserves the pixel size, clipping is an integer-pixel window operation. Computing a single pixel offset between the clipped grid and the source array guarantees the source and destination windows always have matching sizes, for both yorigin conventions. As a bonus, extents that extend beyond the source domain can then be padded with the zero value instead of raising.

I will open a PR with this fix and regression tests.

Environment

  • pysteps: master
  • numpy: 2.3.5

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Fields

    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions