Skip to content

fill_d8(): numpy and cupy backends have no memory guard #1334

Description

@brendancol

Description

fill_d8() on the numpy and cupy backends has no memory guard.

_fill_tile_kernel (xrspatial/hydro/fill_d8.py:48) plus the public dispatcher (line 485) allocate, beyond the caller's input:

  • dem_f64 = data.astype(np.float64) -> 8 B/pixel (always copies)
  • ring: np.full((h+2, w+2), ...) -> ~8 B/pixel
  • fill = np.empty((h, w), float64) inside the kernel -> 8 B/pixel (returned as output)
  • with z_limit: extra np.where(...) -> 8 B/pixel

That is ~24 B/pixel of working memory plus the caller's input, or ~32 B/pixel with z_limit. A 50000x50000 numpy DEM asks for ~60 GB of host memory before anything errors out.

_fill_cupy (line 211) is the same shape on the device:

  • dem_f64 = dem_data.astype(cp.float64) -> 8 B/pixel
  • fill = cp.empty((H, W), float64) -> 8 B/pixel
  • cp.where(fill > 1e307, dem_f64, fill) returns a new array -> 8 B/pixel (output)
  • with z_limit: extra cp.where(...) -> 8 B/pixel

~24 B/pixel of GPU memory, ~32 B/pixel with z_limit, no check.

The dask paths (_fill_dask_iterative, _fill_dask_cupy) are bounded per-tile by the user's chunk size and stay safe.

Same guard pattern was added in flow_accumulation (#1318/#1319), sieve (#1296), kde (#1287), resample (#1295), focal (#1286), geodesic (#1283), mahalanobis (#1288), true_color (#1291), diffuse (#1267), erode (#1275), emerging_hotspots (#1274), dasymetric (#1261), sky_view_factor (#1299), surface_distance (#1303).

fill_d8 is the depression-filling preprocessing step most hydrology pipelines run first on the full DEM, so it needs the same guard.

Expected behavior

fill_d8() raises MemoryError with a clear message on the eager numpy and cupy backends when the projected working set exceeds available memory. Dask paths skip the guard.

Proposed fix

Add _available_memory_bytes(), _available_gpu_memory_bytes(), _check_memory(rows, cols), and _check_gpu_memory(rows, cols) helpers (32 B/pixel CPU budget, 32 B/pixel GPU budget, 50% threshold). Call them from the public fill_d8() dispatcher before the float64 cast and ring allocation on the numpy and cupy paths.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't workinghigh-priorityoomOut-of-memory risk with large datasets

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions