Repository: NCAS-CMS/cfdm

Issue:
Create indices for Field constructs from metadata values
@JonathanGregory suggested that it would be natural in cfdm to be able to subspace field constructs based on metadata constructs values, rather than just by indices. An excellent idea.

For instance, with
```python
>>> import cfdm
>>> f = cfdm.example_field(0)
>>> print(f)
Field: specific_humidity (ncvar%q)
----------------------------------
Data            : specific_humidity(latitude(5), longitude(8)) 1
Cell methods    : area: mean
Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
                : longitude(8) = [22.5, ..., 337.5] degrees_east
                : time(1) = [2019-01-01 00:00:00]
```

we currently can only create a new field for only longitude=337.5 by:
```python
>>> print(f[:, 7])
Field: specific_humidity (ncvar%q)
----------------------------------
Data            : specific_humidity(latitude(5), longitude(1)) 1
Cell methods    : area: mean
Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
                : longitude(1) = [337.5] degrees_east
                : time(1) = [2019-01-01 00:00:00]
```
However, this requires knowledge of **a)** which position the longitude dimension is and **b)** where along the longitude axis `337.5` is.

It seems a natural use of CF metadata to be able to do this _without_ have to know a priori the answers to **a)** and **b)**.

cf-python does this with its `Field.indices`  and `Field.subspace` methods.

These methods are fully generalised, rely on the `cf.Query` class, and take cyclicity into account.

I don't propose that cfdm goes this far (certainly not yet!), so I suggest that we for cfdm we just introduce a `cfdm.Field.indices` method that behave much like its cf-python counterpart, but is restricted to
* Only 1-d metadata constructs (no domain axis constructs nor N-d constructs)
* Only metadata construct data values can be provided (as oposed to axis indices)
* No fancy "modes" (full, envelope, compress)  - cfdm doesn't have the machinery to back them up

The result would be:
```python
>>> ind = f.indices(longitude=337.5)
>>> print(ind)
(slice(None, None, None), array([False, False, False, False, False, False, False,  True]))
>>> print(f[ind])
Field: specific_humidity (ncvar%q)
----------------------------------
Data            : specific_humidity(latitude(5), longitude(1)) 1
Cell methods    : area: mean
Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
                : longitude(1) = [337.5] degrees_east
                : time(1) = [2019-01-01 00:00:00]
```

The default position is to not burden cfdm with too much bloat (relative to cf-python), and cfdm started off very minimal. Since then there has been a trickle for back-ported functionality from cf-python to cfdm that has become obviously needed by people whilst actully using the library. There is no problem in doing this so long as each case is assessed on its own merits. 

Discussion:

Patch:
diff --git a/Changelog.rst b/Changelog.rst
index 95e4edffa..8732b1f6a 100644
--- a/Changelog.rst
+++ b/Changelog.rst
@@ -12,6 +12,7 @@ Version 1.9.1.0
 * New method: `cfdm.Field.domain_ancillary`
 * New method: `cfdm.Field.domain_axis`
 * New method: `cfdm.Field.field_ancillary`
+* New method: `cfdm.Field.indices`
 * New attribute: `cfdm.Field.array`
 * New attribute: `cfdm.Field.datetime_array`
 * New construct retrieval API methods
diff --git a/README.md b/README.md
index b9ccc6754..575448cc0 100644
--- a/README.md
+++ b/README.md
@@ -81,7 +81,8 @@ The ``cfdm`` package can:
 * inspect field and domain constructs,
 * test whether two constructs are the same,
 * modify field and domain construct metadata and data,
-* create subspaces of field and domain constructs,
+* create subspaces of field and domain constructs, from indices or
+  metadata values,
 * incorporate, and create, metadata stored in external files, and
 * read, write, and create data that have been compressed by convention
   (i.e. ragged or gathered arrays, or coordinate arrays compressed by
diff --git a/cfdm/field.py b/cfdm/field.py
index f126c8078..4115643ec 100644
--- a/cfdm/field.py
+++ b/cfdm/field.py
@@ -1,5 +1,7 @@
 import logging
 
+import numpy as np
+
 from . import Constructs, Count, Domain, Index, List, core, mixin
 from .constants import masked as cfdm_masked
 from .data import (
@@ -1882,6 +1884,194 @@ class Field(
     #
     #        return False
 
+    def indices(self, **kwargs):
+        """Create indices that define a subspace of the field construct.
+
+        The subspace is defined by defining indices based on the
+        positions of the given data values of 1-d metadata constructs.
+
+        The returned tuple of indices may be used to created a subspace by
+        indexing the original field construct with them.
+
+        Metadata constructs and the conditions on their data are defined
+        by keyword parameters.
+
+        * Any domain axes that have not been identified remain unchanged.
+
+        * Multiple domain axes may be subspaced simultaneously, and it
+          doesn't matter which order they are specified in.
+
+        * Subspace criteria may be provided for size 1 domain axes that
+          are not spanned by the field construct's data.
+
+        .. versionadded:: 1.9.1.0
+
+        .. seealso:: `__getitem__`, `__setitem__`
+
+        :Parameters:
+
+            kwargs: *optional*
+                Each keyword parameter specifies a value or values
+                whose positions in the 1-d metadata construct's data,
+                identified by the parameter name, define the indices
+                for that dimension.
+
+        :Returns:
+
+            `tuple`
+                The indices meeting the conditions.
+
+        **Examples**
+
+        >>> f = cfdm.example_field(0)
+        >>> print(f)
+        Field: specific_humidity (ncvar%q)
+        ----------------------------------
+        Data            : specific_humidity(latitude(5), longitude(8)) 1
+        Cell methods    : area: mean
+        Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
+                        : longitude(8) = [22.5, ..., 337.5] degrees_east
+                        : time(1) = [2019-01-01 00:00:00]
+
+        >>> indices = f.indices(longitude=112.5)
+        >>> indices
+        (slice(None, None, None),
+         array([False, False,  True, False, False, False, False, False]))
+        >>> print(f[indices])
+        Field: specific_humidity (ncvar%q)
+        ----------------------------------
+        Data            : specific_humidity(latitude(5), longitude(1)) 1
+        Cell methods    : area: mean
+        Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
+                        : longitude(1) = [112.5] degrees_east
+                        : time(1) = [2019-01-01 00:00:00]
+
+        >>> indices = f.indices(longitude=112.5, latitude=[-45, 75])
+        >>> indices
+        (array([False,  True, False, False,  True]),
+         array([False, False,  True, False, False, False, False, False]))
+        >>> print(f[indices])
+        Field: specific_humidity (ncvar%q)
+        ----------------------------------
+        Data            : specific_humidity(latitude(2), longitude(1)) 1
+        Cell methods    : area: mean
+        Dimension coords: latitude(2) = [-45.0, 75.0] degrees_north
+                        : longitude(1) = [112.5] degrees_east
+                        : time(1) = [2019-01-01 00:00:00]
+
+        >>> indices = f.indices(time=31)
+        >>> indices
+        (slice(None, None, None), slice(None, None, None))
+        >>> print(f[indices])
+        Field: specific_humidity (ncvar%q)
+        ----------------------------------
+        Data            : specific_humidity(latitude(5), longitude(8)) 1
+        Cell methods    : area: mean
+        Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
+                        : longitude(8) = [22.5, ..., 337.5] degrees_east
+                        : time(1) = [2019-01-01 00:00:00]
+
+        """
+        # Initialise indices dictionary
+        indices = {axis: slice(None) for axis in self.domain_axes(todict=True)}
+
+        parsed = {}
+        unique_axes = set()
+        n_axes = 0
+
+        for identity, value in kwargs.items():
+            key, construct = self.construct(
+                identity,
+                filter_by_data=True,
+                item=True,
+                default=(None, None),
+                filter_by_naxes=(1,),
+            )
+            if construct is not None:
+                axes = self.get_data_axes(key)
+            else:
+                da_key = self.domain_axis(identity, key=True, default=None)
+                if da_key is not None:
+                    axes = (da_key,)
+                    key = None
+                    construct = None
+                else:
+                    raise ValueError(
+                        f"Can't find indices. Ambiguous axis or axes "
+                        f"defined by {identity!r}"
+                    )
+
+            if axes in parsed:
+                # The axes are the same as an existing key
+                parsed[axes].append((axes, key, construct, value, identity))
+            else:
+                new_key = True
+                y = set(axes)
+                for x in parsed:
+                    if set(x) == set(y):
+                        # The axes are the same but in a different
+                        # order, so we don't need a new key.
+                        parsed[x].append(
+                            (axes, key, construct, value, identity)
+                        )
+                        new_key = False
+                        break
+
+                if new_key:
+                    # The axes, taken in any order, are not the same
+                    # as any keys, so create an new key.
+                    n_axes += len(axes)
+                    parsed[axes] = [(axes, key, construct, value, identity)]
+
+            unique_axes.update(axes)
+
+        if len(unique_axes) < n_axes:
+            raise ValueError(
+                "Can't find indices: Multiple constructs with incompatible "
+                "domain axes"
+            )
+
+        for canonical_axes, axes_key_construct_value in parsed.items():
+            axes, keys, constructs, points, identities = list(
+                zip(*axes_key_construct_value)
+            )
+
+            n_items = len(constructs)
+            n_axes = len(canonical_axes)
+
+            if n_axes != 1:
+                raise ValueError("TODO")
+
+            if n_items > n_axes:
+                raise ValueError(
+                    "Can't specify the same axis more than once. Got: "
+                    f"{identities}"
+                )
+
+            axis = axes[0][0]
+            item = constructs[0]
+            value = points[0]
+
+            if item is None:
+                raise ValueError(
+                    "Can only specify 1-d metadata constructs from which "
+                    f"to create indices. Got: {identities[0]!r}"
+                )
+
+            index = np.isin(item, np.asanyarray(value).astype(item.dtype))
+            if np.ma.is_masked(index):
+                index = index.filled(False)
+
+            if not index.any():
+                raise ValueError(
+                    f"{value!r} does not match any {item!r} data values"
+                )
+
+            indices[axis] = index
+
+        # Return indices tuple
+        return tuple([indices[axis] for axis in self.get_data_axes()])
+
     @_inplace_enabled(default=False)
     def insert_dimension(self, axis, position=0, inplace=False):
         """Expand the shape of the data array.
diff --git a/docs/Makefile b/docs/Makefile
index 08e966dea..1abc4a23a 100644
--- a/docs/Makefile
+++ b/docs/Makefile
@@ -11,7 +11,7 @@ PAPER         =
 PAPEROPT_a4     = -D latex_paper_size=a4
 PAPEROPT_letter = -D latex_paper_size=letter
 #ALLSPHINXOPTS   = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source
-ALLSPHINXOPTS   = -j 3 source
+ALLSPHINXOPTS   = -j 2 source
 #ALLSPHINXOPTS   = source
 BUILDDIR        = $(filter-out $@, $(MAKECMDGOALS))
 
diff --git a/docs/source/class/cfdm.Field.rst b/docs/source/class/cfdm.Field.rst
index d2aceccf5..90de5c914 100644
--- a/docs/source/class/cfdm.Field.rst
+++ b/docs/source/class/cfdm.Field.rst
@@ -194,6 +194,7 @@ Miscellaneous
    ~cfdm.Field.convert
    ~cfdm.Field.has_bounds
    ~cfdm.Field.has_geometry
+   ~cfdm.Field.indices
    ~cfdm.Field.uncompress
    ~cfdm.Field.get_filenames
    ~cfdm.Field.to_memory
diff --git a/docs/source/tutorial.rst b/docs/source/tutorial.rst
index 934f482ea..5811cddf8 100644
--- a/docs/source/tutorial.rst
+++ b/docs/source/tutorial.rst
@@ -1164,6 +1164,63 @@ that apply to the `Data` class.
                    : latitude(5) = [75.0, ..., -75.0] degrees_north
                    : longitude(1) = [22.5] degrees_east
 
+.. _Subspacing-by-metadata:
+
+Subspacing by metadata
+^^^^^^^^^^^^^^^^^^^^^^
+
+Indices that correspond to the locations of particular metadata
+construct values are easily created with the `~Field.indices` method
+of the field construct.
+
+Metadata constructs and the conditions on their data are defined by
+keyword parameters to the `~Field.indices` method. A keyword name is
+:ref:`an identity <Metadata-construct-properties>` of a 1-d metadata
+construct, and its values are one or more construct array
+values. Indices are created that select every location for which the
+metadata construct's data equals any of the given values.
+
+	
+.. code-block:: python
+   :caption: *Create a new field constructs based on metadata-defined
+             subspaces.*
+	     
+   >>> print(q)
+   Field: specific_humidity (ncvar%q)
+   ----------------------------------
+   Data            : specific_humidity(latitude(5), longitude(8)) 1
+   Cell methods    : area: mean
+   Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
+                   : longitude(8) = [22.5, ..., 337.5] degrees_east
+                   : time(1) = [2019-01-01T00:00:00Z]
+   >>> print(q.construct('longitude').data.array)
+   [ 22.5  67.5 112.5 157.5 202.5 247.5 292.5 337.5]
+   >>> ind = q.indices(longitude=[112.5, 67.5])
+   >>> print(ind)
+   (slice(None, None, None),
+    array([False,  True,  True, False, False, False, False, False]))
+   >>> print(q[ind])
+   Field: specific_humidity (ncvar%q)
+   ----------------------------------
+   Data            : specific_humidity(latitude(5), longitude(2)) 1
+   Cell methods    : area: mean
+   Dimension coords: latitude(5) = [-75.0, ..., 75.0] degrees_north
+                   : longitude(2) = [67.5, 112.5] degrees_east
+                   : time(1) = [2019-01-01 00:00:00]
+   >>> print(q[q.indices(longitude=[112.5, 67.5], latitude=75)])
+   Field: specific_humidity (ncvar%q)
+   ----------------------------------
+   Data            : specific_humidity(latitude(1), longitude(2)) 1
+   Cell methods    : area: mean
+   Dimension coords: latitude(1) = [75.0] degrees_north
+	           : longitude(2) = [67.5, 112.5] degrees_east
+                   : time(1) = [2019-01-01 00:00:00]
+
+Any domain axes that have not been identified remain unchanged.
+
+Multiple domain axes may be subspaced simultaneously, and it doesn't
+matter which order they are specified in the `~Field.indices` call.
+				   
 ----
 
 .. _Selecting-metadata-constructs:
@@ -1599,7 +1656,7 @@ returning a default value or raising a customised exception:
 .. _Metadata-construct-properties:
 
 Metadata construct properties
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Metadata constructs share the :ref:`same API as the field construct
 <Properties>` for accessing their properties:


Test patch:
diff --git a/cfdm/test/test_Field.py b/cfdm/test/test_Field.py
index 993fc41c4..95df66640 100644
--- a/cfdm/test/test_Field.py
+++ b/cfdm/test/test_Field.py
@@ -6,7 +6,7 @@ import re
 import tempfile
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -130,7 +130,7 @@ class FieldTest(unittest.TestCase):
                 indices = list(indices)
                 for axis, i in enumerate(indices):
                     if isinstance(i, list):
-                        e = numpy.take(e, indices=i, axis=axis)
+                        e = np.take(e, indices=i, axis=axis)
                         indices[axis] = slice(None)
 
                 e = e[tuple(indices)]
@@ -177,7 +177,7 @@ class FieldTest(unittest.TestCase):
     #            array = f[indices].data.array
     #            self.assertTrue((array == -1).all())
     #
-    #            values, counts = numpy.unique(f.data.array, return_counts=True)
+    #            values, counts = np.unique(f.data.array, return_counts=True)
     #            self.assertEqual(counts[0], array.size)
 
     def test_Field_get_filenames(self):
@@ -790,6 +790,42 @@ class FieldTest(unittest.TestCase):
         with self.assertRaises(ValueError):
             f.domain_axis("qwerty")
 
+    def test_Field_indices(self):
+        """Test that Field.indices."""
+        f = cfdm.example_field(0)
+
+        g = f[f.indices(longitude=112.5)]
+        self.assertEqual(g.shape, (5, 1))
+        x = g.dimension_coordinate("longitude").data.array
+        self.assertTrue((x == 112.5).all())
+
+        g = f[f.indices(longitude=112.5, latitude=[-45, 75])]
+        self.assertEqual(g.shape, (2, 1))
+        x = g.dimension_coordinate("longitude").data.array
+        y = g.dimension_coordinate("latitude").data.array
+        self.assertTrue((x == 112.5).all())
+        self.assertTrue((y == [-45, 75]).all())
+
+        g = f[f.indices(time=31)]
+        self.assertTrue(g.equals(f))
+
+        g = f[f.indices(time=np.array([31, 9999]))]
+        self.assertTrue(g.equals(f))
+
+        with self.assertRaises(ValueError):
+            f.indices(bad_name=23)
+
+        with self.assertRaises(ValueError):
+            f.indices(longitude=-999)
+
+        with self.assertRaises(ValueError):
+            f.indices(longitude="bad_value_type")
+
+        # Test for same axis specified twice
+        key = f.construct("longitude", key=True)
+        with self.assertRaises(ValueError):
+            f.indices(**{"longitude": 112.5, key: 22.5})
+
 
 if __name__ == "__main__":
     print("Run date:", datetime.datetime.now())