diff --git a/R/CFAxisCharacter.R b/R/CFAxisCharacter.R
index c2856af..0a3bf8d 100644
--- a/R/CFAxisCharacter.R
+++ b/R/CFAxisCharacter.R
@@ -130,8 +130,9 @@ CFAxisCharacter <- R6::R6Class("CFAxisCharacter",
       else {
         rng <- range(rng)
         if (self$has_resource) {
-          ax <- CFAxisCharacter$new(private$.NCobj, group = group, start = private$.NC_map$start + rng[1L] -1L,
-                                    count = rng[2L] - rng[1L] + 1L, attributes = self$attributes)
+          nc <- private$to_nc_indices(rng[1L], rng[2L] - rng[1L] + 1L)
+          ax <- CFAxisCharacter$new(private$.NCobj, group = group, start = nc$start,
+                                    count = nc$count, attributes = self$attributes)
           if (nzchar(name))
             ax$name <- name
         } else {
diff --git a/R/CFAxisLatitude.R b/R/CFAxisLatitude.R
index 21d7523..608b78e 100644
--- a/R/CFAxisLatitude.R
+++ b/R/CFAxisLatitude.R
@@ -95,8 +95,9 @@ CFAxisLatitude <- R6::R6Class("CFAxisLatitude",
       else {
         rng <- range(rng)
         if (self$has_resource) {
-          ax <- CFAxisLatitude$new(private$.NCobj, group = group, start = private$.NC_map$start + rng[1L] - 1L,
-                                   count = rng[2L] - rng[1L] + 1L, attributes = self$attributes)
+          nc <- private$to_nc_indices(rng[1L], rng[2L] - rng[1L] + 1L)
+          ax <- CFAxisLatitude$new(private$.NCobj, group = group, start = nc$start,
+                                   count = nc$count, attributes = self$attributes)
           if (nzchar(name))
             ax$name <- name
         } else {
diff --git a/R/CFAxisLongitude.R b/R/CFAxisLongitude.R
index cfd819c..b527fd8 100644
--- a/R/CFAxisLongitude.R
+++ b/R/CFAxisLongitude.R
@@ -97,8 +97,9 @@ CFAxisLongitude <- R6::R6Class("CFAxisLongitude",
       else {
         rng <- range(rng)
         if (self$has_resource) {
-          ax <- CFAxisLongitude$new(private$.NCobj, group = group, start = private$.NC_map$start + rng[1L] - 1L,
-                                    count = rng[2L] - rng[1L] + 1L, attributes = self$attributes)
+          nc <- private$to_nc_indices(rng[1L], rng[2L] - rng[1L] + 1L)
+          ax <- CFAxisLongitude$new(private$.NCobj, group = group, start = nc$start,
+                                    count = nc$count, attributes = self$attributes)
           if (nzchar(name))
             ax$name <- name
         } else {
diff --git a/R/CFAxisNumeric.R b/R/CFAxisNumeric.R
index fbbecfa..f0023d2 100644
--- a/R/CFAxisNumeric.R
+++ b/R/CFAxisNumeric.R
@@ -285,8 +285,9 @@ CFAxisNumeric <- R6::R6Class("CFAxisNumeric",
       else {
         rng <- range(rng)
         if (self$has_resource) {
-          ax <- CFAxisNumeric$new(private$.NCobj, group = group, start = private$.NC_map$start + rng[1L] -1L,
-                                  count = rng[2L] - rng[1L] + 1L, orientation = private$.orient,
+          nc <- private$to_nc_indices(rng[1L], rng[2L] - rng[1L] + 1L)
+          ax <- CFAxisNumeric$new(private$.NCobj, group = group, start = nc$start,
+                                  count = nc$count, orientation = private$.orient,
                                   attributes = self$attributes)
           if (nzchar(name))
             ax$name <- name
diff --git a/R/CFAxisTime.R b/R/CFAxisTime.R
index ca64c3b..16d0e2e 100644
--- a/R/CFAxisTime.R
+++ b/R/CFAxisTime.R
@@ -379,9 +379,9 @@ CFAxisTime <- R6::R6Class("CFAxisTime",
       else {
         rng <- as.integer(range(rng))
         if (self$has_resource) {
+          nc <- private$to_nc_indices(rng[1L], rng[2L] - rng[1L] + 1L)
           ax <- CFAxisTime$new(private$.NCobj, group = group, values = self$values[rng[1L]:rng[2L]],
-                               start = private$.NC_map$start + rng[1L] - 1L,
-                               count = rng[2L] - rng[1L] + 1L, attributes = self$attributes)
+                               start = nc$start, count = nc$count, attributes = self$attributes)
           if (nzchar(name))
             ax$name <- name
         } else {
diff --git a/R/CFAxisVertical.R b/R/CFAxisVertical.R
index a150a31..6745c22 100644
--- a/R/CFAxisVertical.R
+++ b/R/CFAxisVertical.R
@@ -426,8 +426,9 @@ CFAxisVertical <- R6::R6Class("CFAxisVertical",
       else {
         rng <- range(rng)
         if (self$has_resource) {
-          ax <- CFAxisVertical$new(private$.NCobj, group = group, start = private$.NC_map$start + rng[1L] - 1L,
-                                   count = rng[2L] - rng[1L] + 1L, attributes = self$attributes)
+          nc <- private$to_nc_indices(rng[1L], rng[2L] - rng[1L] + 1L)
+          ax <- CFAxisVertical$new(private$.NCobj, group = group, start = nc$start,
+                                   count = nc$count, attributes = self$attributes)
           if (nzchar(name))
             ax$name <- name
         } else {
diff --git a/R/CFData.R b/R/CFData.R
index d6b02f4..ed442ed 100644
--- a/R/CFData.R
+++ b/R/CFData.R
@@ -64,8 +64,9 @@ CFData <- R6::R6Class("CFData",
     },
 
     # Sanitize the start and count values. Called during initialization and by
-    # read_chunk(). NAs are converted to numbers and values have to agree with
-    # .dims. Returns the sanitized start and count vectors as a list.
+    # read_chunk()/read_window(). NAs are converted to numbers and values have
+    # to agree with .dims. Returns the sanitized start and count vectors as a
+    # list.
     check_start_count = function(start, count) {
       d <- private$.dims
       len <- length(d)
@@ -76,9 +77,9 @@ CFData <- R6::R6Class("CFData",
       else if (length(start) == len) {
         start[is.na(start)] <- 1L
         if (any(start > d))
-          stop("Start values cannot be larger than the dimensions of the data.", call. = FALSE) # nocov
+          stop("Start values cannot be larger than the dimensions of the data", call. = FALSE) # nocov
       } else
-        stop("`start` vector is not the length of the object dimensions.", call. = FALSE)
+        stop("`start` vector is not the length of the object dimensions", call. = FALSE)
 
       if (length(count) == 1L && is.na(count))
         count <- d - start + 1L
@@ -86,13 +87,77 @@ CFData <- R6::R6Class("CFData",
         ndx <- which(is.na(count))
         count[ndx] <- d[ndx] - start[ndx] + 1L
         if (any(count > d - start + 1L))
-          stop("Count values cannot extend beyond the dimensions of the data.", call. = FALSE) # nocov
+          stop("Count values cannot extend beyond the dimensions of the data", call. = FALSE) # nocov
       } else
-        stop("`count` vector is not the length of the object dimensions.", call. = FALSE)
+        stop("`count` vector is not the length of the object dimensions", call. = FALSE)
 
       list(start = start, count = count)
     },
 
+    # Translate local (validated, .dims-relative) start/count into the
+    # file-absolute (NC-space) indices NCobj$get_data()/write_data()
+    # require, applying this object's .NC_map offset. Only meaningful when
+    # .NC_map is non-empty; callers with an empty .NC_map are already in
+    # NC-space (nothing to translate) and must not call this.
+    to_nc_indices = function(local_start, local_count) {
+      len <- length(private$.NC_map$start)
+      list(start = private$.NC_map$start + local_start[seq_len(len)] - 1L,
+           count = local_count[seq_len(len)])
+    },
+
+    # Warn if reading `count` elements at the object's current data type
+    # would allocate more than CF.options$memory_cell_limit bytes. For
+    # packed data, RNetCDF's unpack=TRUE typically widens to double
+    # regardless of the on-disk (packed) type -- fitnum may narrow it back
+    # down afterwards, but that's not knowable in advance, so packed
+    # variables are sized conservatively at 8 bytes/element here. Returns
+    # the estimated byte count, invisibly, or NA if the type is unknown.
+    check_memory_limit = function(count) {
+      itemsize <- if (!is.null(private$.NCobj) && private$.NCobj$is_packed) 8L
+      else .nc_type_size(private$.data_type)
+      if (is.na(itemsize)) return(invisible(NA_real_))
+
+      bytes <- prod(count) * itemsize
+      if (bytes > CF.options$memory_cell_limit)
+        warning(sprintf(
+          "Reading %s elements will allocate approximately %s, exceeding the memory limit (%s). Consider reading a smaller extent, or raise the limit if you have the memory to spare",
+          paste(count, collapse = "x"), .format_bytes(bytes), .format_bytes(CF.options$memory_cell_limit)),
+          call. = FALSE)
+      invisible(bytes)
+    },
+
+    # Core chunk-planned read + reassembly for resolved start/count,
+    # in whatever coordinate space the caller has established is correct for
+    # them (file-relative for a .NC_map read, local for a validated one). No
+    # NA resolution, no bounds-checking against .dims -- that's the caller's
+    # job, done once before calling this.
+    read_window_core = function(start, count, reader) {
+      bytes <- private$check_memory_limit(count)
+
+      if (!is.na(bytes) && bytes <= CF.options$memory_cell_limit)
+        return(reader(start, count))  # fits outright: skip .chunk_plan() entirely
+
+      chunks <- if (!is.null(private$.NCobj)) private$.NCobj$netcdf4$chunksizes else NULL
+      itemsize <- if (!is.null(private$.NCobj) && private$.NCobj$is_packed) 8L
+      else (.nc_type_size(private$.data_type) %||% 8L)
+
+      plan <- .chunk_plan(count, chunks, start, count, itemsize, CF.options$memory_cell_limit)
+      if (length(plan) == 1L)
+        return(reader(plan[[1L]]$start, plan[[1L]]$count))
+
+      out <- NULL
+      for (p in plan) {
+        block <- reader(p$start, p$count)
+        if (is.null(out)) out <- array(vector(storage.mode(block), prod(count)), dim = count)
+        idx <- lapply(seq_along(count), function(d) {
+          off <- p$start[d] - start[d]
+          seq.int(off + 1L, off + p$count[d])
+        })
+        out <- do.call(`[<-`, c(list(out), idx, list(value = block)))
+      }
+      out
+    },
+
     # Set the values of the object. Perform some basic checks when set programmatically.
     # Values may be NULL
     set_values = function(values) {
@@ -127,7 +192,7 @@ CFData <- R6::R6Class("CFData",
       # Set the actual_range attribute for the values
       if (is.null(values))
         self$delete_attribute("actual_range")
-      else if (prod(dim(values)) <= CF.options$memory_cell_limit) {
+      else if (prod(dim(values)) * (.nc_type_size(private$.data_type) %||% 8L) <= CF.options$memory_cell_limit) {
         rng <- suppressWarnings(range(values, na.rm = TRUE))
         if (is.infinite(rng[1L]) || is.na(rng[1L]))
           self$delete_attribute("actual_range")
@@ -165,9 +230,8 @@ CFData <- R6::R6Class("CFData",
           # are trimmed to that length (noting that there may be "scalar" axes
           # in a variable backed by a netCDF resource).
           sc <- private$check_start_count(start, count)
-          len <- length(private$.NC_map$start)
-          start <- private$.NC_map$start + sc$start[1L:len] - 1L
-          private$.NCobj$write_data(d = dt, start = start, count = sc$count[1L:len], ...)
+          nc <- private$to_nc_indices(sc$start, sc$count)
+          private$.NCobj$write_data(d = dt, start = nc$start, count = nc$count, ...)
         } else {
           # .NC_map is an empty list for private$.values being a complete array.
           private$.NCobj$write_data(d = dt, start = NA, count = NA, ...)
@@ -252,7 +316,8 @@ CFData <- R6::R6Class("CFData",
     #' argument `refresh` is `TRUE`. This method will not assess how big the
     #' data is before reading it so there is a chance that memory will be
     #' exhausted. The calling code should check for this possibility and break
-    #' up the reading of data into chunks.
+    #' up the reading of data into chunks using `read_chunk()` or
+    #' `read_window()` - these methods will not cache the data, however.
     #' @param refresh Should the data be read from file if the object is linked?
     #'   This will replace current values, if previously loaded. Default
     #'   `FALSE`.
@@ -262,9 +327,16 @@ CFData <- R6::R6Class("CFData",
     read_data = function(refresh = FALSE) {
       if ((!is.null(private$.NCobj)) && (is.null(private$.values) || refresh)) {
         if (!length(private$.NC_map))
-          private$set_values(private$.NCobj$get_data())
+          # No .NC_map: local space and NC-space coincide, so a full-extent read
+          # against .dims is valid: read directly.
+          private$set_values(
+            private$read_window_core(rep(1L, length(private$.dims)), private$.dims,
+                                     private$.NCobj$get_data))
         else
-          private$set_values(private$.NCobj$get_data(private$.NC_map$start, private$.NC_map$count))
+          # .NC_map is in NC-space: read directly via NCobj$get_data().
+          private$set_values(
+            private$read_window_core(private$.NC_map$start, private$.NC_map$count,
+                                     private$.NCobj$get_data))
       }
       invisible(private$.values)
     },
@@ -273,7 +345,8 @@ CFData <- R6::R6Class("CFData",
     #' `count` vectors. Note that these vectors are relative to any subset of
     #' the data variable that this CF object refers to. The data read by this
     #' method will not be stored in `self` so the calling code must take a
-    #' reference to it.
+    #' reference to it. If the chunk is (potentially) large, use `read_window()`
+    #' for additional safeguards against memory exhaustion.
     #' @param start Vector of indices where to start reading data along the
     #'   dimensions of the array. The vector must be `NA` to read all data,
     #'   otherwise it must agree with the dimensions of the array.
@@ -288,17 +361,30 @@ CFData <- R6::R6Class("CFData",
       if (!length(sc)) return (NULL)
 
       if (!is.null(private$.values)) {
-        # Extract from loaded data
         cll <- paste0("private$.values[", paste(sc$start, ":", sc$start + sc$count - 1L, sep = "", collapse = ", "), "]")
         eval(parse(text = cll))
       } else {
-        # Read from the netCDF resource. .NC_map always refers to the initial
-        # dimensions, so the arguments are trimmed to that length (noting that
-        # there may be "scalar" axes in a variable backed by a netCDF resource).
-        len <- length(private$.NC_map$start)
-        start <- private$.NC_map$start + sc$start[1L:len] - 1L
-        private$.NCobj$get_data(start, sc$count[1L:len])
+        # Read from the netCDF resource: translate validated LOCAL start/count
+        # to NC-space via this object's .NC_map offset, then read directly.
+        nc <- private$to_nc_indices(sc$start, sc$count)
+        private$.NCobj$get_data(nc$start, nc$count)
       }
+    },
+
+    #' @description Read the array for `(start, count)`, transparently splitting
+    #'   the request into budget-bounded, chunk-aligned sub-reads when the
+    #'   window is large, and reassembling into one array via indexed
+    #'   assignment. This materializes the full `(start, count)` window: it
+    #'   bounds each individual read, and warns if the reassembled whole exceeds
+    #'   `CF.options$memory_cell_limit`, but does not itself reduce memory below
+    #'   the size of that window. Use this in preference to `read_chunk()` for
+    #'   any read that might be large.
+    #' @param start,count As for `read_chunk()`.
+    #' @return An array, as `read_chunk()` would return.
+    read_window = function(start, count) {
+      sc <- private$check_start_count(start, count)
+      if (!length(sc)) return(NULL)
+      private$read_window_core(sc$start, sc$count, self$read_chunk)
     }
   ),
   active = list(
diff --git a/R/CFVariable.R b/R/CFVariable.R
index fcb17d0..753ef05 100644
--- a/R/CFVariable.R
+++ b/R/CFVariable.R
@@ -178,6 +178,157 @@ CFVariable <- R6::R6Class("CFVariable",
       list(index = index, X = c(ry[1L], rows), Y = c(rx[1L], cols), aoi = private$.llgrid$aoi, box = dim_index)
     },
 
+    # === Reducer functions === These functions implement base functions sum,
+    # mean, min, max, range over chunks of data such that they can be applied
+    # over variables whose data exceeds available memory. The reducer functions
+    # are used by process_level().
+    .reducer_sum = list(
+      init = function(shape) array(0, dim = shape),
+      update = function(acc, block, tdim, na_rm)
+        acc + .process.data(block, tdim, FUN = sum, na.rm = na_rm)[[1L]],
+      finalize = function(acc, na_rm) list(acc)
+    ),
+
+    .reducer_mean = list(
+      init = function(shape) list(sum = array(0, dim = shape), n = array(0, dim = shape)),
+      update = function(acc, block, tdim, na_rm) {
+        list(sum = acc$sum + .process.data(block, tdim, FUN = sum, na.rm = na_rm)[[1L]],
+             n   = acc$n   + .process.data(block, tdim, FUN = function(x, ...)
+               if (na_rm) sum(!is.na(x)) else length(x))[[1L]])
+      },
+      finalize = function(acc, na_rm) list(acc$sum / acc$n)
+    ),
+
+    # min/max/range: per-block partials come from min()/max() themselves, which
+    # already turn an all-NA block-at-location into Inf/-Inf (with a warning)
+    # under na.rm = TRUE -- that warning is suppressed per block and re-issued
+    # exactly once at finalize(), conditioned on the FINAL accumulator, matching
+    # base R's own single warning for a whole-vector call.
+    .reducer_min = list(
+      init = function(shape) array(Inf, dim = shape),
+      update = function(acc, block, tdim, na_rm)
+        pmin(acc, suppressWarnings(.process.data(block, tdim, FUN = min, na.rm = na_rm)[[1L]]), na.rm = na_rm),
+      finalize = function(acc, na_rm) {
+        if (na_rm && any(is.infinite(acc) & acc > 0))
+          warning("no non-missing arguments to min; returning Inf", call. = FALSE)
+        list(acc)
+      }
+    ),
+
+    .reducer_max = list(
+      init = function(shape) array(-Inf, dim = shape),
+      update = function(acc, block, tdim, na_rm)
+        pmax(acc, suppressWarnings(.process.data(block, tdim, FUN = max, na.rm = na_rm)[[1L]]), na.rm = na_rm),
+      finalize = function(acc, na_rm) {
+        if (na_rm && any(is.infinite(acc) & acc < 0))
+          warning("no non-missing arguments to max; returning -Inf", call. = FALSE)
+        list(acc)
+      }
+    ),
+
+    .reducer_range = list(
+      init = function(shape) list(min = array(Inf, dim = shape), max = array(-Inf, dim = shape)),
+      update = function(acc, block, tdim, na_rm) {
+        list(min = pmin(acc$min, suppressWarnings(.process.data(block, tdim, FUN = min, na.rm = na_rm)[[1L]]), na.rm = na_rm),
+             max = pmax(acc$max, suppressWarnings(.process.data(block, tdim, FUN = max, na.rm = na_rm)[[1L]]), na.rm = na_rm))
+      },
+      finalize = function(acc, na_rm) {
+        if (na_rm) {
+          if (any(is.infinite(acc$min) & acc$min > 0))
+            warning("no non-missing arguments to min; returning Inf", call. = FALSE)
+          if (any(is.infinite(acc$max) & acc$max < 0))
+            warning("no non-missing arguments to max; returning -Inf", call. = FALSE)
+        }
+        list(acc$min, acc$max)  # matches .process.data()'s asplit() order exactly
+      }
+    ),
+
+    # Match by function identity, not name -- a user's own function called
+    # "mean" is never silently intercepted. mean's sum/count decomposition
+    # is only valid for trim = 0 (the default); a trimmed mean isn't
+    # linearly decomposable across blocks, so that case falls through to
+    # the materialized path instead.
+    find_reducer = function(fun, dots) {
+      if (identical(fun, base::sum)) return(private$.reducer_sum)
+      if (identical(fun, base::mean) && (is.null(dots$trim) || dots$trim == 0))
+        return(private$.reducer_mean)
+      if (identical(fun, base::min)) return(private$.reducer_min)
+      if (identical(fun, base::max)) return(private$.reducer_max)
+      if (identical(fun, base::range) && !isTRUE(dots$finite)) return(private$.reducer_range)
+      NULL
+    },
+
+    # Plan blocks for streaming that vary only along a time dimension, at the
+    # full requested extent on every other dimension so every block's
+    # per-location reduction lines up with the same accumulator locations and
+    # can be merged with plain element-wise arithmetic. Deliberately not
+    # .chunk_plan(): that planner is free to also split the other dimensions to
+    # hit budget, which breaks a location-keyed accumulator (blocks would cover
+    # different, non-overlapping spatial regions rather than more time at the
+    # same locations).
+    plan_tdim_only = function(start, count, tdim, tdim_native, itemsize, budget) {
+      other_elems <- prod(as.numeric(count[-tdim]))
+      max_mult <- max(1L, floor(budget / (itemsize * other_elems * tdim_native)))
+      block_tdim <- min(max_mult * tdim_native, count[tdim])
+
+      tstart <- seq(start[tdim], start[tdim] + count[tdim] - 1L, by = block_tdim)
+      lapply(tstart, function(s) {
+        st <- start; st[tdim] <- s
+        ct <- count; ct[tdim] <- min(block_tdim, start[tdim] + count[tdim] - s)
+        list(start = st, count = ct)
+      })
+    },
+
+    # Compute one factor level's contribution from one or more contiguous
+    # tdim runs (already split at disparate-index boundaries by the caller).
+    # Streams via .chunk_plan() when the level exceeds budget and `fun` has
+    # a registered reducer; otherwise materializes and reduces as before.
+    process_level = function(runs, tdim, num_dims, fun, ...) {
+      runs <- lapply(runs, function(r) {
+        sc <- private$check_start_count(r$start, r$count)
+        list(start = sc$start, count = sc$count)
+      })
+
+      itemsize <- .nc_type_size(self$data_type) %||% 8L
+      total_elems <- sum(sapply(runs, function(r) prod(r$count)))
+
+      materialize <- function() {
+        if (length(runs) == 1L)
+          self$read_window(runs[[1L]]$start, runs[[1L]]$count)
+        else
+          abind::abind(lapply(runs, function(r) self$read_window(r$start, r$count)), along = num_dims)
+      }
+
+      limit <- CF.options$memory_cell_limit
+      if (total_elems * itemsize <= limit)
+        return(.process.data(materialize(), tdim, FUN = fun, ...))
+
+      dots <- list(...)
+      reducer <- private$find_reducer(fun, dots)
+      chunks <- if (!is.null(private$.NCobj)) private$.NCobj$netcdf4$chunksizes else NULL
+      tdim_native <- if (!is.null(chunks)) chunks[tdim] else 1L
+      feasible <- !is.null(reducer) && all(sapply(runs, function(r)
+        prod(as.numeric(r$count[-tdim])) * tdim_native * itemsize <= limit))
+      if (feasible) {
+        na_rm <- isTRUE(dots$na.rm)
+        acc <- NULL
+        for (r in runs) {
+          plan <- private$plan_tdim_only(r$start, r$count, tdim, tdim_native, itemsize, limit)
+          for (p in plan) {
+            block <- self$read_chunk(p$start, p$count)
+            if (is.null(acc)) {
+              bd <- dim(block)
+              acc <- reducer$init(if (is.null(bd)) 1L else bd[-tdim])
+            }
+            acc <- reducer$update(acc, block, tdim, na_rm)
+          }
+        }
+        return(reducer$finalize(acc, na_rm))
+      }
+
+      .process.data(materialize(), tdim, FUN = fun, ...)
+    },
+
     # Internal apply/tapply method for this class. If the size of the data
     # variable is below a certain threshold, read the data and process in one
     # go. Otherwise processing goes per factor level. In other words, for each
@@ -187,7 +338,8 @@ CFVariable <- R6::R6Class("CFVariable",
     process_data = function(tdim, fac, fun, ...) {
       if (!is.null(private$.values))
         return(.process.data(self$values, tdim, fac, fun, ...))
-      else if (prod(sapply(private$.axes, function(x) x$length)) < CF.options$memory_cell_limit)
+      itemsize <- .nc_type_size(self$data_type) %||% 8L
+      if (prod(sapply(private$.axes, function(x) x$length)) * itemsize < CF.options$memory_cell_limit)
         # Read the whole data array because size is manageable
         return(.process.data(self$read_data(), tdim, fac, fun, ...))
 
@@ -203,21 +355,20 @@ CFVariable <- R6::R6Class("CFVariable",
       for (l in 1L:lvls) {
         indices <- which(ndx == l)
         dff <- diff(indices)
-        if (all(dff == 1L)) {       # Data is contiguous per factor level
+        if (all(dff == 1L)) {
           rng <- range(indices)
-          start[tdim] <- rng[1L]
-          count[tdim] <- rng[2L] - rng[1L] + 1L
-          values <- self$read_chunk(start, count)
-        } else {                    # Era factors have disparate indices
+          st <- start; st[tdim] <- rng[1L]
+          ct <- count; ct[tdim] <- rng[2L] - rng[1L] + 1L
+          runs <- list(list(start = st, count = ct))
+        } else {
           cutoffs <- c(0L, which(c(dff, 2L) > 1L))
-          values <- lapply(2L:length(cutoffs), function(i) {
-            start[tdim] <- indices[cutoffs[i - 1L] + 1L]
-            count[tdim] <- cutoffs[i] - cutoffs[i - 1L]
-            self$read_chunk(start, count)
+          runs <- lapply(2L:length(cutoffs), function(i) {
+            st <- start; st[tdim] <- indices[cutoffs[i - 1L] + 1L]
+            ct <- count; ct[tdim] <- cutoffs[i] - cutoffs[i - 1L]
+            list(start = st, count = ct)
           })
-          values <- abind::abind(values, along = num_dims)
         }
-        d[[l]] <- .process.data(values, tdim, FUN = fun, ...)
+        d[[l]] <- private$process_level(runs, tdim, num_dims, fun, ...)
         # d is a list with lvls elements, each element a list with elements for
         # the number of function results, possibly 1; each element having an
         # array of dimensions from private$values that are not tdim.
@@ -599,12 +750,10 @@ CFVariable <- R6::R6Class("CFVariable",
       # Get the data for the result CFVariable
       d <- NULL
       if (is.null(aux)) {
-        # Regular axes selected so stay virtual if data has not been loaded yet
         if (!is.null(private$.values))
-          d <- self$read_chunk(start, count)
+          d <- self$read_window(start, count)
       } else {
-        # Auxiliary grids selected, index the data
-        d <- self$read_chunk(start, count)
+        d <- self$read_window(start, count)
 
         lon_idx <- which(sapply(out_axes, inherits, "CFAxisLongitude"))
         lat_idx <- which(sapply(out_axes, inherits, "CFAxisLatitude"))
@@ -640,9 +789,9 @@ CFVariable <- R6::R6Class("CFVariable",
       if (is.null(aux)) {
         v <- if (self$has_resource) {
           NCdims <- length(private$.NC_map$start)
+          nc <- private$to_nc_indices(start[seq_len(NCdims)], count[seq_len(NCdims)])
           CFVariable$new(private$.NCobj, group = grp, values = d, axes = out_axes,
-                         start = start[1:NCdims] + private$.NC_map$start - 1L, count = count[1:NCdims],
-                         attributes = atts)
+                         start = nc$start, count = nc$count, attributes = atts)
         } else
           CFVariable$new(self$name, group = grp, values = d, axes = out_axes, attributes = atts)
         v$crs <- private$.crs
@@ -1555,7 +1704,7 @@ dimnames.CFVariable <- function(x) {
       }
     }
   }
-  data <- x$read_chunk(start, count)
+  data <- x$read_window(start, count)
 
   # Apply dimension data and other attributes
   if (length(x$axes) && length(dim(data)) == length(dnames)) { # dimensions may have been dropped automatically, e.g. NC_CHAR to character string
diff --git a/R/utils.R b/R/utils.R
index 425a01f..d58b3a3 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -27,6 +27,104 @@ netcdf_data_types <- c("NC_BYTE", "NC_UBYTE", "NC_CHAR", "NC_SHORT",
          "NC_NAT")
 }
 
+# Approximate byte size of a single element of the given netCDF data type.
+# Returns NA for an unrecognized type (the caller should then skip any
+# byte-based guard rather than block on an unknown quantity).
+.nc_type_size <- function(nc_type) {
+  switch(nc_type,
+         NC_BYTE = , NC_UBYTE = , NC_CHAR = 1L,
+         NC_SHORT = , NC_USHORT = 2L,
+         NC_INT = , NC_UINT = , NC_FLOAT = 4L,
+         NC_INT64 = , NC_UINT64 = , NC_DOUBLE = 8L,
+         NC_STRING = 8L,  # a floor: actual size is per-string and unknowable in advance
+         NA_integer_)
+}
+
+# Human-readable byte count for warning messages.
+.format_bytes <- function(bytes) {
+  units <- c("B", "KB", "MB", "GB", "TB")
+  e <- min(length(units) - 1L, max(0L, floor(log(max(bytes, 1), 1024))))
+  sprintf("%.1f %s", bytes / 1024^e, units[e + 1L])
+}
+
+# Plan a sequence of (start, count) sub-reads tiling the requested window,
+# each within budget_bytes, aligned to the native chunk grid where that
+# alignment is actually meaningful.
+#
+# shape, chunk_shape: integer vectors, same order as start/count.
+#   chunk_shape may be NULL for contiguous storage.
+# start, count: the requested window.
+# itemsize: bytes per element.
+# budget_bytes: cap on any single returned block, in bytes.
+#
+# Returns a list of `list(start = ..., count = ...)`.
+.chunk_plan <- function(shape, chunk_shape, start, count, itemsize, budget_bytes) {
+  if (itemsize > budget_bytes)
+    stop("Argument `budget_bytes` is smaller than a single element; cannot plan a read", call. = FALSE)
+
+  nd <- length(shape)
+  if (is.null(chunk_shape)) {
+    chunk_shape <- count
+    chunk_shape[1L] <- 1L  # contiguous storage: slab along the first dimension
+  }
+  chunk_shape <- as.integer(pmin(chunk_shape, shape))
+
+  last_idx     <- start + count - 1L
+  first_native <- (start - 1L) %/% chunk_shape
+  last_native  <- (last_idx - 1L) %/% chunk_shape
+  n_native     <- last_native - first_native + 1L
+  native_bytes <- prod(as.numeric(chunk_shape)) * itemsize
+
+  if (native_bytes <= budget_bytes) {
+    # One native chunk fits: try to GROUP several whole native chunks into
+    # one request, up to budget, growing whichever dimension has the most
+    # native chunks left to absorb.
+    group <- rep(1L, nd)
+    repeat {
+      grown <- FALSE
+      for (d in order(-(n_native / group))) {
+        if (group[d] >= n_native[d]) next
+        trial <- group; trial[d] <- group[d] + 1L
+        extent <- pmin(trial * chunk_shape, count)
+        if (prod(as.numeric(extent)) * itemsize > budget_bytes) next
+        group <- trial; grown <- TRUE
+        break
+      }
+      if (!grown) break
+    }
+    tile        <- pmin(group * chunk_shape, count)
+    grid_origin <- first_native * chunk_shape + 1L
+  } else {
+    # A single native chunk already exceeds budget: SHRINK below it, reducing
+    # the largest dimension(s) first. Every request that touches a given native
+    # chunk under this regime causes that chunk to be decompressed again
+    # internally -- unavoidable once the chunk itself doesn't fit the budget,
+    # not a defect of the tiling. This finds a tile that fits, via simple
+    # halving; not necessarily the largest one that would still fit.
+    tile <- chunk_shape
+    ord  <- order(-chunk_shape)
+    i <- 1L
+    while (prod(as.numeric(tile)) * itemsize > budget_bytes && i <= nd) {
+      d <- ord[i]
+      while (tile[d] > 1L && prod(as.numeric(tile)) * itemsize > budget_bytes)
+        tile[d] <- max(1L, tile[d] %/% 2L)
+      i <- i + 1L
+    }
+    tile        <- pmin(tile, count)
+    grid_origin <- start  # no alignment benefit within a single native chunk
+  }
+
+  seqs   <- lapply(seq_len(nd), function(d) seq(grid_origin[d], last_idx[d], by = tile[d]))
+  combos <- expand.grid(seqs, KEEP.OUT.ATTRS = FALSE)
+
+  lapply(seq_len(nrow(combos)), function(i) {
+    s_grid <- as.integer(combos[i, ])
+    e      <- pmin(s_grid + tile - 1L, last_idx)
+    s      <- pmax(s_grid, start)
+    list(start = s, count = as.integer(e - s + 1L))
+  })
+}
+
 # This function is a bare-bones implementation of `apply(X, MARGIN, tapply, INDEX, FUN, ...)`,
 # i.e. apply a factor over a dimension of an array. There are several restrictions
 # compared to the base::apply/tapply pair (but note that function arguments are
diff --git a/R/zzz.R b/R/zzz.R
index 19fa22a..8d8ccdb 100644
--- a/R/zzz.R
+++ b/R/zzz.R
@@ -20,7 +20,7 @@ CF.options <- new.env(parent = emptyenv())
   assign("newVarId", function() {CF$currentVarId <- CF$currentVarId - 1L; CF$currentVarId}, envir = CF)
 
   # User-modifiable options
-  assign("memory_cell_limit", 1e8, envir = CF.options)
+  assign("memory_cell_limit", 1e9, envir = CF.options)
   assign("digits", 6L, envir = CF.options)
   assign("cache_stale_days", 90, envir = CF.options)
 }
diff --git a/tests/testthat/test-CFobjects.R b/tests/testthat/test-CFobjects.R
index f84b2e8..3a37862 100644
--- a/tests/testthat/test-CFobjects.R
+++ b/tests/testthat/test-CFobjects.R
@@ -169,3 +169,55 @@ test_that("Math and Ops functions", {
   dvabove0 <- dvcos > 0
   expect_true(is.logical(dvabove0$raw()))
 })
+
+test_that("subset() on a resource-backed variable and its axes computes correct nc-space indices", {
+  skip_if_not_installed("RNetCDF")
+
+  fn <- tempfile(fileext = ".nc")
+  nc <- RNetCDF::create.nc(fn, format = "netcdf4")
+  RNetCDF::dim.def.nc(nc, "lon", 10L)
+  RNetCDF::dim.def.nc(nc, "lat", 8L)
+  RNetCDF::dim.def.nc(nc, "time", 6L)
+
+  RNetCDF::var.def.nc(nc, "lon", "NC_DOUBLE", "lon")
+  RNetCDF::att.put.nc(nc, "lon", "units", "NC_CHAR", "degrees_east")
+  RNetCDF::att.put.nc(nc, "lon", "standard_name", "NC_CHAR", "longitude")
+  RNetCDF::att.put.nc(nc, "lon", "axis", "NC_CHAR", "X")
+  RNetCDF::var.put.nc(nc, "lon", seq(0, 9))
+
+  RNetCDF::var.def.nc(nc, "lat", "NC_DOUBLE", "lat")
+  RNetCDF::att.put.nc(nc, "lat", "units", "NC_CHAR", "degrees_north")
+  RNetCDF::att.put.nc(nc, "lat", "standard_name", "NC_CHAR", "latitude")
+  RNetCDF::att.put.nc(nc, "lat", "axis", "NC_CHAR", "Y")
+  RNetCDF::var.put.nc(nc, "lat", seq(0, 7))
+
+  RNetCDF::var.def.nc(nc, "time", "NC_DOUBLE", "time")
+  RNetCDF::att.put.nc(nc, "time", "units", "NC_CHAR", "days since 2000-01-01")
+  RNetCDF::att.put.nc(nc, "time", "standard_name", "NC_CHAR", "time")
+  RNetCDF::att.put.nc(nc, "time", "axis", "NC_CHAR", "T")
+  RNetCDF::att.put.nc(nc, "time", "calendar", "NC_CHAR", "standard")
+  RNetCDF::var.put.nc(nc, "time", seq(0, 5))
+
+  RNetCDF::var.def.nc(nc, "temp", "NC_DOUBLE", c("lon", "lat", "time"),
+                      chunking = TRUE, chunksizes = c(4L, 3L, 2L))
+  data <- array(as.double(seq_len(10L * 8L * 6L)), dim = c(10L, 8L, 6L))
+  RNetCDF::var.put.nc(nc, "temp", data)
+  RNetCDF::close.nc(nc)
+  on.exit(unlink(fn))
+
+  ds <- open_ncdf(fn)
+  v <- ds[[ds$var_names[1L]]]
+
+  sub <- v$subset(X = c(3, 5), Y = c(2, 4))  # coordinate values, per subset()'s contract
+
+  is_x <- vapply(sub$axes, function(a) a$orientation, character(1L)) == "X"
+  is_y <- vapply(sub$axes, function(a) a$orientation, character(1L)) == "Y"
+  lon_idx <- match(sub$axes[[which(is_x)]]$values, seq(0, 9))
+  lat_idx <- match(sub$axes[[which(is_y)]]$values, seq(0, 7))
+
+  # Self-consistent check: whatever indices the subsetted axes themselves
+  # report having selected, the variable's data must match exactly those
+  # indices in the original array -- doesn't depend on hand-computing the
+  # expected index range separately.
+  expect_equal(sub$values, data[lon_idx, lat_idx, ])
+})
diff --git a/tests/testthat/test-chunking.R b/tests/testthat/test-chunking.R
new file mode 100644
index 0000000..9ef1386
--- /dev/null
+++ b/tests/testthat/test-chunking.R
@@ -0,0 +1,266 @@
+test_that(".chunk_plan: whole array fits in budget as a single block", {
+  plan <- .chunk_plan(shape = 100L, chunk_shape = 10L, start = 1L, count = 100L,
+                      itemsize = 4L, budget_bytes = 1000L)
+  expect_length(plan, 1L)
+  expect_equal(plan[[1L]]$start, 1L)
+  expect_equal(plan[[1L]]$count, 100L)
+})
+
+test_that(".chunk_plan: grouping, full window, tight budget", {
+  # native chunk = 10 elements = 40 bytes; budget covers 2 native chunks (80B)
+  plan <- .chunk_plan(shape = 100L, chunk_shape = 10L, start = 1L, count = 100L,
+                      itemsize = 4L, budget_bytes = 80L)
+  counts <- sapply(plan, function(p) p$count)
+  starts <- sapply(plan, function(p) p$start)
+  expect_true(all(counts * 4L <= 80L))
+  expect_equal(sum(counts), 100L)
+  # contiguous, gap-free, non-overlapping coverage
+  ord <- order(starts)
+  ends <- starts[ord] + counts[ord] - 1L
+  expect_equal(starts[ord][1L], 1L)
+  expect_equal(ends[length(ends)], 100L)
+  expect_true(all(starts[ord][-1L] == ends[-length(ends)] + 1L))
+})
+
+test_that(".chunk_plan: grouping, off-boundary partial window", {
+  # window [15,64] within a 100-length, chunk-10 array; budget = 2 native chunks
+  plan <- .chunk_plan(shape = 100L, chunk_shape = 10L, start = 15L, count = 50L,
+                      itemsize = 4L, budget_bytes = 80L)
+  counts <- sapply(plan, function(p) p$count)
+  starts <- sapply(plan, function(p) p$start)
+  expect_true(all(counts * 4L <= 80L))
+  expect_equal(sum(counts), 50L)  # exact coverage, no double-counted overlap
+  ord <- order(starts)
+  ends <- starts[ord] + counts[ord] - 1L
+  expect_equal(starts[ord][1L], 15L)
+  expect_equal(ends[length(ends)], 64L)
+  expect_true(all(starts[ord][-1L] == ends[-length(ends)] + 1L))
+})
+
+test_that(".chunk_plan: multi-dimensional grouping covers the window exactly once", {
+  shape <- c(288L, 180L, 9645L)        # (lon, lat, time), as RNetCDF reports it
+  chunk <- c(58L, 36L, 1929L)          # this file's actual native chunk
+  start <- c(1L, 1L, 1L)
+  count <- c(288L, 180L, 1929L)        # one full time-chunk's worth, full space
+  itemsize <- 4L
+  budget <- 8L * 1024^2                # 8MB -- forces splitting below the ~16MB full slab
+
+  plan <- .chunk_plan(shape, chunk, start, count, itemsize, budget)
+  total_elems <- sum(sapply(plan, function(p) prod(p$count)))
+  expect_equal(total_elems, prod(count))
+  expect_true(all(sapply(plan, function(p) prod(p$count) * itemsize <= budget)))
+})
+
+test_that(".chunk_plan: single native chunk exceeds budget, must shrink below it", {
+  shape <- c(288L, 180L, 9645L)
+  chunk <- c(58L, 36L, 1929L)          # ~16MB uncompressed
+  plan <- .chunk_plan(shape, chunk, start = c(1L, 1L, 1L), count = c(58L, 36L, 1929L),
+                      itemsize = 4L, budget_bytes = 100L)  # absurdly tight on purpose
+  expect_true(all(sapply(plan, function(p) prod(p$count) * 4L <= 100L)))
+  total_elems <- sum(sapply(plan, function(p) prod(p$count)))
+  expect_equal(total_elems, 58L * 36L * 1929L)
+})
+
+test_that(".chunk_plan: contiguous storage falls back to slabbing on the first dimension", {
+  plan <- .chunk_plan(shape = c(50L, 20L), chunk_shape = NULL,
+                      start = c(1L, 1L), count = c(50L, 20L),
+                      itemsize = 8L, budget_bytes = 8L * 20L * 5L)  # 5 rows at a time
+  counts <- t(sapply(plan, function(p) p$count))
+  expect_true(all(counts[, 2L] == 20L))       # second dim never split
+  expect_equal(sum(counts[, 1L]), 50L)
+})
+
+test_that(".chunk_plan: errors when budget can't hold a single element", {
+  expect_error(.chunk_plan(10L, 5L, 1L, 10L, itemsize = 8L, budget_bytes = 4L))
+})
+
+test_that(".chunk_plan: contiguous storage shrinks non-first dimensions when they dominate", {
+  # First dimension is small; the fallback's chunk_shape[1] <- 1L must not be
+  # treated as a floor -- the shrink logic still has to reduce dims 2/3.
+  shape <- c(5L, 20000L, 20000L)
+  plan <- .chunk_plan(shape, chunk_shape = NULL,
+                      start = c(1L, 1L, 1L), count = shape,
+                      itemsize = 8L, budget_bytes = 1e6)
+  expect_true(all(sapply(plan, function(p) prod(p$count) * 8L <= 1e6)))
+  total <- sum(sapply(plan, function(p) prod(p$count)))
+  expect_equal(total, prod(shape))
+})
+
+test_that("read_window splits large reads via .chunk_plan and reassembles correctly", {
+  skip_if_not_installed("RNetCDF")
+
+  fn <- tempfile(fileext = ".nc")
+  nc <- RNetCDF::create.nc(fn, format = "netcdf4")
+  RNetCDF::dim.def.nc(nc, "x", 20L)
+  RNetCDF::dim.def.nc(nc, "y", 15L)
+  RNetCDF::var.def.nc(nc, "v", "NC_DOUBLE", c("x", "y"),
+                      chunking = TRUE, chunksizes = c(4L, 3L))
+  data <- matrix(as.double(seq_len(20L * 15L)), nrow = 20L, ncol = 15L)
+  RNetCDF::var.put.nc(nc, "v", data)
+  RNetCDF::close.nc(nc)
+  on.exit(unlink(fn))
+
+  old_limit <- CF.options$memory_cell_limit
+  assign("memory_cell_limit", 200L * 8L, envir = CF.options)  # forces splitting
+  on.exit(assign("memory_cell_limit", old_limit, envir = CF.options), add = TRUE)
+
+  ds <- open_ncdf(fn)
+  v <- ds[[ds$var_names[1L]]]
+
+  expect_warning(full <- v$read_window(c(1L, 1L), c(20L, 15L)))
+  slice <- v$read_window(c(3L, 2L), c(10L, 8L))
+
+  expect_equal(full, data)
+  expect_equal(slice, data[3:12, 2:9])
+})
+
+test_that("process_data() streams sum/mean via .chunk_plan() and matches the materialized result", {
+  skip_if_not_installed("RNetCDF")
+
+  fn <- tempfile(fileext = ".nc")
+  nc <- RNetCDF::create.nc(fn, format = "netcdf4")
+  RNetCDF::dim.def.nc(nc, "lon", 6L)
+  RNetCDF::dim.def.nc(nc, "lat", 5L)
+  RNetCDF::dim.def.nc(nc, "time", 60L)
+
+  RNetCDF::var.def.nc(nc, "lon", "NC_DOUBLE", "lon")
+  RNetCDF::att.put.nc(nc, "lon", "units", "NC_CHAR", "degrees_east")
+  RNetCDF::att.put.nc(nc, "lon", "standard_name", "NC_CHAR", "longitude")
+  RNetCDF::att.put.nc(nc, "lon", "axis", "NC_CHAR", "X")
+  RNetCDF::var.put.nc(nc, "lon", seq(0, 5))
+
+  RNetCDF::var.def.nc(nc, "lat", "NC_DOUBLE", "lat")
+  RNetCDF::att.put.nc(nc, "lat", "units", "NC_CHAR", "degrees_north")
+  RNetCDF::att.put.nc(nc, "lat", "standard_name", "NC_CHAR", "latitude")
+  RNetCDF::att.put.nc(nc, "lat", "axis", "NC_CHAR", "Y")
+  RNetCDF::var.put.nc(nc, "lat", seq(0, 4))
+
+  RNetCDF::var.def.nc(nc, "time", "NC_DOUBLE", "time")
+  RNetCDF::att.put.nc(nc, "time", "units", "NC_CHAR", "days since 2000-01-01")
+  RNetCDF::att.put.nc(nc, "time", "standard_name", "NC_CHAR", "time")
+  RNetCDF::att.put.nc(nc, "time", "axis", "NC_CHAR", "T")
+  RNetCDF::att.put.nc(nc, "time", "calendar", "NC_CHAR", "standard")
+  RNetCDF::var.put.nc(nc, "time", 0:59)
+
+  RNetCDF::var.def.nc(nc, "temp", "NC_DOUBLE", c("lon", "lat", "time"),
+                      chunking = TRUE, chunksizes = c(2L, 2L, 5L))
+  set.seed(1)
+  data <- array(rnorm(6L * 5L * 60L), dim = c(6L, 5L, 60L))
+  data[1L, 1L, 3L] <- NA_real_  # exercises na.rm handling at one location
+  RNetCDF::var.put.nc(nc, "temp", data)
+  RNetCDF::close.nc(nc)
+  on.exit(unlink(fn))
+
+  old_limit <- CF.options$memory_cell_limit
+  on.exit(assign("memory_cell_limit", old_limit, envir = CF.options), add = TRUE)
+
+  ds <- open_ncdf(fn)
+  v <- ds[[ds$var_names[1L]]]
+
+  assign("memory_cell_limit", 1e9, envir = CF.options)  # generous: materialize path
+  sum_full  <- v$summarise("total", sum,  "month", na.rm = TRUE)
+  mean_full <- v$summarise("avg",   mean, "month", na.rm = TRUE)
+
+  assign("memory_cell_limit", 6L * 5L * 3L * 8L, envir = CF.options)  # forces streaming
+  sum_stream  <- v$summarise("total", sum,  "month", na.rm = TRUE)
+  mean_stream <- v$summarise("avg",   mean, "month", na.rm = TRUE)
+
+  expect_equal(sum_stream$values, sum_full$values)
+  expect_equal(mean_stream$values, mean_full$values)
+
+  # na.rm = FALSE: the NA must propagate identically in both paths
+  assign("memory_cell_limit", 1e9, envir = CF.options)
+  sum_full_narm_false <- v$summarise("total", sum, "month")
+  assign("memory_cell_limit", 6L * 5L * 3L * 8L, envir = CF.options)
+  sum_stream_narm_false <- v$summarise("total", sum, "month")
+  expect_equal(sum_stream_narm_false$values, sum_full_narm_false$values)
+  expect_true(is.na(sum_stream_narm_false$values[1L, 1L, 1L]))
+})
+
+test_that("process_data() streams min/max/range, replicating base R's Inf/warning for all-NA locations", {
+  skip_if_not_installed("RNetCDF")
+
+  lon  <- seq(0, 5)
+  lat  <- seq(0, 4)
+  time <- as.character(as.Date("2000-01-01") + 0:59)
+  set.seed(2)
+  data <- array(rnorm(6L * 5L * 60L), dim = c(6L, 5L, 60L))
+  data[1L, 1L, ] <- NA_real_  # this location is entirely NA, every level
+  dimnames(data) <- list(as.character(lon), as.character(lat), time)
+  names(dimnames(data)) <- c("lon", "lat", "time")
+
+  # Ground truth: built entirely in memory via as_CF(), no file I/O, no
+  # chunking question -- process_data() takes the already-materialized
+  # branch immediately since .values is already set.
+  v_mem <- as_CF("temp", data)
+
+  # The streamed path needs an actual on-disk, explicitly-chunked variable;
+  # as_CF()$save() can't set chunksizes (NCVariable's own header comment:
+  # "chunking - .ncdf4 - Not used"), so this part still needs raw RNetCDF.
+  fn <- tempfile(fileext = ".nc")
+  nc <- RNetCDF::create.nc(fn, format = "netcdf4")
+  RNetCDF::dim.def.nc(nc, "lon", 6L)
+  RNetCDF::dim.def.nc(nc, "lat", 5L)
+  RNetCDF::dim.def.nc(nc, "time", 60L)
+
+  RNetCDF::var.def.nc(nc, "lon", "NC_DOUBLE", "lon")
+  RNetCDF::att.put.nc(nc, "lon", "units", "NC_CHAR", "degrees_east")
+  RNetCDF::att.put.nc(nc, "lon", "standard_name", "NC_CHAR", "longitude")
+  RNetCDF::att.put.nc(nc, "lon", "axis", "NC_CHAR", "X")
+  RNetCDF::var.put.nc(nc, "lon", lon)
+
+  RNetCDF::var.def.nc(nc, "lat", "NC_DOUBLE", "lat")
+  RNetCDF::att.put.nc(nc, "lat", "units", "NC_CHAR", "degrees_north")
+  RNetCDF::att.put.nc(nc, "lat", "standard_name", "NC_CHAR", "latitude")
+  RNetCDF::att.put.nc(nc, "lat", "axis", "NC_CHAR", "Y")
+  RNetCDF::var.put.nc(nc, "lat", lat)
+
+  RNetCDF::var.def.nc(nc, "time", "NC_DOUBLE", "time")
+  RNetCDF::att.put.nc(nc, "time", "units", "NC_CHAR", "days since 2000-01-01")
+  RNetCDF::att.put.nc(nc, "time", "standard_name", "NC_CHAR", "time")
+  RNetCDF::att.put.nc(nc, "time", "axis", "NC_CHAR", "T")
+  RNetCDF::att.put.nc(nc, "time", "calendar", "NC_CHAR", "standard")
+  RNetCDF::var.put.nc(nc, "time", 0:59)
+
+  RNetCDF::var.def.nc(nc, "temp", "NC_DOUBLE", c("lon", "lat", "time"),
+                      chunking = TRUE, chunksizes = c(2L, 2L, 5L))
+  RNetCDF::var.put.nc(nc, "temp", data)  # same poisoned array, written verbatim
+  RNetCDF::close.nc(nc)
+  on.exit(unlink(fn))
+
+  old_limit <- CF.options$memory_cell_limit
+  on.exit(assign("memory_cell_limit", old_limit, envir = CF.options), add = TRUE)
+  # Feasibility floor for this shape: 6*5 (full space) * 5 (native time-chunk)
+  # * 8 bytes = 1200B. Must exceed that for streaming to engage at all.
+  assign("memory_cell_limit", 6L * 5L * 10L * 8L, envir = CF.options)
+
+  ds <- open_ncdf(fn)
+  v  <- ds[[ds$var_names[1L]]]
+
+  min_full   <- suppressWarnings(v_mem$summarise("mn", min, "month", na.rm = TRUE))
+  max_full   <- suppressWarnings(v_mem$summarise("mx", max, "month", na.rm = TRUE))
+  range_full <- suppressWarnings(v_mem$summarise(c("rmin", "rmax"), range, "month", na.rm = TRUE))
+
+  min_stream   <- suppressWarnings(v$summarise("mn", min, "month", na.rm = TRUE))
+  max_stream   <- suppressWarnings(v$summarise("mx", max, "month", na.rm = TRUE))
+  range_stream <- suppressWarnings(v$summarise(c("rmin", "rmax"), range, "month", na.rm = TRUE))
+
+  expect_equal(min_stream$values, min_full$values)
+  expect_equal(max_stream$values, max_full$values)
+  expect_equal(range_stream$rmin$values, range_full$rmin$values)
+  expect_equal(range_stream$rmax$values, range_full$rmax$values)
+
+  # Output shape is (month, lon, lat) -- the new time-like axis first, per
+  # process_data()'s aperm(x, c(num_dims, 1:(num_dims-1))) and summarise()'s
+  # axes <- c(new_ax, other_axes). Location (lon=1, lat=1) is [, 1L, 1L].
+  expect_true(all(is.infinite(min_stream$values[, 1L, 1L]) & min_stream$values[, 1L, 1L] > 0))
+  expect_true(all(is.infinite(max_stream$values[, 1L, 1L]) & max_stream$values[, 1L, 1L] < 0))
+  expect_false(any(is.infinite(min_stream$values[, -1L, ])))
+
+  ws <- character(0L)
+  withCallingHandlers(
+    v$summarise("mn2", min, "month", na.rm = TRUE),
+    warning = function(w) { ws <<- c(ws, conditionMessage(w)); invokeRestart("muffleWarning") }
+  )
+  expect_true(any(grepl("no non-missing arguments to min", ws, fixed = TRUE)))
+})