Repository: dcs4cop/xcube
Pull request #260 (2020-01-15 14:29:13)

Spatial subset 
Integrated new rectification in default input processor. It is now the default reprojection method in  `xcube.core.gen.iproc.XYInputProcessor`, if ground controll points are not specified, i.e. the input processor  is configured with `xy_gcp_step` being None.
Introduced new function `xcube.core.select.select_spatial_subset`.

Closes issue #206 .


Patch:
diff --git a/CHANGES.md b/CHANGES.md
index 32bdb8f9..23c4839c 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -2,7 +2,10 @@
 
 ### New
 
-
+* Integrated new rectification in default input processor. It is now the default reprojection method in 
+  `xcube.core.gen.iproc.XYInputProcessor`, if ground control points are not specified, i.e. the input processor 
+  is configured with `xy_gcp_step` being None.
+* Introduced new function `xcube.core.select.select_spatial_subset`.
 * Introduced new (ortho-)rectification algorithm allowing reprojection of 
   satellite images that come with (terrain-corrected) geo-locations for every pixel.
 
diff --git a/xcube/core/evaluate.py b/xcube/core/evaluate.py
index 8d3c2ca5..ebfb2e4c 100644
--- a/xcube/core/evaluate.py
+++ b/xcube/core/evaluate.py
@@ -89,6 +89,8 @@ def evaluate_dataset(dataset: xr.Dataset,
             if var_props is None:
                 var_props = dict()
 
+        do_load = var_props.get('load', False)
+
         expression = var_props.get('expression')
         if expression:
             # Compute new variable
@@ -100,6 +102,8 @@ def evaluate_dataset(dataset: xr.Dataset,
                 if hasattr(computed_array, 'attrs'):
                     var = computed_array
                     var.attrs.update(var_props)
+                if do_load:
+                    computed_array.load()
                 namespace[var_name] = computed_array
 
         valid_pixel_expression = var_props.get('valid_pixel_expression')
@@ -115,6 +119,8 @@ def evaluate_dataset(dataset: xr.Dataset,
                 masked_var = var.where(valid_mask)
                 if hasattr(masked_var, 'attrs'):
                     masked_var.attrs.update(var_props)
+                if do_load:
+                    masked_var.load()
                 namespace[var_name] = masked_var
 
     computed_dataset = dataset.copy()
diff --git a/xcube/core/gen/gen.py b/xcube/core/gen/gen.py
index 7944a6f8..c44331e1 100644
--- a/xcube/core/gen/gen.py
+++ b/xcube/core/gen/gen.py
@@ -30,12 +30,13 @@ import warnings
 from typing import Any, Callable, Dict, Sequence, Tuple
 
 import xarray as xr
-
 from xcube.core.dsio import DatasetIO, find_dataset_io, guess_dataset_format, rimraf
 from xcube.core.evaluate import evaluate_dataset
 from xcube.core.gen.defaults import DEFAULT_OUTPUT_PATH, DEFAULT_OUTPUT_RESAMPLING, DEFAULT_OUTPUT_SIZE
 from xcube.core.gen.iproc import InputProcessor, find_input_processor_class
-from xcube.core.select import select_vars
+from xcube.core.geocoding import GeoCoding
+from xcube.core.imgeom import ImageGeom
+from xcube.core.select import select_spatial_subset, select_variables_subset
 from xcube.core.timecoord import add_time_coords, from_time_in_days_since_1970
 from xcube.core.timeslice import find_time_slice
 from xcube.core.update import update_dataset_attrs, update_dataset_temporal_attrs, update_dataset_var_attrs
@@ -216,6 +217,11 @@ def _process_input(input_processor: InputProcessor,
     time_index, update_mode = find_time_slice(output_path,
                                               from_time_in_days_since_1970((time_range[0] + time_range[1]) / 2))
 
+    width, height = output_size
+    x_min, y_min, x_max, y_max = output_region
+    xy_res = max((x_max - x_min) / width, (y_max - y_min) / height)
+    output_geom = ImageGeom(size=output_size, x_min=x_min, y_min=y_min, xy_res=xy_res, is_geo_crs=True)
+
     steps = []
 
     # noinspection PyShadowingNames
@@ -224,6 +230,25 @@ def _process_input(input_processor: InputProcessor,
 
     steps.append((step1, 'pre-processing input slice'))
 
+    geo_coding = None
+
+    # noinspection PyShadowingNames
+    def step1a(input_slice):
+        nonlocal geo_coding
+        geo_coding = GeoCoding.from_dataset(input_slice)
+        subset = select_spatial_subset(input_slice,
+                                       xy_bbox=output_geom.xy_bbox,
+                                       xy_border=output_geom.xy_res,
+                                       ij_border=1,
+                                       geo_coding=geo_coding)
+        if subset is None:
+            monitor('no spatial overlap with input')
+        elif subset is not input_slice:
+            geo_coding = GeoCoding.from_dataset(subset)
+        return subset
+
+    steps.append((step1a, 'spatial subsetting'))
+
     # noinspection PyShadowingNames
     def step2(input_slice):
         return evaluate_dataset(input_slice, processed_variables=processed_variables)
@@ -235,16 +260,17 @@ def _process_input(input_processor: InputProcessor,
         extra_vars = input_processor.get_extra_vars(input_slice)
         selected_variables = set([var_name for var_name, _ in output_variables])
         selected_variables.update(extra_vars or set())
-        return select_vars(input_slice, selected_variables)
+        return select_variables_subset(input_slice, selected_variables)
 
     steps.append((step3, 'selecting input slice variables'))
 
     # noinspection PyShadowingNames
     def step4(input_slice):
+        # noinspection PyTypeChecker
         return input_processor.process(input_slice,
-                                       dst_size=output_size,
-                                       dst_region=output_region,
-                                       dst_resampling=output_resampling,
+                                       geo_coding=geo_coding,
+                                       output_geom=output_geom,
+                                       output_resampling=output_resampling,
                                        include_non_spatial_vars=False)
 
     steps.append((step4, 'transforming input slice'))
diff --git a/xcube/core/gen/iproc.py b/xcube/core/gen/iproc.py
index 90c39832..977b8ea2 100644
--- a/xcube/core/gen/iproc.py
+++ b/xcube/core/gen/iproc.py
@@ -21,11 +21,10 @@
 
 import warnings
 from abc import ABCMeta, abstractmethod
-from typing import Collection, Optional, Tuple, Union, Any, Mapping, Dict
+from typing import Any, Collection, Dict, Mapping, Optional, Tuple, Union
 
 import numpy as np
 import xarray as xr
-
 from xcube.constants import CRS_WKT_EPSG_4326
 from xcube.constants import EXTENSION_POINT_INPUT_PROCESSORS
 from xcube.core.geocoding import GeoCoding
@@ -208,9 +207,9 @@ class InputProcessor(ExtensionComponent, metaclass=ABCMeta):
     @abstractmethod
     def process(self,
                 dataset: xr.Dataset,
-                dst_size: Tuple[int, int],
-                dst_region: Tuple[float, float, float, float],
-                dst_resampling: str,
+                geo_coding: GeoCoding,
+                output_geom: ImageGeom,
+                output_resampling: str,
                 include_non_spatial_vars=False) -> xr.Dataset:
         """
         Perform spatial transformation into the cube's WGS84 SRS such that all variables in the output dataset
@@ -219,9 +218,9 @@ class InputProcessor(ExtensionComponent, metaclass=ABCMeta):
         * must have *dst_region* as their bounding box in geographic coordinates.
 
         :param dataset: The input dataset.
-        :param dst_size: The output size in pixels as tuple (width ,height).
-        :param dst_region: The output region in coordinates of the target CRS. A tuple (x_min, y_min, x_max, y_max).
-        :param dst_resampling: The spatial resampling method.
+        :param geo_coding: The input's geo-coding.
+        :param output_geom: The output's spatial image geometry.
+        :param output_resampling: The output's spatial resampling method.
         :param include_non_spatial_vars: Whether to include non-spatial variables in the output.
         :return: The transformed output dataset or the original one, if no transformation is required.
         """
@@ -299,9 +298,9 @@ class XYInputProcessor(InputProcessor, metaclass=ABCMeta):
 
     def process(self,
                 dataset: xr.Dataset,
-                dst_size: Tuple[int, int],
-                dst_region: Tuple[float, float, float, float],
-                dst_resampling: str,
+                geo_coding: GeoCoding,
+                output_geom: ImageGeom,
+                output_resampling: str,
                 include_non_spatial_vars=False) -> xr.Dataset:
         """
         Perform reprojection using tie-points / ground control points.
@@ -314,26 +313,17 @@ class XYInputProcessor(InputProcessor, metaclass=ABCMeta):
             if reprojection_info.xy_tp_gcp_step is not None:
                 warnings.warn(f'{warn_prefix}: ignoring '
                               f'reprojection_info.xy_tp_gcp_step = {reprojection_info.xy_tp_gcp_step!r}')
-            if dst_resampling != 'Nearest':
+            if output_resampling != 'Nearest':
                 warnings.warn(f'{warn_prefix}: ignoring '
-                              f'dst_resampling = {dst_resampling!r}')
+                              f'dst_resampling = {output_resampling!r}')
             if include_non_spatial_vars:
                 warnings.warn(f'{warn_prefix}: ignoring '
                               f'include_non_spatial_vars = {include_non_spatial_vars!r}')
 
-            geo_coding = GeoCoding.from_dataset(dataset, xy_names=reprojection_info.xy_names)
-
-            dst_width, dst_height = dst_size
-            dst_x1, dst_y1, dst_x2, dst_y2 = dst_region
-            dst_res = max((dst_x2 - dst_x1) / dst_width, (dst_y2 - dst_y1) / dst_height)
-
-            output_geom = ImageGeom(size=(dst_width, dst_height),
-                                    x_min=dst_x1,
-                                    y_min=dst_y1,
-                                    xy_res=dst_res,
-                                    is_geo_crs=geo_coding.is_geo_crs)
+            geo_coding = geo_coding.derive(x_name=reprojection_info.xy_names[0], y_name=reprojection_info.xy_names[1])
 
             dataset = rectify_dataset(dataset,
+                                      compute_subset=False,
                                       geo_coding=geo_coding,
                                       output_geom=output_geom,
                                       is_y_axis_inverted=True)
@@ -350,9 +340,9 @@ class XYInputProcessor(InputProcessor, metaclass=ABCMeta):
                                          src_xy_crs=reprojection_info.xy_crs,
                                          src_xy_gcp_step=reprojection_info.xy_gcp_step or 1,
                                          src_xy_tp_gcp_step=reprojection_info.xy_tp_gcp_step or 1,
-                                         dst_size=dst_size,
-                                         dst_region=dst_region,
-                                         dst_resampling=dst_resampling,
+                                         dst_size=output_geom.size,
+                                         dst_region=output_geom.xy_bbox,
+                                         dst_resampling=output_resampling,
                                          include_non_spatial_vars=include_non_spatial_vars)
 
 
diff --git a/xcube/core/rectify.py b/xcube/core/rectify.py
index f7c28a98..c8aeb5dc 100644
--- a/xcube/core/rectify.py
+++ b/xcube/core/rectify.py
@@ -19,7 +19,7 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 # SOFTWARE.
 
-from typing import Sequence, Tuple, Optional, Union, Mapping
+from typing import Mapping, Optional, Sequence, Tuple, Union
 
 import dask.array as da
 import numba as nb
@@ -28,6 +28,7 @@ import xarray as xr
 
 from xcube.core.geocoding import GeoCoding
 from xcube.core.imgeom import ImageGeom
+from xcube.core.select import select_spatial_subset
 from xcube.util.dask import compute_array_from_func
 
 
@@ -40,6 +41,7 @@ def rectify_dataset(dataset: xr.Dataset,
                     tile_size: Union[int, Tuple[int, int]] = None,
                     compute_xy: bool = False,
                     compute_vars: bool = False,
+                    compute_subset: bool = True,
                     uv_delta: float = 1e-3) -> Optional[xr.Dataset]:
     """
     Reproject *dataset* using its per-pixel x,y coordinates or the given *geo_coding*.
@@ -80,11 +82,14 @@ def rectify_dataset(dataset: xr.Dataset,
 
     if output_geom is None:
         output_geom = ImageGeom.from_dataset(dataset, geo_coding=src_geo_coding)
-    else:
-        src_bbox = src_geo_coding.ij_bbox(output_geom.xy_bbox, ij_border=1, xy_border=output_geom.xy_res)
-        if src_bbox[0] == -1:
+    elif compute_subset:
+        dataset_subset = select_spatial_subset(dataset,
+                                               xy_bbox=output_geom.xy_bbox,
+                                               ij_border=1,
+                                               xy_border=output_geom.xy_res,
+                                               geo_coding=src_geo_coding)
+        if dataset_subset is None:
             return None
-        dataset_subset = select_spatial_subset(dataset, src_bbox, geo_coding=src_geo_coding)
         if dataset_subset is not dataset:
             src_geo_coding = GeoCoding.from_dataset(dataset_subset)
             src_x, src_y = src_geo_coding.x, src_geo_coding.y
@@ -271,30 +276,6 @@ def select_variables(dataset,
     return src_vars
 
 
-def select_spatial_subset(dataset: xr.Dataset,
-                          bbox: Tuple[int, int, int, int],
-                          geo_coding: GeoCoding = None,
-                          xy_names: Tuple[str, str] = None) -> Optional[xr.Dataset]:
-    """
-    Select a spatial subset of *dataset* for the bounding box *bbox*.
-
-    :param dataset: Source dataset.
-    :param bbox: Bounding box (i_min, i_min, j_max, j_max) in pixel coordinates.
-    :param geo_coding: Optional dataset geo-coding.
-    :param xy_names: Optional tuple of the x- and y-coordinate variables in *dataset*. Ignored if *geo_coding* is given.
-    :return: Spatial dataset subset
-    """
-    geo_coding = geo_coding if geo_coding is not None else GeoCoding.from_dataset(dataset, xy_names=xy_names)
-    width, height = geo_coding.size
-    i_min, j_min, i_max, j_max = bbox
-    if i_min > 0 or j_min > 0 or i_max < width - 1 or j_max < height - 1:
-        x_dim, y_dim = geo_coding.dims
-        i_slice = slice(i_min, i_max + 1)
-        j_slice = slice(j_min, j_max + 1)
-        return dataset.isel({x_dim: i_slice, y_dim: j_slice})
-    return dataset
-
-
 def _is_2d_var(var: xr.DataArray, two_d_coord_var: xr.DataArray) -> bool:
     return var.ndim >= 2 and var.shape[-2:] == two_d_coord_var.shape and var.dims[-2:] == two_d_coord_var.dims
 
diff --git a/xcube/core/resample.py b/xcube/core/resample.py
index 9f353b1a..125cabd1 100644
--- a/xcube/core/resample.py
+++ b/xcube/core/resample.py
@@ -24,7 +24,7 @@ from typing import Dict, Any, Sequence, Union
 import xarray as xr
 
 from xcube.core.schema import CubeSchema
-from xcube.core.select import select_vars
+from xcube.core.select import select_variables_subset
 from xcube.core.verify import assert_cube
 
 
@@ -62,7 +62,7 @@ def resample_in_time(cube: xr.Dataset,
         assert_cube(cube)
 
     if var_names:
-        cube = select_vars(cube, var_names)
+        cube = select_variables_subset(cube, var_names)
 
     resampler = cube.resample(skipna=True,
                               closed='left',
diff --git a/xcube/core/select.py b/xcube/core/select.py
index 98111fe3..5aaafe9a 100644
--- a/xcube/core/select.py
+++ b/xcube/core/select.py
@@ -1,9 +1,10 @@
-from typing import Collection
+from typing import Collection, Optional, Tuple
 
 import xarray as xr
+from xcube.core.geocoding import GeoCoding
 
 
-def select_vars(dataset: xr.Dataset, var_names: Collection[str] = None) -> xr.Dataset:
+def select_variables_subset(dataset: xr.Dataset, var_names: Collection[str] = None) -> xr.Dataset:
     """
     Select data variable from given *dataset* and create new dataset.
 
@@ -17,3 +18,39 @@ def select_vars(dataset: xr.Dataset, var_names: Collection[str] = None) -> xr.Da
     if not dropped_variables:
         return dataset
     return dataset.drop_vars(dropped_variables)
+
+
+def select_spatial_subset(dataset: xr.Dataset,
+                          ij_bbox: Tuple[int, int, int, int] = None,
+                          ij_border: int = 0,
+                          xy_bbox: Tuple[float, float, float, float] = None,
+                          xy_border: float = 0.,
+                          geo_coding: GeoCoding = None,
+                          xy_names: Tuple[str, str] = None) -> Optional[xr.Dataset]:
+    """
+    Select a spatial subset of *dataset* for the bounding box *ij_bbox* or *xy_bbox*.
+
+    :param dataset: Source dataset.
+    :param ij_bbox: Bounding box (i_min, i_min, j_max, j_max) in pixel coordinates.
+    :param geo_coding: Optional dataset geo-coding.
+    :param xy_names: Optional tuple of the x- and y-coordinate variables in *dataset*. Ignored if *geo_coding* is given.
+    :return: Spatial dataset subset
+    """
+
+    if ij_bbox is None and xy_bbox is None:
+        raise ValueError('One of ij_bbox and xy_bbox must be given')
+    if ij_bbox and xy_bbox:
+        raise ValueError('Only one of ij_bbox and xy_bbox can be given')
+    geo_coding = geo_coding if geo_coding is not None else GeoCoding.from_dataset(dataset, xy_names=xy_names)
+    if xy_bbox:
+        ij_bbox = geo_coding.ij_bbox(xy_bbox, ij_border=ij_border, xy_border=xy_border)
+        if ij_bbox[0] == -1:
+            return None
+    width, height = geo_coding.size
+    i_min, j_min, i_max, j_max = ij_bbox
+    if i_min > 0 or j_min > 0 or i_max < width - 1 or j_max < height - 1:
+        x_dim, y_dim = geo_coding.dims
+        i_slice = slice(i_min, i_max + 1)
+        j_slice = slice(j_min, j_max + 1)
+        return dataset.isel({x_dim: i_slice, y_dim: j_slice})
+    return dataset
diff --git a/xcube/core/timeseries.py b/xcube/core/timeseries.py
index 19c33a1e..60b62003 100644
--- a/xcube/core/timeseries.py
+++ b/xcube/core/timeseries.py
@@ -27,7 +27,7 @@ import shapely.wkt
 import xarray as xr
 
 from xcube.core.geom import mask_dataset_by_geometry, convert_geometry, GeometryLike, get_dataset_geometry
-from xcube.core.select import select_vars
+from xcube.core.select import select_variables_subset
 from xcube.core.verify import assert_cube
 
 Date = Union[np.datetime64, str]
@@ -78,7 +78,7 @@ def get_time_series(cube: xr.Dataset,
 
     geometry = convert_geometry(geometry)
 
-    dataset = select_vars(cube, var_names)
+    dataset = select_variables_subset(cube, var_names)
     if len(dataset.data_vars) == 0:
         return None
 
diff --git a/xcube/core/xarray.py b/xcube/core/xarray.py
index 08abb4cc..565c3c34 100644
--- a/xcube/core/xarray.py
+++ b/xcube/core/xarray.py
@@ -15,7 +15,7 @@ from xcube.core.extract import get_cube_values_for_points, get_cube_point_indexe
 from xcube.core.level import compute_levels, read_levels, write_levels
 from xcube.core.new import new_cube
 from xcube.core.schema import CubeSchema, get_cube_schema
-from xcube.core.select import select_vars
+from xcube.core.select import select_variables_subset
 from xcube.core.vars2dim import vars_to_dim
 from xcube.core.verify import verify_cube
 
@@ -219,14 +219,14 @@ class DatasetAccessor:
         """
         return verify_cube(self._dataset)
 
-    def select_vars(self, var_names: Sequence[str] = None):
+    def select_variables_subset(self, var_names: Sequence[str] = None):
         """
         Select data variable from given *dataset* and create new dataset.
 
         :param var_names: The names of data variables to select.
         :return: A new dataset. It is empty, if *var_names* is empty. It is *dataset*, if *var_names* is None.
         """
-        return select_vars(self._dataset, var_names)
+        return select_variables_subset(self._dataset, var_names)
 
     def levels(self, **kwargs) -> List[xr.Dataset]:
         """


Test patch:
diff --git a/test/core/gen/test_gen.py b/test/core/gen/test_gen.py
index ef5e4e63..4e73f4b1 100644
--- a/test/core/gen/test_gen.py
+++ b/test/core/gen/test_gen.py
@@ -31,7 +31,7 @@ class DefaultProcessTest(unittest.TestCase):
         status, output = gen_cube_wrapper(
             [get_inputdata_path('20170101-IFR-L4_GHRSST-SSTfnd-ODYSSEA-NWE_002-v2.0-fv1.0.nc')], 'l2c-single.nc')
         self.assertEqual(True, status)
-        self.assertTrue('\nstep 8 of 8: creating input slice in l2c-single.nc...\n' in output)
+        self.assertTrue('\nstep 9 of 9: creating input slice in l2c-single.nc...\n' in output)
         self.assert_cube_ok(xr.open_dataset('l2c-single.nc', autoclose=True),
                             expected_time_dim=1,
                             expected_extra_attrs=dict(date_modified=None,
@@ -46,7 +46,7 @@ class DefaultProcessTest(unittest.TestCase):
                 ('analysed_sst', dict(valid_pixel_expression=None)),
                 ('analysis_error', dict(valid_pixel_expression=None)),
                 ('sea_ice_fraction', dict(valid_pixel_expression=None)),
-                ('water_mask', dict(expression='(mask.sea or mask.lake) and not mask.ice')),
+                ('water_mask', dict(expression='(mask.sea or mask.lake) and not mask.ice', load=True)),
                 ('ice_mask', dict(expression='mask.sea and mask.ice')),
                 ('analysed_sst', dict(valid_pixel_expression='water_mask')),
                 ('analysis_error', dict(valid_pixel_expression='water_mask')),
@@ -59,7 +59,7 @@ class DefaultProcessTest(unittest.TestCase):
             ),
         )
         self.assertEqual(True, status)
-        self.assertTrue('\nstep 8 of 8: creating input slice in l2c-single.nc...\n' in output)
+        self.assertTrue('\nstep 9 of 9: creating input slice in l2c-single.nc...\n' in output)
         self.assert_cube_ok(xr.open_dataset('l2c-single.nc', autoclose=True),
                             expected_time_dim=1,
                             expected_output_vars=('SST', 'SST_uncertainty', 'sea_ice_fraction'),
@@ -72,8 +72,8 @@ class DefaultProcessTest(unittest.TestCase):
             [get_inputdata_path('201701??-IFR-L4_GHRSST-SSTfnd-ODYSSEA-NWE_002-v2.0-fv1.0.nc')], 'l2c.nc',
             no_sort_mode=False)
         self.assertEqual(True, status)
-        self.assertTrue('\nstep 8 of 8: creating input slice in l2c.nc...\n' in output)
-        self.assertTrue('\nstep 8 of 8: appending input slice to l2c.nc...\n' in output)
+        self.assertTrue('\nstep 9 of 9: creating input slice in l2c.nc...\n' in output)
+        self.assertTrue('\nstep 9 of 9: appending input slice to l2c.nc...\n' in output)
         self.assert_cube_ok(xr.open_dataset('l2c.nc', autoclose=True),
                             expected_time_dim=3,
                             expected_extra_attrs=dict(date_modified=None,
@@ -84,7 +84,7 @@ class DefaultProcessTest(unittest.TestCase):
         status, output = gen_cube_wrapper(
             [get_inputdata_path('20170101-IFR-L4_GHRSST-SSTfnd-ODYSSEA-NWE_002-v2.0-fv1.0.nc')], 'l2c-single.zarr')
         self.assertEqual(True, status)
-        self.assertTrue('\nstep 8 of 8: creating input slice in l2c-single.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: creating input slice in l2c-single.zarr...\n' in output)
         self.assert_cube_ok(xr.open_zarr('l2c-single.zarr'),
                             expected_time_dim=1,
                             expected_extra_attrs=dict(date_modified=None,
@@ -96,8 +96,8 @@ class DefaultProcessTest(unittest.TestCase):
             [get_inputdata_path('201701??-IFR-L4_GHRSST-SSTfnd-ODYSSEA-NWE_002-v2.0-fv1.0.nc')], 'l2c.zarr',
             no_sort_mode=False)
         self.assertEqual(True, status)
-        self.assertTrue('\nstep 8 of 8: creating input slice in l2c.zarr...\n' in output)
-        self.assertTrue('\nstep 8 of 8: appending input slice to l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: creating input slice in l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: appending input slice to l2c.zarr...\n' in output)
         self.assert_cube_ok(xr.open_zarr('l2c.zarr'),
                             expected_time_dim=3,
                             expected_extra_attrs=dict(date_modified=None,
@@ -111,9 +111,9 @@ class DefaultProcessTest(unittest.TestCase):
              get_inputdata_path('20170101-IFR-L4_GHRSST-SSTfnd-ODYSSEA-NWE_002-v2.0-fv1.0.nc')], 'l2c.zarr',
             no_sort_mode=True)
         self.assertEqual(True, status)
-        self.assertTrue('\nstep 8 of 8: creating input slice in l2c.zarr...\n' in output)
-        self.assertTrue('\nstep 8 of 8: appending input slice to l2c.zarr...\n' in output)
-        self.assertTrue('\nstep 8 of 8: inserting input slice before index 0 in l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: creating input slice in l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: appending input slice to l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: inserting input slice before index 0 in l2c.zarr...\n' in output)
         self.assert_cube_ok(xr.open_zarr('l2c.zarr'),
                             expected_time_dim=3,
                             expected_extra_attrs=dict(date_modified=None,
@@ -128,9 +128,9 @@ class DefaultProcessTest(unittest.TestCase):
              get_inputdata_path('20170102-IFR-L4_GHRSST-SSTfnd-ODYSSEA-NWE_002-v2.0-fv1.0.nc')], 'l2c.zarr',
             no_sort_mode=True)
         self.assertEqual(True, status)
-        self.assertTrue('\nstep 8 of 8: creating input slice in l2c.zarr...\n' in output)
-        self.assertTrue('\nstep 8 of 8: appending input slice to l2c.zarr...\n' in output)
-        self.assertTrue('\nstep 8 of 8: replacing input slice at index 1 in l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: creating input slice in l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: appending input slice to l2c.zarr...\n' in output)
+        self.assertTrue('\nstep 9 of 9: replacing input slice at index 1 in l2c.zarr...\n' in output)
         self.assert_cube_ok(xr.open_zarr('l2c.zarr'),
                             expected_time_dim=3,
                             expected_extra_attrs=dict(date_modified=None,
diff --git a/test/core/test_evaluate.py b/test/core/test_evaluate.py
index 7b86da82..0900da45 100644
--- a/test/core/test_evaluate.py
+++ b/test/core/test_evaluate.py
@@ -1,8 +1,8 @@
 import unittest
 
+import dask.array as da
 import numpy as np
 import xarray as xr
-
 from xcube.core.evaluate import evaluate_dataset
 
 nan = float('nan')
@@ -17,7 +17,7 @@ class EvaluateDatasetTest(unittest.TestCase):
                                c=((), 1.0),
                                d=((), 1.0, dict(expression='a * b'))),
                           coords=dict(x=(('x',), [1, 2, 3, 4]), y=(('y',), [1, 2])),
-                          attrs=dict(title='test_compute_dataset'))
+                          attrs=dict(title='test_compute_dataset')).chunk(dict(x=2, y=2))
 
     def test_compute_dataset_without_processed_variables(self):
         dataset = self.get_test_dataset()
@@ -52,9 +52,11 @@ class EvaluateDatasetTest(unittest.TestCase):
         dataset = self.get_test_dataset()
         computed_dataset = evaluate_dataset(dataset,
                                             processed_variables=[('a', None),
-                                                                ('b', dict(valid_pixel_expression=None)),
-                                                                ('c', dict(expression='a + b')),
-                                                                ('d', dict(valid_pixel_expression='c > 0.4'))])
+                                                                 ('b', dict(valid_pixel_expression=None)),
+                                                                 ('c', dict(expression='a + b',
+                                                                            load=True)),
+                                                                 ('d', dict(valid_pixel_expression='c > 0.4',
+                                                                            load=True))])
         self.assertIsNot(computed_dataset, dataset)
         self.assertIn('x', computed_dataset)
         self.assertIn('y', computed_dataset)
@@ -68,6 +70,10 @@ class EvaluateDatasetTest(unittest.TestCase):
         self.assertEqual((2, 4), computed_dataset.a.shape)
         self.assertEqual((2, 4), computed_dataset.b.shape)
         self.assertEqual((2, 4), computed_dataset.c.shape)
+        self.assertIsInstance(computed_dataset.a.data, da.Array)
+        self.assertIsInstance(computed_dataset.b.data, da.Array)
+        self.assertIsInstance(computed_dataset.c.data, np.ndarray)  # load=True --> load c as numpy array
+        self.assertIsInstance(computed_dataset.d.data, np.ndarray)  # load=True --> load d as numpy array
         self.assertIn('expression', computed_dataset.c.attrs)
         self.assertEqual((2, 4), computed_dataset.d.shape)
         self.assertIn('expression', computed_dataset.d.attrs)
diff --git a/test/core/test_select.py b/test/core/test_select.py
index 5ec77a73..5fe41ee7 100644
--- a/test/core/test_select.py
+++ b/test/core/test_select.py
@@ -1,27 +1,74 @@
 import unittest
 
 from test.sampledata import create_highroc_dataset
-from xcube.core.select import select_vars
+from xcube.core.select import select_spatial_subset, select_variables_subset
 
 
-class SelectVariablesTest(unittest.TestCase):
-    def test_select_all(self):
+class SelectVariablesSubsetTest(unittest.TestCase):
+    def test_select_variables_subset_all(self):
         ds1 = create_highroc_dataset()
         # noinspection PyTypeChecker
-        ds2 = select_vars(ds1, None)
+        ds2 = select_variables_subset(ds1, None)
         self.assertIs(ds2, ds1)
-        ds2 = select_vars(ds1, ds1.data_vars.keys())
+        ds2 = select_variables_subset(ds1, ds1.data_vars.keys())
         self.assertIs(ds2, ds1)
 
-    def test_select_none(self):
+    def test_select_variables_subset_none(self):
         ds1 = create_highroc_dataset()
-        ds2 = select_vars(ds1, [])
+        ds2 = select_variables_subset(ds1, [])
         self.assertEqual(0, len(ds2.data_vars))
-        ds2 = select_vars(ds1, ['bibo'])
+        ds2 = select_variables_subset(ds1, ['bibo'])
         self.assertEqual(0, len(ds2.data_vars))
 
-    def test_select_variables_for_some(self):
+    def test_select_variables_subset_some(self):
         ds1 = create_highroc_dataset()
         self.assertEqual(36, len(ds1.data_vars))
-        ds2 = select_vars(ds1, ['conc_chl', 'c2rcc_flags', 'rtoa_10'])
+        ds2 = select_variables_subset(ds1, ['conc_chl', 'c2rcc_flags', 'rtoa_10'])
         self.assertEqual(3, len(ds2.data_vars))
+
+
+class SelectSpatialSubsetTest(unittest.TestCase):
+    def test_select_spatial_subset_all_ij_bbox(self):
+        ds1 = create_highroc_dataset()
+        ds2 = select_spatial_subset(ds1, ij_bbox=(0, 0, 4, 3))
+        self.assertIs(ds2, ds1)
+
+    def test_select_spatial_subset_some_ij_bbox(self):
+        ds1 = create_highroc_dataset()
+        ds2 = select_spatial_subset(ds1, ij_bbox=(1, 1, 4, 3))
+        self.assertEqual((2, 3), ds2.conc_chl.shape)
+
+    def test_select_spatial_subset_none_ij_bbox(self):
+        ds1 = create_highroc_dataset()
+        ds2 = select_spatial_subset(ds1, ij_bbox=(5, 6, 7, 8))
+        self.assertEqual(None, ds2)
+        ds2 = select_spatial_subset(ds1, ij_bbox=(-6, -4, 2, 2))
+        self.assertEqual(None, ds2)
+
+    def test_select_spatial_subset_all_xy_bbox(self):
+        ds1 = create_highroc_dataset()
+        ds2 = select_spatial_subset(ds1, xy_bbox=(7.9, 53.9, 12., 56.4))
+        self.assertIs(ds2, ds1)
+
+    def test_select_spatial_subset_some_xy_bbox(self):
+        ds1 = create_highroc_dataset()
+        ds2 = select_spatial_subset(ds1, xy_bbox=(8., 55, 10., 56.))
+        self.assertEqual((2, 2), ds2.conc_chl.shape)
+
+    def test_select_spatial_subset_none_xy_bbox(self):
+        ds1 = create_highroc_dataset()
+        ds2 = select_spatial_subset(ds1, xy_bbox=(13., 57., 15., 60.))
+        self.assertEqual(None, ds2)
+        ds2 = select_spatial_subset(ds1, xy_bbox=(5.5, 55, 6.5, 56))
+        self.assertEqual(None, ds2)
+
+    def test_select_spatial_subset_invalid_params(self):
+        ds1 = create_highroc_dataset()
+        with self.assertRaises(ValueError) as cm:
+            select_spatial_subset(ds1, ij_bbox=(5, 6, 7, 8), xy_bbox=(0., 0., 1., 2.))
+        self.assertEqual("Only one of ij_bbox and xy_bbox can be given", f'{cm.exception}')
+
+        with self.assertRaises(ValueError) as cm:
+            select_spatial_subset(ds1)
+        self.assertEqual("One of ij_bbox and xy_bbox must be given", f'{cm.exception}')
+
diff --git a/test/core/test_xarray.py b/test/core/test_xarray.py
index 295c914e..19db6389 100644
--- a/test/core/test_xarray.py
+++ b/test/core/test_xarray.py
@@ -20,7 +20,7 @@ class XCubeDatasetAccessorTest(unittest.TestCase):
         self.assertIsInstance(xr.Dataset.xcube.new(), xr.Dataset)
 
     def test_select_vars(self):
-        self.assertIsInstance(xr.Dataset().xcube.select_vars(), xr.Dataset)
+        self.assertIsInstance(xr.Dataset().xcube.select_variables_subset(), xr.Dataset)
 
     def test_dump(self):
         self.assertIsInstance(xr.Dataset().xcube.dump(), str)