Repository: dcs4cop/xcube
Pull request #108 (2019-06-21 16:47:48)

Forman 102 improve extract
To test it:

    $ conda activate xcube
    $ xcube extract --help
    $ xcube extract xcube/webapi/res/demo/cube.nc test/cli/extract-points/all-valid-points.csv 
    $ xcube extract xcube/webapi/res/demo/cube.nc test/cli/extract-points/all-valid-points.csv -CIRB
    $ xcube extract xcube/webapi/res/demo/cube.nc test/cli/extract-points/all-valid-points.csv -CIRB -o test.csv


Patch:
diff --git a/CHANGES.md b/CHANGES.md
index 98c7eb11..59eb8efb 100644
--- a/CHANGES.md
+++ b/CHANGES.md
@@ -2,6 +2,8 @@
 
 ### Enhancements
 
+* Reworked `xcube extract` to be finally useful and effective for point data extraction. (#102) 
+* `xcube server`can now filter datasets by point coordinate, e.g. `/datasets?point=12.5,52.8`. (#50) 
 * `xcube server`can now limit time series to a maximum number of 
   valid (not NaN) values. To activate, pass optional query parameter `maxValids` to the various `/ts`
   functions. The special value `-1` will restrict the result to contain only valid values. (#113) 
diff --git a/README.md b/README.md
index 5b36bc78..f99c7232 100644
--- a/README.md
+++ b/README.md
@@ -123,7 +123,8 @@ To start a demo using docker use the following commands
     
 
 ## `xcube chunk`
-Can be used for changing the chunks of an existing data cube. 
+
+(Re-)chunk dataset.
 
     $ xcube chunk --help
     Usage: xcube chunk [OPTIONS] <input> <output>
@@ -146,7 +147,8 @@ Example:
     $ xcube chunk input_not_chunked.zarr output_rechunked.zarr --chunks "time=1,lat=270,lon=270"
 
 ## `xcube dump`
-Can be used for printing out the metadata of a data cube. 
+
+Dump contents of a dataset.
 
     $ xcube dump --help
     Usage: xcube dump [OPTIONS] <path>
@@ -164,30 +166,44 @@ Example:
     $ xcube dump xcube_cube.zarr 
 
 ## `xcube extract`
-A time series for a specific location can be requested using `xcube extract`.
+
+Extract cube points.
 
     $ xcube dump --help
-    Usage: xcube extract [OPTIONS] <cube> <coords>
-    
-      Extract cube time series. Extracts data from <cube> at points given by
-      coordinates <coords> and writes the resulting time series to <output>.
-    
+    Usage: xcube extract [OPTIONS] CUBE POINTS
+
+      Extract cube points.
+
+      Extracts data cells from CUBE at coordinates given in each POINTS record
+      and writes the resulting values to given output path and format.
+
+      <points> must be a CSV file that provides at least the columns "lon",
+      "lat", and "time". The "lon" and "lat" columns provide a point's location
+      in decimal degrees. The "time" column provides a point's date or date-
+      time. Its format should preferably be ISO, but other formats may work as
+      well.
+
     Options:
-      -i, --indexes          Include indexes in output.
-      -o, --output <output>  Output file.
-      -f, --format <format>  Format of the output. If not given, guessed from
-                             <output>, otherwise <stdout> is used.
-      -p, --params <params>  Parameters specific for the output format. Comma-
-                             separated list of <key>=<value> pairs.
-      --help                 Show this message and exit.
+      -o, --output TEXT             Output file. If omitted, output is written to
+                                    stdout.
+      -f, --format [csv|json|xlsx]  Output format. Currently, only 'csv' is
+                                    supported.
+      -C, --coords                  Include cube cell coordinates in output.
+      -B, --bounds                  Include cube cell coordinate boundaries (if
+                                    any) in output.
+      -I, --indexes                 Include cube cell indexes in output.
+      -R, --refs                    Include point values as reference in output.
+      --help                        Show this message and exit.
+
 
-Example: # TODO: Help is needed here - how are the coords passed by the user? 
+Example:  
     
-    $ xcube extract xcube_cube.zarr 
+    $ xcube extract xcube_cube.zarr point_data.csv -CBIR
     
     
 ## `xcube gen`
-Is used to generate data cubes.
+
+Generate data cube.
 
     $ xcube gen --help
     Usage: xcube gen [OPTIONS] INPUT_FILES
@@ -275,7 +291,8 @@ Available xcube input processors within xcube's organisation:
 
 
 ## `xcube grid`
-For choosing a suitable gridding for a data cube, `xbube grid` is a useful tool.
+
+Find spatial data cube resolutions and adjust bounding boxes.
 
     $ xcube grid --help
     Usage: xcube grid [OPTIONS] COMMAND [ARGS]...
@@ -449,6 +466,8 @@ Example:
 
 ## `xcube level`
 
+Generate multi-resolution levels.
+
     $ xcube level --help
     Usage: xcube level [OPTIONS] <input>
     
@@ -483,6 +502,8 @@ Example:
 
 ## `xcube prune`
 
+Delete empty chunks.
+
     $ xcube prune --help
     Usage: xcube prune [OPTIONS] INPUT
     
@@ -497,6 +518,8 @@ Example:
 
 ## `xcube vars2dim`
 
+Convert cube variables into new dimension.
+
     $ xcube vars2dim --help
     Usage: xcube vars2dim [OPTIONS] <cube>
     
@@ -515,9 +538,10 @@ Example:
       --help                     Show this message and exit.
 
 
-
 ## `xcube resample`
 
+Resample data along the time dimension.
+
     $ xcube resample --help
     Usage: xcube resample [OPTIONS] INPUT
     
@@ -589,8 +613,7 @@ Downsampling example:
 
 ## `xcube serve`
 
-Is a light-weight web server that provides various services based on 
-xcube data cubes. 
+Serve data cubes via web service. 
 
     $ xcube serve --help
     Usage: xcube serve [OPTIONS] CUBE...
diff --git a/xcube/api/__init__.py b/xcube/api/__init__.py
index 7d449db6..a3d4c726 100644
--- a/xcube/api/__init__.py
+++ b/xcube/api/__init__.py
@@ -1,10 +1,13 @@
 from .api import XCubeAPI
 from .chunk import chunk_dataset
 from .dump import dump_dataset
-from .extract import get_cube_point_indexes, get_cube_values_for_indexes, get_cube_values_for_points, \
-    get_dataset_indexes
+from .extract import get_cube_values_for_points, get_cube_point_indexes, get_cube_values_for_indexes, \
+    get_dataset_indexes, DEFAULT_INDEX_NAME_PATTERN, DEFAULT_REF_NAME_PATTERN, DEFAULT_INTERP_POINT_METHOD
+from .levels import compute_levels, read_levels, write_levels
 from .new import new_cube
 from .readwrite import read_cube, open_cube, open_dataset, read_dataset, write_dataset
+from .resample import resample_in_time
+from .select import select_vars
 from .ts import get_time_series
 from .vars_to_dim import vars_to_dim
 from .verify import assert_cube, verify_cube
diff --git a/xcube/api/api.py b/xcube/api/api.py
index e426f731..a903f1cc 100644
--- a/xcube/api/api.py
+++ b/xcube/api/api.py
@@ -4,16 +4,19 @@ from typing import Dict, List, Mapping, Any, Union, Sequence
 import numpy as np
 import pandas as pd
 import xarray as xr
-from xcube.api.levels import compute_levels
 
-from .vars_to_dim import vars_to_dim
 from .chunk import chunk_dataset
+# noinspection PyUnresolvedReferences
+from .compute import compute_dataset
 from .dump import dump_dataset
 from .extract import get_cube_values_for_points, get_cube_point_indexes, get_cube_values_for_indexes, \
-    get_dataset_indexes, INDEX_NAME_PATTERN
+    get_dataset_indexes, DEFAULT_INDEX_NAME_PATTERN
+# noinspection PyUnresolvedReferences
+from .levels import compute_levels, read_levels, write_levels
 from .new import new_cube
 from .readwrite import read_cube, open_cube, write_cube
 from .select import select_vars
+from .vars_to_dim import vars_to_dim
 from .verify import verify_cube
 
 
@@ -39,6 +42,7 @@ class XCubeAPI:
             time_periods=5,
             time_freq="D",
             time_start='2010-01-01T00:00:00',
+            inverse_lat=False,
             drop_bounds=False,
             variables=None) -> xr.Dataset:
         """
@@ -53,6 +57,7 @@ class XCubeAPI:
         :param time_periods: Number of time steps
         :param time_freq: Duration of each time step
         :param time_start: First time value
+        :param inverse_lat: Whether to create an inverse latitude axis
         :param drop_bounds: If True, coordinate bounds variables are not created.
         :param variables: Dictionary of data variables to be added.
         :return: A cube instance
@@ -66,6 +71,7 @@ class XCubeAPI:
                         time_periods=time_periods,
                         time_freq=time_freq,
                         time_start=time_start,
+                        inverse_lat=inverse_lat,
                         drop_bounds=drop_bounds,
                         variables=variables)
 
@@ -110,7 +116,7 @@ class XCubeAPI:
     def values_for_points(self,
                           points: Union[xr.Dataset, pd.DataFrame, Mapping[str, Any]],
                           var_names: Sequence[str] = None,
-                          index_name_pattern: str = INDEX_NAME_PATTERN,
+                          index_name_pattern: str = DEFAULT_INDEX_NAME_PATTERN,
                           include_indexes: bool = False,
                           method: str = "nearest",
                           cube_asserted: bool = False):
@@ -137,7 +143,7 @@ class XCubeAPI:
     def values_for_indexes(self,
                            indexes: Union[xr.Dataset, pd.DataFrame, Mapping[str, Any]],
                            var_names: Sequence[str] = None,
-                           index_name_pattern: str = INDEX_NAME_PATTERN,
+                           index_name_pattern: str = DEFAULT_INDEX_NAME_PATTERN,
                            method: str = "nearest",
                            cube_asserted: bool = False) -> xr.Dataset:
         """
@@ -153,7 +159,7 @@ class XCubeAPI:
         """
         return get_cube_values_for_indexes(self._dataset,
                                            indexes,
-                                           var_names=var_names,
+                                           data_var_names=var_names,
                                            index_name_pattern=index_name_pattern,
                                            method=method,
                                            cube_asserted=cube_asserted)
@@ -161,8 +167,8 @@ class XCubeAPI:
     def point_indexes(self,
                       points: Union[xr.Dataset, pd.DataFrame, Mapping[str, Any]],
                       dim_name_mapping: Mapping[str, str] = None,
-                      index_name_pattern: str = INDEX_NAME_PATTERN,
-                      dtype=np.float64,
+                      index_name_pattern: str = DEFAULT_INDEX_NAME_PATTERN,
+                      index_dtype=np.float64,
                       cube_asserted: bool = False):
         """
         Get indexes of given coordinates in *points* into this cube.
@@ -171,23 +177,24 @@ class XCubeAPI:
         :param dim_name_mapping: A mapping from dimension names in *cube* to column names in *points*.
         :param index_name_pattern: A naming pattern for the computed indexes columns.
                Must include "{name}" which will be replaced by the dimension name.
+        :param index_dtype: Numpy data type for the indexes. If it is a floating point type (default),
+               then *indexes* will contain fractions, which may be used for interpolation.
+               For out-of-range coordinates in *points*, indexes will be -1 if *index_dtype* is an integer type, and NaN,
+               if *index_dtype* is a floating point types.
         :param cube_asserted: If False, *cube* will be verified, otherwise it is expected to be a valid cube.
-        :param dtype: Numpy data type for the indexes. If it is floating point type (default),
-               then *indexes* contain fractions, which may be used for interpolation. If *dtype* is an integer
-               type out-of-range coordinates are indicated by index -1, and NaN if it is is a floating point type.
         :return: A dataset containing the index columns.
         """
         return get_cube_point_indexes(self._dataset,
                                       points,
                                       dim_name_mapping=dim_name_mapping,
                                       index_name_pattern=index_name_pattern,
-                                      dtype=dtype,
+                                      index_dtype=index_dtype,
                                       cube_asserted=cube_asserted)
 
     def indexes(self,
                 coord_var_name: str,
                 coord_values: Union[xr.DataArray, np.ndarray],
-                dtype=np.float64) -> Union[xr.DataArray, np.ndarray]:
+                index_dtype=np.float64) -> Union[xr.DataArray, np.ndarray]:
         """
         Compute the indexes into a coordinate variable *coord_var_name* of this cube
         for the given coordinate values *coord_values*.
@@ -196,22 +203,23 @@ class XCubeAPI:
         otherwise the result will be nonsense.
 
         For any value in *coord_values* that is out of the bounds of the coordinate variable's values,
-        the index depends on the value of *dtype*. If *dtype* is an integer type, invalid indexes are
+        the index depends on the value of *index_dtype*. If *index_dtype* is an integer type, invalid indexes are
         encoded as -1 while for floating point types, NaN will be used.
 
         Returns the indexes as an array-like object of type *dtype*.
 
         :param coord_var_name: Name of a coordinate variable.
         :param coord_values: Array-like coordinate values.
-        :param dtype: (Numpy) data type for the indexes. If it is floating point type (default),
-               then *indexes* contain fractions, which may be used for interpolation. If *dtype* is an integer
-               type out-of-range coordinates are indicated by index -1, and NaN if it is is a floating point type.
+        :param index_dtype: Numpy data type for the indexes. If it is a floating point type (default),
+               then *indexes* will contain fractions, which may be used for interpolation.
+               For out-of-range coordinates in *points*, indexes will be -1 if *index_dtype* is an integer type, and NaN,
+               if *index_dtype* is a floating point types.
         :return: The indexes and their fractions as a tuple of numpy int64 and float64 arrays.
         """
         return get_dataset_indexes(self._dataset,
                                    coord_var_name=coord_var_name,
                                    coord_values=coord_values,
-                                   dtype=dtype)
+                                   index_dtype=index_dtype)
 
     def chunk(self, chunk_sizes: Dict[str, int] = None, format_name: str = None) -> xr.Dataset:
         """
@@ -285,3 +293,19 @@ class XCubeAPI:
         :return: A list of dataset instances representing the multi-level pyramid.
         """
         return compute_levels(self._dataset, **kwargs)
+
+    def resample_in_time(self):
+        """
+        Resample a data cube in the time dimension.
+
+        :param frequency: Resampling frequency.
+        :param method: Resampling method or sequence of resampling methods.
+        :param offset: Offset used to adjust the resampled time labels. Some pandas date offset strings are supported.
+        :param base: Resampling method.
+        :param var_names: Variable names to include.
+        :param tolerance: Time tolerance for selective upsampling methods. Defaults to *frequency*.
+        :param interp_kind: Kind of interpolation if *method* is 'interpolation'.
+        :param metadata: Output metadata.
+        :return: A new data cube resampled in time.
+        """
+        return resample_in_time()
\ No newline at end of file
diff --git a/xcube/api/extract.py b/xcube/api/extract.py
index 8e43105e..d342ac06 100644
--- a/xcube/api/extract.py
+++ b/xcube/api/extract.py
@@ -6,8 +6,8 @@ import xarray as xr
 
 from .verify import assert_cube
 
-INDEX_NAME_PATTERN = '{name}_index'
-FRACTION_NAME_PATTERN = '{name}_fraction'
+DEFAULT_INDEX_NAME_PATTERN = '{name}_index'
+DEFAULT_REF_NAME_PATTERN = '{name}_ref'
 INDEX_DIM_NAME = "idx"
 
 POINT_INTERP_METHOD_NEAREST = "nearest"
@@ -18,8 +18,12 @@ DEFAULT_INTERP_POINT_METHOD = POINT_INTERP_METHOD_NEAREST
 def get_cube_values_for_points(cube: xr.Dataset,
                                points: Union[xr.Dataset, pd.DataFrame, Mapping[str, Any]],
                                var_names: Sequence[str] = None,
-                               index_name_pattern: str = INDEX_NAME_PATTERN,
+                               include_coords: bool = False,
+                               include_bounds: bool = False,
                                include_indexes: bool = False,
+                               index_name_pattern: str = DEFAULT_INDEX_NAME_PATTERN,
+                               include_refs: bool = False,
+                               ref_name_pattern: str = DEFAULT_REF_NAME_PATTERN,
                                method: str = DEFAULT_INTERP_POINT_METHOD,
                                cube_asserted: bool = False) -> xr.Dataset:
     """
@@ -28,9 +32,14 @@ def get_cube_values_for_points(cube: xr.Dataset,
     :param cube: The cube dataset.
     :param points: Dictionary that maps dimension name to coordinate arrays.
     :param var_names: An optional list of names of data variables in *cube* whose values shall be extracted.
-    :param index_name_pattern: A naming pattern for the computed indexes columns.
-           Must include "{name}" which will be replaced by the dimension name.
-    :param include_indexes: Weather to include computed indexes in return value.
+    :param include_coords: Weather to include the cube coordinates for each point in return value.
+    :param include_bounds: Weather to include the cube coordinate boundaries (if any) for each point in return value.
+    :param include_indexes: Weather to include computed indexes into the cube for each point in return value.
+    :param index_name_pattern: A naming pattern for the computed index columns.
+           Must include "{name}" which will be replaced by the index' dimension name.
+    :param include_refs: Weather to include point (reference) values in return value.
+    :param ref_name_pattern: A naming pattern for the computed point data columns.
+           Must include "{name}" which will be replaced by the point's attribute name.
     :param method: "nearest" or "linear".
     :param cube_asserted: If False, *cube* will be verified, otherwise it is expected to be a valid cube.
     :return: A new data frame whose columns are values from *cube* variables at given *points*.
@@ -41,12 +50,14 @@ def get_cube_values_for_points(cube: xr.Dataset,
     point_indexes = get_cube_point_indexes(cube,
                                            points,
                                            index_name_pattern=index_name_pattern,
-                                           dtype=np.int64 if method == POINT_INTERP_METHOD_NEAREST else np.float64,
+                                           index_dtype=np.int64 if method == POINT_INTERP_METHOD_NEAREST else np.float64,
                                            cube_asserted=True)
 
     cube_values = get_cube_values_for_indexes(cube,
                                               point_indexes,
-                                              var_names=var_names,
+                                              include_coords,
+                                              include_bounds,
+                                              data_var_names=var_names,
                                               index_name_pattern=index_name_pattern,
                                               method=method,
                                               cube_asserted=True)
@@ -54,13 +65,20 @@ def get_cube_values_for_points(cube: xr.Dataset,
     if include_indexes:
         cube_values.update(point_indexes)
 
+    if include_refs:
+        point_refs = xr.Dataset({ref_name_pattern.format(name=name): xr.DataArray(points[name], dims=[INDEX_DIM_NAME])
+                                 for name in points.keys()})
+        cube_values.update(point_refs)
+
     return cube_values
 
 
 def get_cube_values_for_indexes(cube: xr.Dataset,
                                 indexes: Union[xr.Dataset, pd.DataFrame, Mapping[str, Any]],
-                                var_names: Sequence[str] = None,
-                                index_name_pattern: str = INDEX_NAME_PATTERN,
+                                include_coords: bool = False,
+                                include_bounds: bool = False,
+                                data_var_names: Sequence[str] = None,
+                                index_name_pattern: str = DEFAULT_INDEX_NAME_PATTERN,
                                 method: str = DEFAULT_INTERP_POINT_METHOD,
                                 cube_asserted: bool = False) -> xr.Dataset:
     """
@@ -68,7 +86,9 @@ def get_cube_values_for_indexes(cube: xr.Dataset,
 
     :param cube: A cube dataset.
     :param indexes: A mapping from column names to index and fraction arrays for all cube dimensions.
-    :param var_names: An optional list of names of data variables in *cube* whose values shall be extracted.
+    :param include_coords: Weather to include the cube coordinates for each point in return value.
+    :param include_bounds: Weather to include the cube coordinate boundaries (if any) for each point in return value.
+    :param data_var_names: An optional list of names of data variables in *cube* whose values shall be extracted.
     :param index_name_pattern: A naming pattern for the computed indexes columns.
            Must include "{name}" which will be replaced by the dimension name.
     :param method: "nearest" or "linear".
@@ -83,35 +103,44 @@ def get_cube_values_for_indexes(cube: xr.Dataset,
     if method != POINT_INTERP_METHOD_NEAREST:
         raise NotImplementedError(f"method {method!r} not yet implemented")
 
-    all_var_names = tuple(cube.data_vars.keys())
-    if len(all_var_names) == 0:
+    all_data_var_names = tuple(cube.data_vars.keys())
+    if len(all_data_var_names) == 0:
         raise ValueError("cube is empty")
 
-    if var_names is not None:
-        if len(var_names) == 0:
+    if data_var_names is not None:
+        if len(data_var_names) == 0:
             return xr.Dataset(coords=indexes.coords if hasattr(indexes, "coords") else None)
-        for var_name in var_names:
+        for var_name in data_var_names:
             if var_name not in cube.data_vars:
                 raise ValueError(f"variable {var_name!r} not found in cube")
     else:
-        var_names = all_var_names
+        data_var_names = all_data_var_names
 
-    dim_names = cube[var_names[0]].dims
+    dim_names = cube[data_var_names[0]].dims
     num_dims = len(dim_names)
     index_names = [index_name_pattern.format(name=dim_name) for dim_name in dim_names]
     num_points = _validate_points(indexes, index_names, param_name="indexes")
     indexes = _normalize_points(indexes)
 
-    # Flatten any coordinate bounds variables
-    cube = xr.Dataset({var_name: cube[var_name] for var_name in var_names}, coords=cube.coords)
+    cube = xr.Dataset({var_name: cube[var_name] for var_name in data_var_names}, coords=cube.coords)
+
+    new_bounds_vars = {}
     bounds_var_names = _get_coord_bounds_var_names(cube)
-    new_vars = {}
-    for var_name, bnds_var_name in bounds_var_names.items():
-        bnds_var = cube[bnds_var_name]
-        new_vars[var_name + "_start"] = bnds_var[:, 0]
-        new_vars[var_name + "_end"] = bnds_var[:, 1]
-    cube = cube.drop(list(bounds_var_names.keys()))
-    cube = cube.assign_coords(**new_vars)
+    drop_coords = None
+    if bounds_var_names:
+        if include_bounds:
+            # Flatten any coordinate bounds variables
+            for var_name, bnds_var_name in bounds_var_names.items():
+                bnds_var = cube[bnds_var_name]
+                new_bounds_vars[var_name + "_lower"] = bnds_var[:, 0]
+                new_bounds_vars[var_name + "_upper"] = bnds_var[:, 1]
+            cube = cube.assign_coords(**new_bounds_vars)
+        cube = cube.drop(bounds_var_names.values())
+        if not include_coords:
+            drop_coords = set(cube.coords).difference(new_bounds_vars.keys())
+    else:
+        if not include_coords:
+            drop_coords = set(cube.coords)
 
     # Generate a validation condition so we can filter out invalid rows (where any index == -1)
     is_valid_point = None
@@ -130,10 +159,10 @@ def get_cube_values_for_indexes(cube: xr.Dataset,
         cube_values = cube.isel(cube_selector)
     elif num_valid_points == 0:
         # All indexes are invalid
-        new_vars = {}
-        for var_name in var_names:
-            new_vars[var_name] = _empty_points_var(cube[var_name], num_points)
-        cube_values = xr.Dataset(new_vars)
+        new_bounds_vars = {}
+        for var_name in cube.variables:
+            new_bounds_vars[var_name] = _empty_points_var(cube[var_name], num_points)
+        cube_values = xr.Dataset(new_bounds_vars)
     else:
         # Some invalid indexes
         idx = np.arange(num_points)
@@ -144,14 +173,17 @@ def get_cube_values_for_indexes(cube: xr.Dataset,
         cube_selector = {dim_names[i]: good_indexes[index_names[i]] for i in range(num_dims)}
         cube_values = cube.isel(cube_selector)
 
-        new_vars = {}
-        for var_name in var_names:
+        new_bounds_vars = {}
+        for var_name in cube.variables:
             var = cube_values[var_name]
             new_var = _empty_points_var(var, num_points)
             new_var[good_idx] = var
-            new_vars[var_name] = new_var
+            new_bounds_vars[var_name] = new_var
+
+        cube_values = xr.Dataset(new_bounds_vars)
 
-        cube_values = xr.Dataset(new_vars)
+    if drop_coords:
+        cube_values = cube_values.drop(drop_coords)
 
     return cube_values
 
@@ -159,8 +191,8 @@ def get_cube_values_for_indexes(cube: xr.Dataset,
 def get_cube_point_indexes(cube: xr.Dataset,
                            points: Union[xr.Dataset, pd.DataFrame, Mapping[str, Any]],
                            dim_name_mapping: Mapping[str, str] = None,
-                           index_name_pattern: str = INDEX_NAME_PATTERN,
-                           dtype=np.float64,
+                           index_name_pattern: str = DEFAULT_INDEX_NAME_PATTERN,
+                           index_dtype=np.float64,
                            cube_asserted: bool = False) -> xr.Dataset:
     """
     Get indexes of given point coordinates *points* into the given *dataset*.
@@ -170,9 +202,10 @@ def get_cube_point_indexes(cube: xr.Dataset,
     :param dim_name_mapping: A mapping from dimension names in *cube* to column names in *points*.
     :param index_name_pattern: A naming pattern for the computed indexes columns.
            Must include "{name}" which will be replaced by the dimension name.
-    :param dtype: Numpy data type for the indexes. If it is floating point type (default),
-           then *indexes* contain fractions, which may be used for interpolation. If *dtype* is an integer
-           type out-of-range coordinates are indicated by index -1, and NaN if it is is a floating point type.
+    :param index_dtype: Numpy data type for the indexes. If it is a floating point type (default),
+           then *indexes* will contain fractions, which may be used for interpolation.
+           For out-of-range coordinates in *points*, indexes will be -1 if *index_dtype* is an integer type, and NaN,
+           if *index_dtype* is a floating point types.
     :param cube_asserted: If False, *cube* will be verified, otherwise it is expected to be a valid cube.
     :return: A dataset containing the index columns.
     """
@@ -188,7 +221,7 @@ def get_cube_point_indexes(cube: xr.Dataset,
     indexes = []
     for dim_name, col_name in zip(dim_names, col_names):
         col = points[col_name]
-        coord_indexes = get_dataset_indexes(cube, dim_name, col, dtype=dtype)
+        coord_indexes = get_dataset_indexes(cube, dim_name, col, index_dtype=index_dtype)
         indexes.append((index_name_pattern.format(name=dim_name),
                         xr.DataArray(coord_indexes, dims=[INDEX_DIM_NAME])))
 
@@ -198,7 +231,7 @@ def get_cube_point_indexes(cube: xr.Dataset,
 def get_dataset_indexes(dataset: xr.Dataset,
                         coord_var_name: str,
                         coord_values: Union[xr.DataArray, np.ndarray],
-                        dtype=np.float64) -> Union[xr.DataArray, np.ndarray]:
+                        index_dtype=np.float64) -> Union[xr.DataArray, np.ndarray]:
     """
     Compute the indexes and their fractions into a coordinate variable *coord_var_name* of a *dataset*
     for the given coordinate values *coord_values*.
@@ -207,14 +240,15 @@ def get_dataset_indexes(dataset: xr.Dataset,
     otherwise the result will be nonsense.
 
     For any value in *coord_values* that is out of the bounds of the coordinate variable's values,
-    the index will be -1 and the fraction will be NaN.
+    the index depends on the value of *index_dtype*. If *index_dtype* is an integer type, invalid indexes are
+    encoded as -1 while for floating point types, NaN will be used.
 
     Returns a tuple of indexes as int64 array and fractions as float64 array.
 
     :param dataset: A cube dataset.
     :param coord_var_name: Name of a coordinate variable.
     :param coord_values: Array-like coordinate values.
-    :param dtype: Numpy data type for the indexes. If it is floating point type (default),
+    :param index_dtype: Numpy data type for the indexes. If it is floating point type (default),
            then *indexes* contain fractions, which may be used for interpolation. If *dtype* is an integer
            type out-of-range coordinates are indicated by index -1, and NaN if it is is a floating point type.
     :return: The indexes and their fractions as a tuple of numpy int64 and float64 arrays.
@@ -228,8 +262,8 @@ def get_dataset_indexes(dataset: xr.Dataset,
         coord_bounds = coord_bounds_var.values
         if np.issubdtype(coord_bounds.dtype, np.datetime64):
             coord_bounds = coord_bounds.astype(np.uint64)
-        is_reversed = (coord_bounds[0, 1] - coord_bounds[0, 0]) < 0
-        if is_reversed:
+        is_inverse = (coord_bounds[0, 1] - coord_bounds[0, 0]) < 0
+        if is_inverse:
             coord_bounds = coord_bounds[::-1, ::-1]
         coords = np.zeros(n2, dtype=coord_bounds.dtype)
         coords[0:-1] = coord_bounds[:, 0]
@@ -238,8 +272,8 @@ def get_dataset_indexes(dataset: xr.Dataset,
         center_coords = coord_var.values
         if np.issubdtype(center_coords.dtype, np.datetime64):
             center_coords = center_coords.astype(np.uint64)
-        is_reversed = (center_coords[-1] - center_coords[0]) < 0
-        if is_reversed:
+        is_inverse = (center_coords[-1] - center_coords[0]) < 0
+        if is_inverse:
             center_coords = center_coords[::-1]
         deltas = np.zeros(n2, dtype=center_coords.dtype)
         deltas[0:-2] = np.diff(center_coords)
@@ -265,16 +299,15 @@ def get_dataset_indexes(dataset: xr.Dataset,
 
     indexes = np.linspace(0.0, n1, n2, dtype=np.float64)
     interp_indexes = np.interp(coord_values, coords, indexes, left=-1, right=-1)
-    if np.issubdtype(dtype, np.integer):
-        interp_indexes = interp_indexes.astype(dtype)
-        interp_indexes[interp_indexes >= n1] = n1 - 1
+    if is_inverse:
+        i = interp_indexes >= 0
+        interp_indexes[i] = n1 - interp_indexes[i]
+    interp_indexes[interp_indexes >= n1] = n1 - 1e-9
+    if np.issubdtype(index_dtype, np.integer):
+        interp_indexes = interp_indexes.astype(index_dtype)
     else:
-        interp_indexes[interp_indexes >= n1] = n1 - 1e-9
         interp_indexes[interp_indexes < 0] = np.nan
 
-    if is_reversed:
-        interp_indexes = interp_indexes[::-1]
-
     return interp_indexes
 
 
diff --git a/xcube/api/new.py b/xcube/api/new.py
index 6a70ffe5..d0ee3f1b 100644
--- a/xcube/api/new.py
+++ b/xcube/api/new.py
@@ -16,6 +16,7 @@ def new_cube(title="Test Cube",
              time_periods=5,
              time_freq="D",
              time_start='2010-01-01T00:00:00',
+             inverse_lat=False,
              drop_bounds=False,
              variables=None):
     """
@@ -30,6 +31,7 @@ def new_cube(title="Test Cube",
     :param time_periods: Number of time steps
     :param time_freq: Duration of each time step
     :param time_start: First time value
+    :param inverse_lat: Whether to create an inverse latitude axis
     :param drop_bounds: If True, coordinate bounds variables are not created.
     :param variables: Dictionary of data variables to be added.
     :return: A cube instance
@@ -43,12 +45,16 @@ def new_cube(title="Test Cube",
     if time_periods < 0:
         raise ValueError()
 
-    lon_data = np.linspace(lon_start + 0.5 * spatial_res, lon_end - 0.5 * spatial_res, width)
+    spatial_res_05 = 0.5 * spatial_res
+
+    lon_data = np.linspace(lon_start + spatial_res_05, lon_end - spatial_res_05, width)
     lon = xr.DataArray(lon_data, dims="lon")
     lon.attrs["units"] = "degrees_east"
 
     lat_data = np.linspace(lat_start + 0.5 * spatial_res, lat_end - 0.5 * spatial_res, height)
     lat = xr.DataArray(lat_data, dims="lat")
+    if inverse_lat:
+        lat = lat[::-1]
     lat.attrs["units"] = "degrees_north"
 
     time_data_2 = pd.date_range(start=time_start, periods=time_periods + 1, freq=time_freq).values
@@ -72,6 +78,8 @@ def new_cube(title="Test Cube",
         lat_bnds_data = np.zeros((height, 2), dtype=np.float64)
         lat_bnds_data[:, 0] = np.linspace(lat_start, lat_end - spatial_res, height)
         lat_bnds_data[:, 1] = np.linspace(lat_start + spatial_res, lat_end, height)
+        if inverse_lat:
+            lat_bnds_data = lat_bnds_data[::-1, ::-1]
         lat_bnds = xr.DataArray(lat_bnds_data, dims=("lat", "bnds"))
         lat_bnds.attrs["units"] = "degrees_north"
 
diff --git a/xcube/cli/cli.py b/xcube/cli/cli.py
index 16fc3600..d0c92a0e 100644
--- a/xcube/cli/cli.py
+++ b/xcube/cli/cli.py
@@ -3,6 +3,7 @@ import sys
 import click
 
 from xcube.cli.apply import apply
+from xcube.cli.extract import extract
 from xcube.cli.gen import gen
 from xcube.cli.grid import grid
 from xcube.cli.prune import prune
@@ -13,47 +14,6 @@ from xcube.util.cliutil import parse_cli_kwargs, new_cli_ctx_obj, handle_cli_exc
     cli_option_scheduler
 from xcube.version import version
 
-DEFAULT_TILE_SIZE = 512
-
-
-@click.command(name="extract")
-@click.argument('cube', metavar='<cube>')
-@click.argument('coords', metavar='<coords>')
-@click.option('--indexes', '-i', is_flag=True,
-              help="Include indexes in output.")
-@click.option('--output', '-o', metavar='<output>',
-              help="Output file.")
-# @click.option('--format', '-f', metavar='<format>', type=click.Choice(['csv', 'stdout']),
-#               help="Format of the output. If not given, guessed from <output>, otherwise <stdout> is used.")
-# @click.option('--params', '-p', metavar='<params>',
-#               help="Parameters specific for the output format."
-#                    " Comma-separated list of <key>=<value> pairs.")
-def extract(cube, coords, indexes=False, output=None,
-            # format=None, params=None
-            ):
-    """
-    Extract cube time series.
-    Extracts data from <cube> at points given by coordinates <coords> and writes the resulting
-    time series to <output>.
-    """
-    import pandas as pd
-
-    cube_path = cube
-    coords_path = coords
-    output_path = output
-    include_indexes = indexes
-
-    from xcube.api import open_dataset, get_cube_values_for_points
-
-    coords = pd.read_csv(coords_path, parse_dates=["time"], infer_datetime_format=True)
-    print(coords, [coords[c].values.dtype for c in coords])
-    with open_dataset(cube_path) as cube:
-        values = get_cube_values_for_points(cube, coords, include_indexes=include_indexes)
-        if output_path:
-            values.to_csv(output_path)
-        else:
-            print(values)
-
 
 # noinspection PyShadowingBuiltins
 @click.command(name="chunk")
@@ -102,6 +62,9 @@ def chunk(input, output, format=None, params=None, chunks=None):
         write_dataset(chunked_dataset, output_path=output, format_name=format_name, **write_kwargs)
 
 
+DEFAULT_TILE_SIZE = 512
+
+
 # noinspection PyShadowingBuiltins
 @click.command(name="level")
 @click.argument('input', metavar='<input>')
diff --git a/xcube/cli/extract.py b/xcube/cli/extract.py
new file mode 100644
index 00000000..80df88b0
--- /dev/null
+++ b/xcube/cli/extract.py
@@ -0,0 +1,76 @@
+import sys
+
+import click
+
+
+# noinspection PyShadowingBuiltins
+@click.command(name="extract")
+@click.argument('cube')
+@click.argument('points')
+@click.option('--output', '-o',
+              help="Output file. If omitted, output is written to stdout.")
+@click.option('--format', '-f', type=click.Choice(['csv', 'json', 'xlsx']),
+              help="Output format. Currently, only 'csv' is supported.", default='csv')
+@click.option('--coords', '-C', is_flag=True,
+              help="Include cube cell coordinates in output.")
+@click.option('--bounds', '-B', is_flag=True,
+              help="Include cube cell coordinate boundaries (if any) in output.")
+@click.option('--indexes', '-I', is_flag=True,
+              help="Include cube cell indexes in output.")
+@click.option('--refs', '-R', is_flag=True,
+              help="Include point values as reference in output.")
+def extract(cube,
+            points,
+            output=None,
+            format=None,
+            coords=False,
+            bounds=False,
+            indexes=False,
+            refs=False,
+):
+    """
+    Extract cube points.
+
+    Extracts data cells from CUBE at coordinates given in each POINTS record and writes the resulting values to given
+    output path and format.
+
+    <points> must be a CSV file that provides at least the columns "lon", "lat", and "time". The "lon" and "lat"
+    columns provide a point's location in decimal degrees. The "time" column provides a point's date or
+    date-time. Its format should preferably be ISO, but other formats may work as well.
+    """
+    if format != 'csv':
+        raise click.ClickException(f'Format {format!r} is not supported.')
+
+    import pandas as pd
+
+    cube_path = cube
+    points_path = points
+    output_path = output
+    include_coords = coords
+    include_bounds = bounds
+    include_indexes = indexes
+    include_refs = refs
+
+    from xcube.api.readwrite import open_dataset
+    from xcube.api.extract import get_cube_values_for_points, DEFAULT_INDEX_NAME_PATTERN, DEFAULT_REF_NAME_PATTERN
+
+    # We may make the following CLI options
+    index_name_pattern = DEFAULT_INDEX_NAME_PATTERN
+    ref_name_pattern = DEFAULT_REF_NAME_PATTERN
+    time_col_names = ["time"]
+
+    points = pd.read_csv(points_path, parse_dates=time_col_names, infer_datetime_format=True)
+    with open_dataset(cube_path) as cube:
+        values = get_cube_values_for_points(cube,
+                                            points,
+                                            include_coords=include_coords,
+                                            include_bounds=include_bounds,
+                                            include_indexes=include_indexes,
+                                            index_name_pattern=index_name_pattern,
+                                            include_refs=include_refs,
+                                            ref_name_pattern=ref_name_pattern).to_dataframe()
+        values.to_csv(output_path if output_path else sys.stdout,
+                      # We may make the following CLI options
+                      sep=',',
+                      date_format='%Y-%m-%dT%H:%M:%SZ',
+                      index=True)


Test patch:
diff --git a/test/api/test_extract.py b/test/api/test_extract.py
index c237dcc1..fc55d8ed 100644
--- a/test/api/test_extract.py
+++ b/test/api/test_extract.py
@@ -9,7 +9,7 @@ from xcube.api.new import new_cube
 
 
 # noinspection PyMethodMayBeStatic
-class ExtractPointsTest(unittest.TestCase):
+class GetCubeValuesForPointsTest(unittest.TestCase):
     def _new_test_cube(self):
         return new_cube(width=2000,
                         height=1000,
@@ -50,8 +50,13 @@ class ExtractPointsTest(unittest.TestCase):
         values = get_cube_values_for_points(cube, points, include_indexes=True)
         self.assertIsInstance(values, xr.Dataset)
         self.assertEqual(5, len(values.data_vars))
-        self.assertEqual(['precipitation', 'temperature',
-                          'time_index', 'lat_index', 'lon_index'], [c for c in values.data_vars])
+        self.assertEqual({
+            'precipitation',
+            'temperature',
+            'lat_index',
+            'lon_index',
+            'time_index'},
+            set(values.data_vars))
         np.testing.assert_array_almost_equal(np.array(expected_prec_values, dtype=np.float64),
                                              values["precipitation"].values)
         np.testing.assert_array_almost_equal(np.array(expected_temp_values, dtype=np.float64),
@@ -73,7 +78,7 @@ class ExtractPointsTest(unittest.TestCase):
         indexes = get_cube_point_indexes(cube, self._new_test_points())
         self.assertIsInstance(indexes, xr.Dataset)
         self.assertEqual(3, len(indexes.data_vars))
-        self.assertEqual(['time_index', 'lat_index', 'lon_index'], [c for c in indexes.data_vars])
+        self.assertEqual({'time_index', 'lat_index', 'lon_index'}, set(indexes.data_vars))
         np.testing.assert_array_almost_equal(np.array(expected_time_index, dtype=np.float64),
                                              indexes["time_index"].values)
         np.testing.assert_array_almost_equal(np.array(expected_lat_index, dtype=np.float64),
@@ -83,9 +88,9 @@ class ExtractPointsTest(unittest.TestCase):
 
 
 # noinspection PyMethodMayBeStatic
-class CoordIndexTest(unittest.TestCase):
+class GetDatasetIndexesTest(unittest.TestCase):
 
-    def test_get_coord_index_for_single_cell(self):
+    def test_get_dataset_indexes_for_single_cell(self):
         dataset = new_cube(width=360, height=180, drop_bounds=True)
         cell = dataset.isel(time=2, lat=20, lon=30)
         with self.assertRaises(ValueError) as cm:
@@ -93,57 +98,59 @@ class CoordIndexTest(unittest.TestCase):
         self.assertEqual("cannot determine cell boundaries for coordinate variable 'lon' of size 1",
                          f"{cm.exception}")
 
-    def test_get_coord_index_with_bounds(self):
+    def test_get_dataset_indexes_with_bounds(self):
         dataset = new_cube(width=360, height=180, drop_bounds=False)
-        self._assert_get_coord_index(dataset)
+        self._assert_get_dataset_indexes_works(dataset)
 
-    def test_get_coord_index_without_bounds(self):
+    def test_get_dataset_indexes_without_bounds(self):
         dataset = new_cube(width=360, height=180, drop_bounds=True)
-        self._assert_get_coord_index(dataset)
+        self._assert_get_dataset_indexes_works(dataset)
 
-    def test_get_coord_index_with_bounds_reverse_lat(self):
-        dataset = new_cube(width=360, height=180, drop_bounds=False)
-        lat = dataset.lat.values[::-1]
-        lat_bnds = dataset.lat_bnds.values[::-1, ::-1]
-        dataset.coords["lat"] = xr.DataArray(lat, dims=("lat",))
-        dataset.coords["lat_bnds"] = xr.DataArray(lat_bnds, dims=("lat", "bnds"))
-        self._assert_get_coord_index(dataset, reverse_lat=True)
+    def test_get_dataset_indexes_with_bounds_inverse_lat(self):
+        dataset = new_cube(width=360, height=180, inverse_lat=True, drop_bounds=False)
+        self._assert_get_dataset_indexes_works(dataset, inverse_lat=True)
 
-    def test_get_coord_index_without_bounds_reverse_lat(self):
-        dataset = new_cube(width=360, height=180, drop_bounds=True)
-        lat = dataset.lat.values[::-1]
-        dataset.coords["lat"] = xr.DataArray(lat, dims=("lat",))
-        self._assert_get_coord_index(dataset, reverse_lat=True)
+    def test_get_dataset_indexes_without_bounds_inverse_lat(self):
+        dataset = new_cube(width=360, height=180, inverse_lat=True, drop_bounds=True)
+        self._assert_get_dataset_indexes_works(dataset, inverse_lat=True)
+
+    def _assert_get_dataset_indexes_works(self, dataset, inverse_lat=False):
+        # - lon -------------------------------------
 
-    def _assert_get_coord_index(self, dataset, reverse_lat=False):
         lon_coords = np.array([-190, -180., -179, -10.4, 0., 10.4, 179., 180.0, 190])
         expected_lon_int64 = np.array([-1, 0, 1, 169, 180, 190, 359, 359, -1], dtype=np.int64)
         expected_lon_float64 = np.array([np.nan, 0., 1., 169.6, 180., 190.4, 359., 360., np.nan], dtype=np.float64)
 
-        lat_coords = np.array([-100, -90., -89, -10.4, 0., 10.4, 89., 90.0, 100])
-        expected_lat_int64 = np.array([-1, 0, 1, 79, 90, 100, 179, 179, -1], dtype=np.int64)
-        expected_lat_float64 = np.array([np.nan, 0., 1., 79.6, 90., 100.4, 179., 180., np.nan], dtype=np.float64)
-        if reverse_lat:
-            expected_lat_int64 = expected_lat_int64[::-1]
-            expected_lat_float64 = expected_lat_float64[::-1]
-
-        time_coords = np.array(
-            ["2010-01-03", "2010-01-02T23:15", "2010-05-15", "2010-01-01T12:00", "2010-01-05", "2009-12-07"],
-            dtype="datetime64[ns]")
-        expected_time_int64 = np.array([2, 1, -1, 0, 4, -1], dtype=np.int64)
-        expected_time_float64 = np.array([2., 1.96875, np.nan, 0.5, 4., np.nan], dtype=np.float64)
-
-        indexes = get_dataset_indexes(dataset, "lon", lon_coords, dtype=np.int32)
+        indexes = get_dataset_indexes(dataset, "lon", lon_coords, index_dtype=np.int32)
         np.testing.assert_array_equal(indexes, expected_lon_int64)
         indexes = get_dataset_indexes(dataset, "lon", lon_coords)
         np.testing.assert_array_almost_equal(indexes, expected_lon_float64)
 
-        indexes = get_dataset_indexes(dataset, "lat", lat_coords, dtype=np.int32)
+        # - lat -------------------------------------
+
+        lat_coords = np.array([-100, -90., -89, -10.4, 0., 10.4, 89., 90.0, 100])
+        if not inverse_lat:
+            expected_lat_int64 = np.array([-1, 0, 1, 79, 90, 100, 179, 179, -1], dtype=np.int64)
+            expected_lat_float64 = np.array([np.nan, 0., 1., 79.6, 90., 100.4, 179., 180., np.nan], dtype=np.float64)
+        else:
+            expected_lat_int64 = np.array([-1, 179, 179, 100, 90, 79, 1, 0, -1], dtype=np.int64)
+            expected_lat_float64 = np.array([np.nan, 180., 179., 100.4, 90., 79.6, 1., 0., np.nan],
+                                            dtype=np.float64)
+
+        indexes = get_dataset_indexes(dataset, "lat", lat_coords, index_dtype=np.int32)
         np.testing.assert_array_equal(indexes, expected_lat_int64)
         indexes = get_dataset_indexes(dataset, "lat", lat_coords)
         np.testing.assert_array_almost_equal(indexes, expected_lat_float64)
 
-        indexes = get_dataset_indexes(dataset, "time", time_coords, dtype=np.int32)
+        # - time ------------------------------------
+
+        time_coords = np.array(
+            ["2010-01-03", "2010-01-02T23:15", "2010-05-15", "2010-01-01T12:00", "2010-01-05", "2009-12-07"],
+            dtype="datetime64[ns]")
+        expected_time_int64 = np.array([2, 1, -1, 0, 4, -1], dtype=np.int64)
+        expected_time_float64 = np.array([2., 1.96875, np.nan, 0.5, 4., np.nan], dtype=np.float64)
+
+        indexes = get_dataset_indexes(dataset, "time", time_coords, index_dtype=np.int32)
         np.testing.assert_array_equal(indexes, expected_time_int64)
         indexes = get_dataset_indexes(dataset, "time", time_coords)
         np.testing.assert_array_almost_equal(indexes, expected_time_float64)
diff --git a/test/cli/extract-points/all-invalid-points.csv b/test/cli/extract-points/all-invalid-points.csv
new file mode 100644
index 00000000..010ce6a7
--- /dev/null
+++ b/test/cli/extract-points/all-invalid-points.csv
@@ -0,0 +1,7 @@
+id,lat,lon,time
+0,51.5,2.0,2017-01-16T08:00:00Z
+1,51.5,2.1,2017-01-16T08:00:00Z
+2,51.5,2.2,2017-01-16T08:00:00Z
+3,51.5,2.3,2017-02-01T12:00:00Z
+4,51.5,2.4,2017-02-01T12:00:00Z
+5,51.5,2.5,2017-02-01T12:00:00Z
diff --git a/test/cli/extract-points/all-valid-points.csv b/test/cli/extract-points/all-valid-points.csv
new file mode 100644
index 00000000..f0554a55
--- /dev/null
+++ b/test/cli/extract-points/all-valid-points.csv
@@ -0,0 +1,7 @@
+id,lat,lon,time
+0,51.5,2.0,2017-01-16T12:00:00Z
+1,51.5,2.1,2017-01-16T12:00:00Z
+2,51.5,2.2,2017-01-16T12:00:00Z
+3,51.5,2.3,2017-01-28T12:00:00Z
+4,51.5,2.4,2017-01-28T12:00:00Z
+5,51.5,2.5,2017-01-28T12:00:00Z
diff --git a/test/cli/extract-points/some-invalid-points.csv b/test/cli/extract-points/some-invalid-points.csv
new file mode 100644
index 00000000..03553772
--- /dev/null
+++ b/test/cli/extract-points/some-invalid-points.csv
@@ -0,0 +1,7 @@
+id,lat,lon,time
+0,51.5,2.0,2017-01-16T08:00:00Z
+1,51.5,2.1,2017-01-16T08:00:00Z
+2,51.5,2.2,2017-01-16T08:00:00Z
+3,51.5,2.3,2017-01-28T12:00:00Z
+4,51.5,2.4,2017-01-28T12:00:00Z
+5,51.5,2.5,2017-01-28T12:00:00Z
diff --git a/test/cli/test_extract.py b/test/cli/test_extract.py
new file mode 100644
index 00000000..d0c2cc39
--- /dev/null
+++ b/test/cli/test_extract.py
@@ -0,0 +1,99 @@
+import os
+import os.path
+from io import StringIO
+
+import pandas as pd
+
+from test.cli.test_cli import CliDataTest
+
+VAR_COLS = [
+    'c2rcc_flags',
+    'conc_chl',
+    'conc_tsm',
+    'kd489',
+    'quality_flags',
+]
+
+COORDS_COLS = [
+    'lon',
+    'lat',
+    'time',
+]
+
+BOUNDS_COLS = [
+    'lon_lower',
+    'lon_upper',
+    'lat_lower',
+    'lat_upper',
+    'time_lower',
+    'time_upper',
+]
+
+INDEX_COLS = [
+    'lon_index',
+    'lat_index',
+    'time_index',
+]
+
+REF_COLS = [
+    'id_ref',
+    'lon_ref',
+    'lat_ref',
+    'time_ref',
+]
+
+
+class ExtractCliTest(CliDataTest):
+
+    def test_help_option(self):
+        result = self.invoke_cli(['extract', '--help'])
+        self.assertEqual(0, result.exit_code)
+
+    def test_all_invalid_points(self):
+        self._do_test_points('all-invalid-points.csv')
+
+    def test_some_invalid_points(self):
+        self._do_test_points('some-invalid-points.csv')
+
+    def test_extract_points_with_refs(self):
+        self._do_test_points('all-valid-points.csv')
+
+    def _do_test_points(self, filename):
+        self._do_test_points_with_options(filename, options=[],
+                                          expected_length=6)
+        self._do_test_points_with_options(filename, options=['--coords'],
+                                          expected_length=6)
+        self._do_test_points_with_options(filename, options=['--bounds'],
+                                          expected_length=6)
+        self._do_test_points_with_options(filename, options=['--refs'],
+                                          expected_length=6)
+        self._do_test_points_with_options(filename, options=['--indexes'],
+                                          expected_length=6)
+        self._do_test_points_with_options(filename, options=['--coords', '--bounds', '--refs', '--indexes'],
+                                          expected_length=6)
+
+    def _do_test_points_with_options(self,
+                                     filename,
+                                     options=None,
+                                     expected_length=0):
+        base_dir = os.path.dirname(__file__)
+        args = ['extract']
+        if options:
+            args += options
+        args += [os.path.join(base_dir, '..', '..', 'xcube', 'webapi', 'res', 'demo', 'cube.nc'),
+                 os.path.join(base_dir, 'extract-points', filename)]
+        result = self.invoke_cli(args)
+        self.assertEqual(0, result.exit_code)
+        self.assertNotEqual('', result.stdout)
+        df = pd.read_csv(StringIO(result.stdout), index_col='idx')
+        self.assertEqual(expected_length, len(df))
+        expected_columns = set(VAR_COLS)
+        if '--coords' in options:
+            expected_columns.update(COORDS_COLS)
+        if '--bounds' in options:
+            expected_columns.update(BOUNDS_COLS)
+        if '--refs' in options:
+            expected_columns.update(REF_COLS)
+        if '--indexes' in options:
+            expected_columns.update(INDEX_COLS)
+        self.assertEqual(expected_columns, set(df.columns))