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
Bug description
pysteps.utils.dimension.clip_domainraises aValueErrorwhen 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
Result:
Root cause
clip_domaincomputes the source window and the destination window independently:dim_x_ = int((right_ - left_) / xpixelsize);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 assignmentfails to broadcast.
Why the tests miss it
The existing
test_clip_domaincases only useextent=(2, 4, 2, 4)on a0..4domain, 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_domainpreserves 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 bothyoriginconventions. 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
master