diff --git a/cfdm/cellmethod.py b/cfdm/cellmethod.py
index f0323de21..bb782e60b 100644
--- a/cfdm/cellmethod.py
+++ b/cfdm/cellmethod.py
@@ -1,7 +1,7 @@
 import logging
 from copy import deepcopy
 
-import numpy
+import numpy as np
 
 from . import core, mixin
 from .data import Data
@@ -529,7 +529,7 @@ def sorted(self, indices=None):
             return new
 
         if indices is None:
-            indices = numpy.argsort(axes)
+            indices = np.argsort(axes)
         elif len(indices) != len(axes):
             raise ValueError(
                 f"Can't sort cell method axes. The given indices ({indices}) "
diff --git a/cfdm/constants.py b/cfdm/constants.py
index c8c940d7e..c04a0d202 100644
--- a/cfdm/constants.py
+++ b/cfdm/constants.py
@@ -2,7 +2,7 @@
 import sys
 from enum import Enum
 
-import numpy
+import numpy as np
 
 """A dictionary of useful constants.
 
@@ -57,4 +57,4 @@ class ValidLogLevels(Enum):
 >>> f[...] = cfdm.masked
 
 """
-masked = numpy.ma.masked
+masked = np.ma.masked
diff --git a/cfdm/core/__init__.py b/cfdm/core/__init__.py
index d89e4aba2..1eee20aa3 100644
--- a/cfdm/core/__init__.py
+++ b/cfdm/core/__init__.py
@@ -31,7 +31,7 @@
     raise ImportError(_error0 + str(error1))
 
 try:
-    import numpy
+    import numpy as np
 except ImportError as error1:
     raise ImportError(_error0 + str(error1))
 
@@ -61,10 +61,10 @@
 
 # Check the version of numpy
 _minimum_vn = "1.15"
-if Version(numpy.__version__) < Version(_minimum_vn):
+if Version(np.__version__) < Version(_minimum_vn):
     raise ValueError(
         f"Bad numpy version: cfdm.core requires numpy>={_minimum_vn}. "
-        f"Got {numpy.__version__} at {numpy.__file__}"
+        f"Got {np.__version__} at {np.__file__}"
     )
 
 from .constructs import Constructs
diff --git a/cfdm/core/data/data.py b/cfdm/core/data/data.py
index b578662fd..a08b5d868 100644
--- a/cfdm/core/data/data.py
+++ b/cfdm/core/data/data.py
@@ -1,4 +1,4 @@
-import numpy
+import numpy as np
 
 from .. import abstract
 from .abstract import Array
@@ -167,7 +167,7 @@ def array(self):
         array = self._get_Array().array
 
         # Set the numpy array fill value
-        if numpy.ma.isMA(array):
+        if np.ma.isMA(array):
             array.set_fill_value(self.get_fill_value(None))
 
         return array
@@ -788,8 +788,8 @@ def _set_Array(self, array, copy=True):
 
         """
         if not isinstance(array, Array):
-            if not isinstance(array, numpy.ndarray):
-                array = numpy.asanyarray(array)
+            if not isinstance(array, np.ndarray):
+                array = np.asanyarray(array)
 
             array = NumpyArray(array)
 
diff --git a/cfdm/core/data/numpyarray.py b/cfdm/core/data/numpyarray.py
index 42c3e1c19..c74160070 100644
--- a/cfdm/core/data/numpyarray.py
+++ b/cfdm/core/data/numpyarray.py
@@ -1,4 +1,4 @@
-import numpy
+import numpy as np
 
 from . import abstract
 
@@ -101,10 +101,10 @@ def array(self):
         """
         array = self._get_component("array")
 
-        if not array.ndim and numpy.ma.isMA(array):
-            # This is because numpy.ma.copy doesn't work for
+        if not array.ndim and np.ma.isMA(array):
+            # This is because np.ma.copy doesn't work for
             # scalar arrays (at the moment, at least)
-            ma_array = numpy.ma.empty((), dtype=array.dtype)
+            ma_array = np.ma.empty((), dtype=array.dtype)
             ma_array[...] = array
             array = ma_array
         else:
diff --git a/cfdm/core/dimensioncoordinate.py b/cfdm/core/dimensioncoordinate.py
index 536c96890..cb0d41055 100644
--- a/cfdm/core/dimensioncoordinate.py
+++ b/cfdm/core/dimensioncoordinate.py
@@ -1,4 +1,4 @@
-import numpy
+import numpy as np
 
 from . import abstract
 
@@ -104,7 +104,7 @@ def set_data(self, data, copy=True, inplace=True):
         None
 
         """
-        if numpy.ndim(data) != 1:
+        if np.ndim(data) != 1:
             raise ValueError(
                 "Dimension coordinate construct must have 1-dimensional data. "
                 f"Got {data!r}"
diff --git a/cfdm/core/field.py b/cfdm/core/field.py
index 9a660d8c3..cc9210480 100644
--- a/cfdm/core/field.py
+++ b/cfdm/core/field.py
@@ -1,4 +1,4 @@
-import numpy
+import numpy as np
 
 from . import Constructs, Domain, abstract, mixin
 
@@ -531,9 +531,9 @@ def set_data(self, data, axes=None, copy=True, inplace=True):
         if axes is None:
             existing_axes = f.get_data_axes(default=None)
             if existing_axes is not None:
-                f.set_data_axes(axes=existing_axes, _shape=numpy.shape(data))
+                f.set_data_axes(axes=existing_axes, _shape=np.shape(data))
         else:
-            f.set_data_axes(axes=axes, _shape=numpy.shape(data))
+            f.set_data_axes(axes=axes, _shape=np.shape(data))
 
         super(Field, f).set_data(data, copy=copy, inplace=True)
 
diff --git a/cfdm/core/functions.py b/cfdm/core/functions.py
index e06c7af9c..4038aa2ec 100644
--- a/cfdm/core/functions.py
+++ b/cfdm/core/functions.py
@@ -4,7 +4,7 @@
 from pickle import dumps, loads
 
 import netCDF4
-import numpy
+import numpy as np
 
 from . import __cf_version__, __file__, __version__
 
@@ -60,7 +60,7 @@ def environment(display=True, paths=True):
         "netcdf library": (netCDF4.__netcdf4libversion__, ""),
         "Python": (platform.python_version(), sys.executable),
         "netCDF4": (netCDF4.__version__, os.path.abspath(netCDF4.__file__)),
-        "numpy": (numpy.__version__, os.path.abspath(numpy.__file__)),
+        "numpy": (np.__version__, os.path.abspath(np.__file__)),
         "cfdm.core": (__version__, os.path.abspath(__file__)),
     }
     string = "{0}: {1!s}"
diff --git a/cfdm/data/netcdfarray.py b/cfdm/data/netcdfarray.py
index af9a5193e..247c2747c 100644
--- a/cfdm/data/netcdfarray.py
+++ b/cfdm/data/netcdfarray.py
@@ -1,5 +1,5 @@
 import netCDF4
-import numpy
+import numpy as np
 
 from . import abstract
 from .mixin import FileArrayMixin
@@ -239,7 +239,7 @@ def __getitem__(self, indices):
             # A netCDF string type scalar variable comes out as Python
             # str object, so convert it to a numpy array.
             # --------------------------------------------------------
-            array = numpy.array(array, dtype=f"U{len(array)}")
+            array = np.array(array, dtype=f"U{len(array)}")
 
         if not self.ndim:
             # Hmm netCDF4 has a thing for making scalar size 1, 1d
@@ -260,9 +260,9 @@ def __getitem__(self, indices):
 
             array = netCDF4.chartostring(array)
             shape = array.shape
-            array = numpy.array([x.rstrip() for x in array.flat], dtype="U")
-            array = numpy.reshape(array, shape)
-            array = numpy.ma.masked_where(array == "", array)
+            array = np.array([x.rstrip() for x in array.flat], dtype="U")
+            array = np.reshape(array, shape)
+            array = np.ma.masked_where(array == "", array)
         elif not string_type and kind == "O":
             # --------------------------------------------------------
             # A netCDF string type N-d (N>=1) variable comes out as a
@@ -273,7 +273,7 @@ def __getitem__(self, indices):
             # --------------------------------------------------------
             # netCDF4 does not auto-mask VLEN variable, so do it here.
             # --------------------------------------------------------
-            array = numpy.ma.where(array == "", numpy.ma.masked, array)
+            array = np.ma.where(array == "", np.ma.masked, array)
 
         return array
 
diff --git a/cfdm/mixin/container.py b/cfdm/mixin/container.py
index 0709914f3..2a4577652 100644
--- a/cfdm/mixin/container.py
+++ b/cfdm/mixin/container.py
@@ -1,6 +1,6 @@
 import logging
 
-import numpy
+import numpy as np
 
 from ..decorators import _manage_log_level_via_verbosity
 from ..docstring import _docstring_substitution_definitions
@@ -143,19 +143,19 @@ def _equals(
             # --------------------------------------------------------
             return eq(x, **kwargs)
 
-        x = numpy.asanyarray(x)
-        y = numpy.asanyarray(y)
+        x = np.asanyarray(x)
+        y = np.asanyarray(y)
         if x.shape != y.shape:
             return False
 
         #        # ------------------------------------------------------------
         #        # Cast x and y as numpy arrays
         #        # ------------------------------------------------------------
-        #        if not isinstance(x, numpy.ndarray):
-        #            x = numpy.asanyarray(x)
+        #        if not isinstance(x, np.ndarray):
+        #            x = np.asanyarray(x)
         #
-        #        if not isinstance(y, numpy.ndarray):
-        #            y = numpy.asanyarray(y)
+        #        if not isinstance(y, np.ndarray):
+        #            y = np.asanyarray(y)
 
         # THIS IS WHERE SOME NUMPY FUTURE WARNINGS ARE COMING FROM
 
@@ -166,14 +166,14 @@ def _equals(
             ):
                 return False
 
-        x_is_masked = numpy.ma.isMA(x)
-        y_is_masked = numpy.ma.isMA(y)
+        x_is_masked = np.ma.isMA(x)
+        y_is_masked = np.ma.isMA(y)
 
         if not (x_is_masked or y_is_masked):
             try:
-                return bool(numpy.allclose(x, y, rtol=rtol, atol=atol))
+                return bool(np.allclose(x, y, rtol=rtol, atol=atol))
             except (IndexError, NotImplementedError, TypeError):
-                return bool(numpy.all(x == y))
+                return bool(np.all(x == y))
         else:
             if x_is_masked and y_is_masked:
                 if (x.mask != y.mask).any():
@@ -184,10 +184,10 @@ def _equals(
                 return False
 
             try:
-                return bool(numpy.ma.allclose(x, y, rtol=rtol, atol=atol))
+                return bool(np.ma.allclose(x, y, rtol=rtol, atol=atol))
             except (IndexError, NotImplementedError, TypeError):
-                out = numpy.ma.all(x == y)
-                if out is numpy.ma.masked:
+                out = np.ma.all(x == y)
+                if out is np.ma.masked:
                     return True
                 else:
                     return bool(out)
diff --git a/cfdm/read_write/netcdf/netcdfread.py b/cfdm/read_write/netcdf/netcdfread.py
index 4633b442b..fe43e867c 100644
--- a/cfdm/read_write/netcdf/netcdfread.py
+++ b/cfdm/read_write/netcdf/netcdfread.py
@@ -16,7 +16,6 @@
 
 import netCDF4
 import netcdf_flattener
-import numpy
 import numpy as np
 from packaging.version import Version
 
diff --git a/cfdm/read_write/netcdf/netcdfwrite.py b/cfdm/read_write/netcdf/netcdfwrite.py
index d8fdd0b56..81190d6ae 100644
--- a/cfdm/read_write/netcdf/netcdfwrite.py
+++ b/cfdm/read_write/netcdf/netcdfwrite.py
@@ -4,7 +4,6 @@
 import re
 
 import netCDF4
-import numpy
 import numpy as np
 from packaging.version import Version
 
@@ -186,7 +185,7 @@ def _numpy_compressed(self, array):
         <type 'numpy.ndarray'>
 
         """
-        if numpy.ma.isMA(array):
+        if np.ma.isMA(array):
             return array.compressed()
 
         return array.flatten()
@@ -237,9 +236,7 @@ def _write_attributes(self, parent, ncvar, extra=None, omit=()):
             data = self.implementation.get_data(parent, None)
             if data is not None:
                 dtype = g["datatype"].get(data.dtype, data.dtype)
-                netcdf_attrs[attr] = numpy.array(
-                    netcdf_attrs[attr], dtype=dtype
-                )
+                netcdf_attrs[attr] = np.array(netcdf_attrs[attr], dtype=dtype)
 
         skip_set_fill_value = False
         if g["post_dry_run"] and parent is not None:
@@ -310,20 +307,20 @@ def _character_array(self, array):
         original_shape = array.shape
         original_size = array.size
 
-        masked = numpy.ma.isMA(array)
+        masked = np.ma.isMA(array)
         if masked:
             fill_value = array.fill_value
-            array = numpy.ma.filled(array, fill_value="")
+            array = np.ma.filled(array, fill_value="")
 
         if array.dtype.kind == "U":
             array = array.astype("S")
 
-        array = numpy.array(tuple(array.tobytes().decode("ascii")), dtype="S1")
+        array = np.array(tuple(array.tobytes().decode("ascii")), dtype="S1")
 
         array.resize(original_shape + (array.size // original_size,))
 
         if masked:
-            array = numpy.ma.masked_where(array == b"", array)
+            array = np.ma.masked_where(array == b"", array)
             array.set_fill_value(fill_value)
 
         if array.dtype.kind != "S":
@@ -371,7 +368,7 @@ def _datatype(self, variable):
         """
         g = self.write_vars
 
-        if not isinstance(variable, numpy.ndarray):
+        if not isinstance(variable, np.ndarray):
             data = self.implementation.get_data(variable, None)
             if data is None:
                 return "S1"
@@ -1854,7 +1851,7 @@ def _write_part_node_count(self, f, coord, bounds, encodings):
         array = self.implementation.get_array(
             self.implementation.get_data(bounds)
         )
-        array = numpy.trim_zeros(numpy.ma.count(array, axis=2).flatten())
+        array = np.trim_zeros(np.ma.count(array, axis=2).flatten())
         array = self._int32(array)
 
         data = self.implementation.initialise_Data(array=array, copy=False)
@@ -2541,10 +2538,10 @@ def _write_grid_mapping(self, f, ref, multiple_grid_mappings):
                     del parameters[term]
                     continue
 
-                if numpy.size(value) == 1:
-                    value = numpy.array(value, copy=False).item()
+                if np.size(value) == 1:
+                    value = np.array(value, copy=False).item()
                 else:
-                    value = numpy.array(value, copy=False).tolist()
+                    value = np.array(value, copy=False).tolist()
 
                 parameters[term] = value
 
@@ -2903,7 +2900,7 @@ def _transform_strings(self, data, ncdimensions):
             # is str, so this conversion does not happen.
             # --------------------------------------------------------
             array = self.implementation.get_array(data)
-            #            if numpy.ma.is_masked(array):
+            #            if np.ma.is_masked(array):
             #                array = array.compressed()
             #            else:
             #                array = array.flatten()
@@ -2986,11 +2983,11 @@ def _write_data(
         # Check that the array doesn't contain any elements
         # which are equal to any of the missing data values
         if unset_values:
-            # if numpy.ma.is_masked(array):
+            # if np.ma.is_masked(array):
             #     temp_array = array.compressed()
             # else:
             #     temp_array = array
-            if numpy.intersect1d(
+            if np.intersect1d(
                 unset_values, self._numpy_compressed(array)
             ).size:
                 raise ValueError(
@@ -3001,7 +2998,7 @@ def _write_data(
         if (
             g["fmt"] == "NETCDF4"
             and array.dtype.kind in "SU"
-            and numpy.ma.isMA(array)
+            and np.ma.isMA(array)
         ):
             # VLEN variables can not be assigned to by masked arrays
             # https://github.com/Unidata/netcdf4-python/pull/465
@@ -3026,7 +3023,7 @@ def _check_valid(self, cfvar, array, attributes):
 
             cfvar: Construct
 
-            array: `numpy.ndarray`
+            array: `np.ndarray`
 
             attributes: `dict`
 
@@ -5076,8 +5073,8 @@ def _file_io_iteration(
         # converted to 64-bit floats.
         # ------------------------------------------------------------
         dtype_conversions = {
-            numpy.dtype(bool): numpy.dtype("int32"),
-            numpy.dtype(object): numpy.dtype(float),
+            np.dtype(bool): np.dtype("int32"),
+            np.dtype(object): np.dtype(float),
         }
         if datatype:
             dtype_conversions.update(datatype)
diff --git a/cfdm/test/setup_create_field.py b/cfdm/test/setup_create_field.py
index 58eb94330..4367110a2 100644
--- a/cfdm/test/setup_create_field.py
+++ b/cfdm/test/setup_create_field.py
@@ -3,7 +3,7 @@
 import os
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -39,11 +39,11 @@ def setUp(self):
     def test_create_field(self):
         """Test ab initio creation of a first variation of field."""
         # Dimension coordinates
-        dim1 = cfdm.DimensionCoordinate(data=cfdm.Data(numpy.arange(10.0)))
+        dim1 = cfdm.DimensionCoordinate(data=cfdm.Data(np.arange(10.0)))
         dim1.set_property("standard_name", "grid_latitude")
         dim1.set_property("units", "degrees")
 
-        data = numpy.arange(9.0) + 20
+        data = np.arange(9.0) + 20
         data[-1] = 34
         dim0 = cfdm.DimensionCoordinate(data=cfdm.Data(data))
         dim0.set_property("standard_name", "grid_longitude")
@@ -52,7 +52,7 @@ def test_create_field(self):
 
         array = dim0.data.array
 
-        array = numpy.array([array - 0.5, array + 0.5]).transpose((1, 0))
+        array = np.array([array - 0.5, array + 0.5]).transpose((1, 0))
         array[-2, 1] = 30
         array[-1, :] = [30, 36]
         dim0.set_bounds(cfdm.Bounds(data=cfdm.Data(array)))
@@ -77,18 +77,18 @@ def test_create_field(self):
         bk.set_bounds(cfdm.Bounds(data=cfdm.Data([[14, 26.0]])))
 
         aux2 = cfdm.AuxiliaryCoordinate(
-            data=cfdm.Data(numpy.arange(-45, 45, dtype="int32").reshape(10, 9))
+            data=cfdm.Data(np.arange(-45, 45, dtype="int32").reshape(10, 9))
         )
         aux2.set_property("units", "degree_N")
         aux2.set_property("standard_name", "latitude")
 
         aux3 = cfdm.AuxiliaryCoordinate(
-            data=cfdm.Data(numpy.arange(60, 150, dtype="int32").reshape(9, 10))
+            data=cfdm.Data(np.arange(60, 150, dtype="int32").reshape(9, 10))
         )
         aux3.set_property("standard_name", "longitude")
         aux3.set_property("units", "degreeE")
 
-        array = numpy.ma.array(
+        array = np.ma.array(
             [
                 "alpha",
                 "beta",
@@ -102,20 +102,20 @@ def test_create_field(self):
                 "kappa",
             ],
         )
-        array[0] = numpy.ma.masked
+        array[0] = np.ma.masked
         aux4 = cfdm.AuxiliaryCoordinate(data=cfdm.Data(array))
         aux4.set_property("long_name", "greek_letters")
 
         # Cell measures
         msr0 = cfdm.CellMeasure(
-            data=cfdm.Data(1 + numpy.arange(90.0).reshape(9, 10) * 1234)
+            data=cfdm.Data(1 + np.arange(90.0).reshape(9, 10) * 1234)
         )
         msr0.set_measure("area")
         msr0.set_property("units", "km2")
         msr0.nc_set_variable("areacella")
 
         # Data
-        data = cfdm.Data(numpy.arange(90.0).reshape(10, 9))
+        data = cfdm.Data(np.arange(90.0).reshape(10, 9))
 
         properties = {"units": "m s-1"}
 
@@ -201,7 +201,7 @@ def test_create_field(self):
         anc.set_property("foo", "bar")
         f.set_construct(anc, axes=[axisY])
 
-        f.set_property("flag_values", numpy.array([1, 2, 4], "int32"))
+        f.set_property("flag_values", np.array([1, 2, 4], "int32"))
         f.set_property("flag_meanings", "a bb ccc")
         f.set_property("flag_masks", [2, 1, 0])
 
diff --git a/cfdm/test/setup_create_field_2.py b/cfdm/test/setup_create_field_2.py
index b0eb96b27..15eb59a3f 100644
--- a/cfdm/test/setup_create_field_2.py
+++ b/cfdm/test/setup_create_field_2.py
@@ -3,7 +3,7 @@
 import os
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -35,11 +35,11 @@ def setUp(self):
     def test_create_field_2(self):
         """Test ab initio creation of a second variation of field."""
         # Dimension coordinates
-        dim1 = cfdm.DimensionCoordinate(data=cfdm.Data(numpy.arange(10.0)))
+        dim1 = cfdm.DimensionCoordinate(data=cfdm.Data(np.arange(10.0)))
         dim1.set_property("standard_name", "projection_y_coordinate")
         dim1.set_property("units", "m")
 
-        data = numpy.arange(9.0) + 20
+        data = np.arange(9.0) + 20
         data[-1] = 34
         dim0 = cfdm.DimensionCoordinate(data=cfdm.Data(data))
         dim0.set_property("standard_name", "projection_x_coordinate")
@@ -47,7 +47,7 @@ def test_create_field_2(self):
 
         array = dim0.data.array
 
-        array = numpy.array([array - 0.5, array + 0.5]).transpose((1, 0))
+        array = np.array([array - 0.5, array + 0.5]).transpose((1, 0))
         array[-2, 1] = 30
         array[-1, :] = [30, 36]
         dim0.set_bounds(cfdm.Bounds(data=cfdm.Data(array)))
@@ -62,18 +62,18 @@ def test_create_field_2(self):
 
         # Auxiliary coordinates
         aux2 = cfdm.AuxiliaryCoordinate(
-            data=cfdm.Data(numpy.arange(-45, 45, dtype="int32").reshape(10, 9))
+            data=cfdm.Data(np.arange(-45, 45, dtype="int32").reshape(10, 9))
         )
         aux2.set_property("units", "degree_N")
         aux2.set_property("standard_name", "latitude")
 
         aux3 = cfdm.AuxiliaryCoordinate(
-            data=cfdm.Data(numpy.arange(60, 150, dtype="int32").reshape(9, 10))
+            data=cfdm.Data(np.arange(60, 150, dtype="int32").reshape(9, 10))
         )
         aux3.set_property("standard_name", "longitude")
         aux3.set_property("units", "degreeE")
 
-        array = numpy.ma.array(
+        array = np.ma.array(
             [
                 "alpha",
                 "beta",
@@ -87,7 +87,7 @@ def test_create_field_2(self):
                 "kappa",
             ],
         )
-        array[0] = numpy.ma.masked
+        array[0] = np.ma.masked
         aux4 = cfdm.AuxiliaryCoordinate(data=cfdm.Data(array))
         aux4.set_property("standard_name", "greek_letters")
 
@@ -101,13 +101,13 @@ def test_create_field_2(self):
 
         # Cell measures
         msr0 = cfdm.CellMeasure(
-            data=cfdm.Data(1 + numpy.arange(90.0).reshape(9, 10) * 1234)
+            data=cfdm.Data(1 + np.arange(90.0).reshape(9, 10) * 1234)
         )
         msr0.set_measure("area")
         msr0.set_property("units", "km2")
 
         # Data
-        data = cfdm.Data(numpy.arange(90.0).reshape(10, 9))
+        data = cfdm.Data(np.arange(90.0).reshape(10, 9))
 
         properties = {"units": "m s-1"}
 
@@ -219,7 +219,7 @@ def test_create_field_2(self):
         anc.standard_name = "ancillaryC"
         f.set_construct(anc, axes=[axisY])
 
-        f.set_property("flag_values", numpy.array([1, 2, 4], "int32"))
+        f.set_property("flag_values", np.array([1, 2, 4], "int32"))
         f.set_property("flag_meanings", "a bb ccc")
         f.set_property("flag_masks", [2, 1, 0])
 
diff --git a/cfdm/test/setup_create_field_3.py b/cfdm/test/setup_create_field_3.py
index 0e9f58695..5b853913d 100644
--- a/cfdm/test/setup_create_field_3.py
+++ b/cfdm/test/setup_create_field_3.py
@@ -3,7 +3,7 @@
 import os
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -36,7 +36,7 @@ def setUp(self):
     def test_create_field_3(self):
         """Test ab initio creation of a third variation of field."""
         # Dimension coordinates
-        data = numpy.arange(9.0) + 20
+        data = np.arange(9.0) + 20
         data[-1] = 34
         dim0 = cfdm.DimensionCoordinate(data=cfdm.Data(data))
         dim0.set_property("standard_name", "grid_longitude")
@@ -44,12 +44,12 @@ def test_create_field_3(self):
 
         array = dim0.data.array
 
-        array = numpy.array([array - 0.5, array + 0.5]).transpose((1, 0))
+        array = np.array([array - 0.5, array + 0.5]).transpose((1, 0))
         array[-2, 1] = 30
         array[-1, :] = [30, 36]
         dim0.set_bounds(cfdm.Bounds(data=cfdm.Data(array)))
 
-        dim1 = cfdm.DimensionCoordinate(data=cfdm.Data(numpy.arange(10.0)))
+        dim1 = cfdm.DimensionCoordinate(data=cfdm.Data(np.arange(10.0)))
         dim1.set_property("standard_name", "grid_latitude")
         dim1.set_property("units", "degrees")
 
@@ -62,7 +62,7 @@ def test_create_field_3(self):
         )
         dim2.set_property("computed_standard_name", "altitude")
 
-        dim3 = cfdm.DimensionCoordinate(data=cfdm.Data(numpy.array([15.0])))
+        dim3 = cfdm.DimensionCoordinate(data=cfdm.Data(np.array([15.0])))
         dim3.set_property("standard_name", "time")
         dim3.set_property("units", "days since 2004-06-01")
 
@@ -79,18 +79,18 @@ def test_create_field_3(self):
         bk.set_bounds(cfdm.Bounds(data=cfdm.Data([[14, 26.0]])))
 
         aux2 = cfdm.AuxiliaryCoordinate(
-            data=cfdm.Data(numpy.arange(-45, 45, dtype="int32").reshape(10, 9))
+            data=cfdm.Data(np.arange(-45, 45, dtype="int32").reshape(10, 9))
         )
         aux2.set_property("units", "degree_N")
         aux2.set_property("standard_name", "latitude")
 
         aux3 = cfdm.AuxiliaryCoordinate(
-            data=cfdm.Data(numpy.arange(60, 150, dtype="int32").reshape(9, 10))
+            data=cfdm.Data(np.arange(60, 150, dtype="int32").reshape(9, 10))
         )
         aux3.set_property("standard_name", "longitude")
         aux3.set_property("units", "degreeE")
 
-        array = numpy.ma.array(
+        array = np.ma.array(
             [
                 "alpha",
                 "beta",
@@ -104,19 +104,19 @@ def test_create_field_3(self):
                 "kappa",
             ],
         )
-        array[0] = numpy.ma.masked
+        array[0] = np.ma.masked
         aux4 = cfdm.AuxiliaryCoordinate(data=cfdm.Data(array))
         aux4.set_property("standard_name", "greek_letters")
 
         # Cell measures
         msr0 = cfdm.CellMeasure(
-            data=cfdm.Data(1 + numpy.arange(90.0).reshape(9, 10) * 1234)
+            data=cfdm.Data(1 + np.arange(90.0).reshape(9, 10) * 1234)
         )
         msr0.set_measure("area")
         msr0.set_property("units", "km2")
 
         # Data
-        data = cfdm.Data(numpy.arange(90.0).reshape(10, 9))
+        data = cfdm.Data(np.arange(90.0).reshape(10, 9))
 
         properties = {"units": "m s-1"}
 
@@ -219,7 +219,7 @@ def test_create_field_3(self):
         anc.standard_name = "ancillaryC"
         f.set_construct(anc, axes=[axisY])
 
-        f.set_property("flag_values", numpy.array([1, 2, 4], "int32"))
+        f.set_property("flag_values", np.array([1, 2, 4], "int32"))
         f.set_property("flag_meanings", "a bb ccc")
         f.set_property("flag_masks", [2, 1, 0])
 
diff --git a/cfdm/test/test_AuxiliaryCoordinate.py b/cfdm/test/test_AuxiliaryCoordinate.py
index abaefe006..fb868b494 100644
--- a/cfdm/test/test_AuxiliaryCoordinate.py
+++ b/cfdm/test/test_AuxiliaryCoordinate.py
@@ -3,7 +3,7 @@
 import os
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -19,7 +19,7 @@ class AuxiliaryCoordinateTest(unittest.TestCase):
 
     aux1 = cfdm.AuxiliaryCoordinate()
     aux1.standard_name = "latitude"
-    a = numpy.array(
+    a = np.array(
         [
             -30,
             -23.5,
@@ -38,7 +38,7 @@ class AuxiliaryCoordinateTest(unittest.TestCase):
     )
     aux1.set_data(cfdm.Data(a, "degrees_north"))
     bounds = cfdm.Bounds()
-    b = numpy.empty(a.shape + (2,))
+    b = np.empty(a.shape + (2,))
     b[:, 0] = a - 0.1
     b[:, 1] = a + 0.1
     bounds.set_data(cfdm.Data(b))
@@ -115,7 +115,7 @@ def test_AuxiliaryCoordinate_transpose(self):
         x = f.auxiliary_coordinates("longitude").value()
 
         bounds = cfdm.Bounds(
-            data=cfdm.Data(numpy.arange(9 * 10 * 4).reshape(9, 10, 4))
+            data=cfdm.Data(np.arange(9 * 10 * 4).reshape(9, 10, 4))
         )
         x.set_bounds(bounds)
 
@@ -136,7 +136,7 @@ def test_AuxiliaryCoordinate_squeeze(self):
         x = f.auxiliary_coordinates("longitude").value()
 
         bounds = cfdm.Bounds(
-            data=cfdm.Data(numpy.arange(9 * 10 * 4).reshape(9, 10, 4))
+            data=cfdm.Data(np.arange(9 * 10 * 4).reshape(9, 10, 4))
         )
         x.set_bounds(bounds)
         x.insert_dimension(1, inplace=True)
@@ -157,7 +157,7 @@ def test_AuxiliaryCoordinate_interior_ring(self):
         """Test the interior ring access AuxiliaryCoordinate methods."""
         c = cfdm.AuxiliaryCoordinate()
 
-        i = cfdm.InteriorRing(data=cfdm.Data(numpy.arange(10).reshape(5, 2)))
+        i = cfdm.InteriorRing(data=cfdm.Data(np.arange(10).reshape(5, 2)))
 
         c.set_interior_ring(i)
         self.assertTrue(c.has_interior_ring())
diff --git a/cfdm/test/test_Data.py b/cfdm/test/test_Data.py
index 31f93a95f..639446a2c 100644
--- a/cfdm/test/test_Data.py
+++ b/cfdm/test/test_Data.py
@@ -5,7 +5,6 @@
 import os
 import unittest
 
-import numpy
 import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
@@ -44,11 +43,11 @@ def test_Data_any(self):
         """Test the any Data method."""
         d = cfdm.Data([[0, 0, 0]])
         self.assertFalse(d.any())
-        d[0, 0] = numpy.ma.masked
+        d[0, 0] = np.ma.masked
         self.assertFalse(d.any())
         d[0, 1] = 3
         self.assertTrue(d.any())
-        d[...] = numpy.ma.masked
+        d[...] = np.ma.masked
         self.assertFalse(d.any())
 
     def test_Data__repr__str(self):
@@ -71,7 +70,7 @@ def test_Data__repr__str(self):
 
         # Test when the data contains date-times with the first
         # element masked
-        dt = numpy.ma.array([10, 20], mask=[True, False])
+        dt = np.ma.array([10, 20], mask=[True, False])
         d = cfdm.Data(dt, units="days since 2000-01-01")
         self.assertTrue(str(d) == "[--, 2000-01-21 00:00:00]")
         self.assertTrue(repr(d) == "<Data(2): [--, 2000-01-21 00:00:00]>")
@@ -230,8 +229,8 @@ def test_Data__setitem__(self):
 
     def test_Data_apply_masking(self):
         """Test the `apply_masking` Data method."""
-        a = numpy.ma.arange(12).reshape(3, 4)
-        a[1, 1] = numpy.ma.masked
+        a = np.ma.arange(12).reshape(3, 4)
+        a[1, 1] = np.ma.masked
 
         d = cfdm.Data(a, units="m")
 
@@ -243,59 +242,59 @@ def test_Data_apply_masking(self):
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
-        b = numpy.ma.where(a == 0, numpy.ma.masked, a)
+        b = np.ma.where(a == 0, np.ma.masked, a)
         e = d.apply_masking(fill_values=[0])
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
-        b = numpy.ma.where((a == 0) | (a == 11), numpy.ma.masked, a)
+        b = np.ma.where((a == 0) | (a == 11), np.ma.masked, a)
         e = d.apply_masking(fill_values=[0, 11])
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
-        b = numpy.ma.where(a < 3, numpy.ma.masked, a)
+        b = np.ma.where(a < 3, np.ma.masked, a)
         e = d.apply_masking(valid_min=3)
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
-        b = numpy.ma.where(a > 6, numpy.ma.masked, a)
+        b = np.ma.where(a > 6, np.ma.masked, a)
         e = d.apply_masking(valid_max=6)
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
-        b = numpy.ma.where((a < 2) | (a > 8), numpy.ma.masked, a)
+        b = np.ma.where((a < 2) | (a > 8), np.ma.masked, a)
         e = d.apply_masking(valid_range=[2, 8])
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
         d.set_fill_value(7)
 
-        b = numpy.ma.where(a == 7, numpy.ma.masked, a)
+        b = np.ma.where(a == 7, np.ma.masked, a)
         e = d.apply_masking(fill_values=True)
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
-        b = numpy.ma.where((a == 7) | (a < 2) | (a > 8), numpy.ma.masked, a)
+        b = np.ma.where((a == 7) | (a < 2) | (a > 8), np.ma.masked, a)
         e = d.apply_masking(fill_values=True, valid_range=[2, 8])
         self.assertTrue((b == e.array).all())
         self.assertTrue((b.mask == e.mask.array).all())
 
     #    def test_Data_astype(self):
-    #        a = numpy.array([1.5, 2, 2.5], dtype=float)
+    #        a = np.array([1.5, 2, 2.5], dtype=float)
     #        d = cfdm.Data(a)
     #
-    #        self.assertTrue(d.dtype == numpy.dtype(float))
-    #        self.assertTrue(d.array.dtype == numpy.dtype(float))
+    #        self.assertTrue(d.dtype == np.dtype(float))
+    #        self.assertTrue(d.array.dtype == np.dtype(float))
     #        self.assertTrue((d.array == a).all())
     #
     #        d.astype('int32')
-    #        self.assertTrue(d.dtype == numpy.dtype('int32'))
-    #        self.assertTrue(d.array.dtype == numpy.dtype('int32'))
+    #        self.assertTrue(d.dtype == np.dtype('int32'))
+    #        self.assertTrue(d.array.dtype == np.dtype('int32'))
     #        self.assertTrue((d.array == [1, 2, 2]).all())
     #
     #        d = cfdm.Data(a)
     #        try:
-    #            d.astype(numpy.dtype(int, casting='safe'))
+    #            d.astype(np.dtype(int, casting='safe'))
     #            self.assertTrue(False)
     #        except TypeError:
     #            pass
@@ -305,19 +304,19 @@ def test_Data_array(self):
         # ------------------------------------------------------------
         # Numpy array interface (__array__)
         # ------------------------------------------------------------
-        a = numpy.arange(12, dtype="int32").reshape(3, 4)
+        a = np.arange(12, dtype="int32").reshape(3, 4)
 
         d = cfdm.Data(a, units="km")
 
-        b = numpy.array(d)
+        b = np.array(d)
 
-        self.assertEqual(b.dtype, numpy.dtype("int32"))
+        self.assertEqual(b.dtype, np.dtype("int32"))
         self.assertEqual(a.shape, b.shape)
         self.assertTrue((a == b).all())
 
-        b = numpy.array(d, dtype="float32")
+        b = np.array(d, dtype="float32")
 
-        self.assertEqual(b.dtype, numpy.dtype("float32"))
+        self.assertEqual(b.dtype, np.dtype("float32"))
         self.assertEqual(a.shape, b.shape)
         self.assertTrue((a == b).all())
 
@@ -325,14 +324,14 @@ def test_Data_array(self):
         d = cfdm.Data(9, units="km")
         a = d.array
         self.assertEqual(a.shape, ())
-        self.assertEqual(a, numpy.array(9))
+        self.assertEqual(a, np.array(9))
         d[...] = cfdm.masked
         a = d.array
         self.assertEqual(a.shape, ())
-        self.assertIs(a[()], numpy.ma.masked)
+        self.assertIs(a[()], np.ma.masked)
 
         # Non-scalar numeric array
-        b = numpy.arange(10 * 15 * 19).reshape(10, 1, 15, 19)
+        b = np.arange(10 * 15 * 19).reshape(10, 1, 15, 19)
         d = cfdm.Data(b, "km")
         a = d.array
         a[0, 0, 0, 0] = -999
@@ -351,7 +350,7 @@ def test_Data_datetime_array(self):
 
         d[0] = cfdm.masked
         dt = d.datetime_array
-        self.assertIs(dt[0], numpy.ma.masked)
+        self.assertIs(dt[0], np.ma.masked)
         self.assertEqual(dt[1], datetime.datetime(2000, 12, 2, 0, 0))
 
         d = cfdm.Data(11292.5, units="days since 1970-1-1")
@@ -360,11 +359,11 @@ def test_Data_datetime_array(self):
 
         d[()] = cfdm.masked
         dt = d.datetime_array
-        self.assertIs(dt[()], numpy.ma.masked)
+        self.assertIs(dt[()], np.ma.masked)
 
     def test_Data_flatten(self):
         """Test the flatten Data method."""
-        ma = numpy.ma.arange(24).reshape(1, 2, 3, 4)
+        ma = np.ma.arange(24).reshape(1, 2, 3, 4)
         ma[0, 1, 1, 2] = cfdm.masked
         ma[0, 0, 2, 1] = cfdm.masked
 
@@ -390,9 +389,7 @@ def test_Data_flatten(self):
                 shape = [n for i, n in enumerate(d.shape) if i not in axes]
                 shape.insert(
                     sorted(axes)[0],
-                    numpy.prod(
-                        [n for i, n in enumerate(d.shape) if i in axes]
-                    ),
+                    np.prod([n for i, n in enumerate(d.shape) if i in axes]),
                 )
 
             self.assertEqual(e.shape, tuple(shape), axes)
@@ -416,12 +413,12 @@ def test_Data_flatten(self):
 
     def test_Data_transpose(self):
         """Test the transpose Data method."""
-        a = numpy.arange(2 * 3 * 5).reshape(2, 1, 3, 5)
+        a = np.arange(2 * 3 * 5).reshape(2, 1, 3, 5)
         d = cfdm.Data(a.copy())
 
         for indices in (list(range(a.ndim)), list(range(-a.ndim, 0))):
             for axes in itertools.permutations(indices):
-                a = numpy.transpose(a, axes)
+                a = np.transpose(a, axes)
                 d = d.transpose(axes)
                 message = (
                     f"cfdm.Data.transpose({axes}) failed: d.shape={d.shape}, "
@@ -455,8 +452,8 @@ def test_Data_unique(self):
 
     def test_Data_equals(self):
         """Test the equality-testing Data method."""
-        a = numpy.ma.arange(10 * 15 * 19).reshape(10, 1, 15, 19)
-        a[0, 0, 2, 3] = numpy.ma.masked
+        a = np.ma.arange(10 * 15 * 19).reshape(10, 1, 15, 19)
+        a[0, 0, 2, 3] = np.ma.masked
 
         d = cfdm.Data(a, units="days since 2000-2-2", calendar="noleap")
         e = copy.deepcopy(d)
@@ -467,10 +464,10 @@ def test_Data_equals(self):
 
     def test_Data_maximum_minimum_sum_squeeze(self):
         """Test the maximum, minimum, sum and squeeze Data methods."""
-        a = numpy.ma.arange(2 * 3 * 5).reshape(2, 1, 3, 5)
-        a[0, 0, 0, 0] = numpy.ma.masked
-        a[-1, -1, -1, -1] = numpy.ma.masked
-        a[1, -1, 1, 3] = numpy.ma.masked
+        a = np.ma.arange(2 * 3 * 5).reshape(2, 1, 3, 5)
+        a[0, 0, 0, 0] = np.ma.masked
+        a[-1, -1, -1, -1] = np.ma.masked
+        a[1, -1, 1, 3] = np.ma.masked
         d = cfdm.Data(a)
 
         b = a.max()
@@ -542,14 +539,14 @@ def test_Data_maximum_minimum_sum_squeeze(self):
 
     def test_Data_dtype_mask(self):
         """Test the dtype and mask Data methods."""
-        a = numpy.ma.array(
+        a = np.ma.array(
             [[280.0, -99, -99, -99], [281.0, 279.0, 278.0, 279.0]],
             dtype=float,
             mask=[[0, 1, 1, 1], [0, 0, 0, 0]],
         )
 
         d = cfdm.Data([[280, -99, -99, -99], [281, 279, 278, 279]])
-        self.assertEqual(d.dtype, numpy.dtype(int))
+        self.assertEqual(d.dtype, np.dtype(int))
 
         d = cfdm.Data(
             [[280, -99, -99, -99], [281, 279, 278, 279]],
@@ -560,12 +557,12 @@ def test_Data_dtype_mask(self):
         self.assertEqual(d.dtype, a.dtype)
         self.assertTrue((d.array == a).all())
         self.assertEqual(d.mask.shape, a.mask.shape)
-        self.assertTrue((d.mask.array == numpy.ma.getmaskarray(a)).all())
+        self.assertTrue((d.mask.array == np.ma.getmaskarray(a)).all())
 
-        a = numpy.array(
+        a = np.array(
             [[280.0, -99, -99, -99], [281.0, 279.0, 278.0, 279.0]], dtype=float
         )
-        mask = numpy.ma.masked_all(a.shape).mask
+        mask = np.ma.masked_all(a.shape).mask
 
         d = cfdm.Data(
             [[280, -99, -99, -99], [281, 279, 278, 279]], dtype=float
@@ -574,7 +571,7 @@ def test_Data_dtype_mask(self):
         self.assertEqual(d.dtype, a.dtype)
         self.assertTrue((d.array == a).all())
         self.assertEqual(d.mask.shape, mask.shape)
-        self.assertTrue((d.mask.array == numpy.ma.getmaskarray(a)).all())
+        self.assertTrue((d.mask.array == np.ma.getmaskarray(a)).all())
 
     def test_Data_get_index(self):
         """Test the `get_index` Data method."""
@@ -644,12 +641,12 @@ def test_Data_insert_dimension(self):
         e = d.squeeze([0, 2])
         self.assertEqual(e.shape, (12, 1))
 
-        array = numpy.arange(12).reshape(1, 4, 3)
+        array = np.arange(12).reshape(1, 4, 3)
         d = cfdm.Data(array)
         e = d.squeeze()
         f = e.insert_dimension(0)
         a = f.array
-        self.assertTrue(numpy.allclose(a, array))
+        self.assertTrue(np.allclose(a, array))
 
         with self.assertRaises(ValueError):
             d.insert_dimension(1000)
diff --git a/cfdm/test/test_DomainAncillary.py b/cfdm/test/test_DomainAncillary.py
index b11869256..981eb98cc 100644
--- a/cfdm/test/test_DomainAncillary.py
+++ b/cfdm/test/test_DomainAncillary.py
@@ -3,7 +3,7 @@
 import os
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -88,7 +88,7 @@ def test_DomainAncillary_transpose(self):
         f = cfdm.read(self.filename)[0]
         a = f.auxiliary_coordinates("longitude").value()
         bounds = cfdm.Bounds(
-            data=cfdm.Data(numpy.arange(9 * 10 * 4).reshape(9, 10, 4))
+            data=cfdm.Data(np.arange(9 * 10 * 4).reshape(9, 10, 4))
         )
         a.set_bounds(bounds)
         x = cfdm.DomainAncillary(source=a)
@@ -109,7 +109,7 @@ def test_DomainAncillary_squeeze(self):
         f = cfdm.read(self.filename)[0]
         a = f.auxiliary_coordinates("longitude").value()
         bounds = cfdm.Bounds(
-            data=cfdm.Data(numpy.arange(9 * 10 * 4).reshape(9, 10, 4))
+            data=cfdm.Data(np.arange(9 * 10 * 4).reshape(9, 10, 4))
         )
         a.set_bounds(bounds)
         x = cfdm.DomainAncillary(source=a)
diff --git a/cfdm/test/test_NumpyArray.py b/cfdm/test/test_NumpyArray.py
index cbab37bdb..3ae3c11ab 100644
--- a/cfdm/test/test_NumpyArray.py
+++ b/cfdm/test/test_NumpyArray.py
@@ -3,7 +3,7 @@
 import faulthandler
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -27,7 +27,7 @@ def setUp(self):
 
     def test_NumpyArray_copy(self):
         """Test the copy module copying behaviour of NumpyArray."""
-        a = numpy.array([1, 2, 3, 4])
+        a = np.array([1, 2, 3, 4])
 
         x = cfdm.NumpyArray(a)
         y = copy.deepcopy(x)
@@ -36,11 +36,11 @@ def test_NumpyArray_copy(self):
 
     def test_NumpyArray__array__(self):
         """Test the NumPy array conversion of NumpyArray."""
-        a = numpy.array([1, 2, 3, 4])
+        a = np.array([1, 2, 3, 4])
 
         x = cfdm.NumpyArray(a)
 
-        b = numpy.array(x)
+        b = np.array(x)
         self.assertTrue((b == a).all())
 
     def test_NumpyArray_get_filename(self):
diff --git a/cfdm/test/test_core_1.py b/cfdm/test/test_core_1.py
index f30fb4fe6..340ae2974 100644
--- a/cfdm/test/test_core_1.py
+++ b/cfdm/test/test_core_1.py
@@ -2,7 +2,7 @@
 import faulthandler
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -31,12 +31,12 @@ def test_core_create_field(self):
         """Test ab initio creation of a `cfdm.core` Field."""
         # Dimension coordinates
         dim1 = cfdm.core.DimensionCoordinate(
-            data=cfdm.core.Data(cfdm.core.NumpyArray(numpy.arange(10.0)))
+            data=cfdm.core.Data(cfdm.core.NumpyArray(np.arange(10.0)))
         )
         dim1.set_property("standard_name", "grid_latitude")
         dim1.set_property("units", "degrees")
 
-        data = numpy.arange(9.0) + 20
+        data = np.arange(9.0) + 20
         data[-1] = 34
         dim0 = cfdm.core.DimensionCoordinate(
             data=cfdm.core.Data(cfdm.core.NumpyArray(data))
@@ -46,18 +46,16 @@ def test_core_create_field(self):
 
         array = dim0.data.array
 
-        array = numpy.array([array - 0.5, array + 0.5]).transpose((1, 0))
+        array = np.array([array - 0.5, array + 0.5]).transpose((1, 0))
         array[-2, 1] = 30
         array[-1, :] = [30, 36]
         data = cfdm.core.Data(cfdm.core.NumpyArray(array))
         dim0.set_bounds(cfdm.core.Bounds(data=data))
 
         dim2 = cfdm.core.DimensionCoordinate(
-            data=cfdm.core.Data(cfdm.core.NumpyArray(numpy.array([1.5]))),
+            data=cfdm.core.Data(cfdm.core.NumpyArray(np.array([1.5]))),
             bounds=cfdm.core.Bounds(
-                data=cfdm.core.Data(
-                    cfdm.core.NumpyArray(numpy.array([[1, 2.0]]))
-                )
+                data=cfdm.core.Data(cfdm.core.NumpyArray(np.array([[1, 2.0]])))
             ),
         )
         dim2.set_property(
@@ -67,24 +65,24 @@ def test_core_create_field(self):
 
         # Auxiliary coordinates
         ak = cfdm.core.DomainAncillary(
-            data=cfdm.core.Data(cfdm.core.NumpyArray(numpy.array([10.0])))
+            data=cfdm.core.Data(cfdm.core.NumpyArray(np.array([10.0])))
         )
         ak.set_property("units", "m")
         ak.set_bounds(
             cfdm.core.Bounds(
                 data=cfdm.core.Data(
-                    cfdm.core.NumpyArray(numpy.array([[5, 15.0]]))
+                    cfdm.core.NumpyArray(np.array([[5, 15.0]]))
                 )
             )
         )
 
         bk = cfdm.core.DomainAncillary(
-            data=cfdm.core.Data(cfdm.core.NumpyArray(numpy.array([20.0])))
+            data=cfdm.core.Data(cfdm.core.NumpyArray(np.array([20.0])))
         )
         bk.set_bounds(
             cfdm.core.Bounds(
                 data=cfdm.core.Data(
-                    cfdm.core.NumpyArray(numpy.array([[14, 26.0]]))
+                    cfdm.core.NumpyArray(np.array([[14, 26.0]]))
                 )
             )
         )
@@ -92,7 +90,7 @@ def test_core_create_field(self):
         aux2 = cfdm.core.AuxiliaryCoordinate(
             data=cfdm.core.Data(
                 cfdm.core.NumpyArray(
-                    numpy.arange(-45, 45, dtype="int32").reshape(10, 9)
+                    np.arange(-45, 45, dtype="int32").reshape(10, 9)
                 )
             )
         )
@@ -102,14 +100,14 @@ def test_core_create_field(self):
         aux3 = cfdm.core.AuxiliaryCoordinate(
             data=cfdm.core.Data(
                 cfdm.core.NumpyArray(
-                    numpy.arange(60, 150, dtype="int32").reshape(9, 10)
+                    np.arange(60, 150, dtype="int32").reshape(9, 10)
                 )
             )
         )
         aux3.set_property("standard_name", "longitude")
         aux3.set_property("units", "degreeE")
 
-        array = numpy.ma.array(
+        array = np.ma.array(
             [
                 "alpha",
                 "beta",
@@ -124,7 +122,7 @@ def test_core_create_field(self):
             ],
             dtype="S",
         )
-        array[0] = numpy.ma.masked
+        array[0] = np.ma.masked
         aux4 = cfdm.core.AuxiliaryCoordinate(
             data=cfdm.core.Data(cfdm.core.NumpyArray(array))
         )
@@ -133,9 +131,7 @@ def test_core_create_field(self):
         # Cell measures
         msr0 = cfdm.core.CellMeasure(
             data=cfdm.core.Data(
-                cfdm.core.NumpyArray(
-                    1 + numpy.arange(90.0).reshape(9, 10) * 1234
-                )
+                cfdm.core.NumpyArray(1 + np.arange(90.0).reshape(9, 10) * 1234)
             )
         )
         msr0.set_measure("area")
@@ -143,7 +139,7 @@ def test_core_create_field(self):
 
         # Data
         data = cfdm.core.Data(
-            cfdm.core.NumpyArray(numpy.arange(90.0).reshape(10, 9))
+            cfdm.core.NumpyArray(np.arange(90.0).reshape(10, 9))
         )
 
         properties = {"units": "m s-1"}
@@ -229,7 +225,7 @@ def test_core_create_field(self):
         )
         f.set_construct(anc, axes=[axisY])
 
-        f.set_property("flag_values", numpy.array([1, 2, 4], "int32"))
+        f.set_property("flag_values", np.array([1, 2, 4], "int32"))
         f.set_property("flag_meanings", "a bb ccc")
         f.set_property("flag_masks", [2, 1, 0])
 
@@ -238,7 +234,7 @@ def test_core_create_field(self):
             method="mean",
             qualifiers={
                 "interval": [
-                    cfdm.core.Data(cfdm.core.NumpyArray(numpy.array(1)), "day")
+                    cfdm.core.Data(cfdm.core.NumpyArray(np.array(1)), "day")
                 ],
                 "comment": "ok",
             },
diff --git a/cfdm/test/test_core_2.py b/cfdm/test/test_core_2.py
index e8cf4886f..4004c14c1 100644
--- a/cfdm/test/test_core_2.py
+++ b/cfdm/test/test_core_2.py
@@ -2,7 +2,7 @@
 import datetime
 import unittest
 
-import numpy
+import numpy as np
 
 import cfdm
 
@@ -25,14 +25,14 @@ def setUp(self):
 
     def test_core_NumpyArray(self):
         """Test cfdm.core.NumpyArray class."""
-        a = cfdm.core.NumpyArray(numpy.array([1, 2, 3]))
+        a = cfdm.core.NumpyArray(np.array([1, 2, 3]))
 
         # __deepcopy__
         copy.deepcopy(a)
 
         # __array__
-        numpy.array(a)
-        numpy.array(a, dtype="float")
+        np.array(a)
+        np.array(a, dtype="float")
 
     def test_core_Container(self):
         """Test cfdm.core.Container class."""
diff --git a/cfdm/test/test_dsg.py b/cfdm/test/test_dsg.py
index 5cfc67f5f..8be4980bb 100644
--- a/cfdm/test/test_dsg.py
+++ b/cfdm/test/test_dsg.py
@@ -5,7 +5,7 @@
 import tempfile
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -50,7 +50,7 @@ class DSGTest(unittest.TestCase):
     i = cfdm.read(indexed)
     ic = cfdm.read(indexed_contiguous)
 
-    a = numpy.ma.masked_all((4, 9), dtype=float)
+    a = np.ma.masked_all((4, 9), dtype=float)
     a[0, 0:3] = [0.0, 1.0, 2.0]
     a[1, 0:7] = [1.0, 11.0, 21.0, 31.0, 41.0, 51.0, 61.0]
     a[2, 0:5] = [2.0, 102.0, 202.0, 302.0, 402.0]
@@ -66,7 +66,7 @@ class DSGTest(unittest.TestCase):
         8003.0,
     ]
 
-    b = numpy.array(
+    b = np.array(
         [
             [
                 [20.7, -99, -99, -99],
@@ -155,7 +155,7 @@ class DSGTest(unittest.TestCase):
         ]
     )
 
-    b = numpy.ma.where(b == -99, numpy.ma.masked, b)
+    b = np.ma.where(b == -99, np.ma.masked, b)
 
     def setUp(self):
         """Preparations called immediately before each test method."""
@@ -197,7 +197,7 @@ def test_DSG_contiguous(self):
         # Test creation
         # ------------------------------------------------------------
         # Define the ragged array values
-        ragged_array = numpy.array(
+        ragged_array = np.array(
             [280, 282.5, 281, 279, 278, 279.5], dtype="float32"
         )
 
@@ -295,7 +295,7 @@ def test_DSG_indexed_contiguous(self):
     def test_DSG_create_contiguous(self):
         """Test the creation of a contiguous ragged array."""
         # Define the ragged array values
-        ragged_array = numpy.array([1, 3, 4, 3, 6], dtype="float32")
+        ragged_array = np.array([1, 3, 4, 3, 6], dtype="float32")
         # Define the count array values
         count_array = [2, 3]
 
@@ -328,7 +328,7 @@ def test_DSG_create_contiguous(self):
         self.assertTrue(
             (
                 z.data.array
-                == numpy.ma.masked_array(
+                == np.ma.masked_array(
                     data=[[1.0, 3.0, 99], [4.0, 3.0, 6.0]],
                     mask=[[False, False, True], [False, False, False]],
                     fill_value=1e20,
@@ -342,12 +342,12 @@ def test_DSG_create_contiguous(self):
         self.assertTrue(
             (
                 z.data.compressed_array
-                == numpy.array([1.0, 3.0, 4.0, 3.0, 6.0], dtype="float32")
+                == np.array([1.0, 3.0, 4.0, 3.0, 6.0], dtype="float32")
             ).all()
         )
 
         self.assertTrue(
-            (z.data.get_count().data.array == numpy.array([2, 3])).all()
+            (z.data.get_count().data.array == np.array([2, 3])).all()
         )
 
 
diff --git a/cfdm/test/test_functions.py b/cfdm/test/test_functions.py
index 8a5552615..c0de26701 100644
--- a/cfdm/test/test_functions.py
+++ b/cfdm/test/test_functions.py
@@ -9,7 +9,7 @@
 import tempfile
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -506,13 +506,13 @@ def test_Constant(self):
         self.assertEqual(20, c)
         self.assertEqual(c, c)
         self.assertEqual(c, copy.deepcopy(c))
-        self.assertEqual(c, numpy.array(20))
+        self.assertEqual(c, np.array(20))
 
         self.assertNotEqual(c, 999)
         self.assertNotEqual(999, c)
         self.assertNotEqual(c, d)
-        self.assertNotEqual(c, numpy.array(999))
-        self.assertNotEqual(numpy.array(999), c)
+        self.assertNotEqual(c, np.array(999))
+        self.assertNotEqual(np.array(999), c)
 
         self.assertLess(c, 999)
         self.assertLessEqual(c, 999)
diff --git a/cfdm/test/test_gathering.py b/cfdm/test/test_gathering.py
index d200c00a4..1492900ca 100644
--- a/cfdm/test/test_gathering.py
+++ b/cfdm/test/test_gathering.py
@@ -5,7 +5,7 @@
 import tempfile
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -49,7 +49,7 @@ def setUp(self):
 
         self.gathered = "gathered.nc"
 
-        a = numpy.ma.masked_all((4, 9), dtype=float)
+        a = np.ma.masked_all((4, 9), dtype=float)
         a[0, 0:3] = [0.0, 1.0, 2.0]
         a[1, 0:7] = [1.0, 11.0, 21.0, 31.0, 41.0, 51.0, 61.0]
         a[2, 0:5] = [2.0, 102.0, 202.0, 302.0, 402.0]
@@ -66,7 +66,7 @@ def setUp(self):
         ]
         self.a = a
 
-        b = numpy.ma.array(
+        b = np.ma.array(
             [
                 [
                     [207.12345561172262, -99, -99, -99],
@@ -250,7 +250,7 @@ def setUp(self):
             ]
         )
 
-        b = numpy.ma.where(b == -99, numpy.ma.masked, b)
+        b = np.ma.where(b == -99, np.ma.masked, b)
         self.b = b
 
     def test_GATHERING(self):
@@ -278,7 +278,7 @@ def test_GATHERING(self):
     def test_GATHERING_create(self):
         """Test the creation of a construct with a gathered array."""
         # Define the gathered values
-        gathered_array = numpy.array(
+        gathered_array = np.array(
             [[280, 282.5, 281], [279, 278, 277.5]], dtype="float32"
         )
         # Define the list array values
@@ -314,7 +314,7 @@ def test_GATHERING_create(self):
         self.assertTrue(
             (
                 tas.data.array
-                == numpy.ma.masked_array(
+                == np.ma.masked_array(
                     data=[
                         [[1, 280.0], [1, 1], [282.5, 281.0]],
                         [[1, 279.0], [1, 1], [278.0, 277.5]],
@@ -334,7 +334,7 @@ def test_GATHERING_create(self):
         self.assertTrue(
             (
                 tas.data.compressed_array
-                == numpy.array(
+                == np.array(
                     [[280.0, 282.5, 281.0], [279.0, 278.0, 277.5]],
                     dtype="float32",
                 )
@@ -342,7 +342,7 @@ def test_GATHERING_create(self):
         )
 
         self.assertTrue(
-            (tas.data.get_list().data.array == numpy.array([1, 4, 5])).all()
+            (tas.data.get_list().data.array == np.array([1, 4, 5])).all()
         )
 
 
diff --git a/cfdm/test/test_string.py b/cfdm/test/test_string.py
index 2681673f8..ff54b4d43 100644
--- a/cfdm/test/test_string.py
+++ b/cfdm/test/test_string.py
@@ -5,7 +5,7 @@
 import tempfile
 import unittest
 
-import numpy
+import numpy as np
 
 faulthandler.enable()  # to debug seg faults and timeouts
 
@@ -50,8 +50,8 @@ def setUp(self):
     def test_STRING(self):
         """Test constructs with underlying string type arrays."""
         for array in (
-            numpy.ma.array(list("abcdefghij"), dtype="S"),
-            numpy.ma.array(
+            np.ma.array(list("abcdefghij"), dtype="S"),
+            np.ma.array(
                 ["a", "b1", "c12", "d123", "e1234", "f", "g", "h", "i", "j"],
                 dtype="S",
             ),
@@ -73,7 +73,7 @@ def test_STRING(self):
 
             # Set the field data
             tas.set_data(
-                cfdm.Data(numpy.arange(90.0).reshape(10, 9)),
+                cfdm.Data(np.arange(90.0).reshape(10, 9)),
                 axes=[axis_Y, axis_X],
             )
 
@@ -83,9 +83,9 @@ def test_STRING(self):
                     "standard_name": "grid_latitude",
                     "units": "degrees",
                 },
-                data=cfdm.Data(numpy.arange(10.0)),
+                data=cfdm.Data(np.arange(10.0)),
                 bounds=cfdm.Bounds(
-                    data=cfdm.Data(numpy.arange(20).reshape(10, 2))
+                    data=cfdm.Data(np.arange(20).reshape(10, 2))
                 ),
             )
 
@@ -94,9 +94,9 @@ def test_STRING(self):
                     "standard_name": "grid_longitude",
                     "units": "degrees",
                 },
-                data=cfdm.Data(numpy.arange(9.0)),
+                data=cfdm.Data(np.arange(9.0)),
                 bounds=cfdm.Bounds(
-                    data=cfdm.Data(numpy.arange(18).reshape(9, 2))
+                    data=cfdm.Data(np.arange(18).reshape(9, 2))
                 ),
             )
 
@@ -104,7 +104,7 @@ def test_STRING(self):
             tas.set_construct(dimension_coordinate_X, axes=[axis_X])
 
             # Create and set the auxiliary coordinates
-            array[0] = numpy.ma.masked
+            array[0] = np.ma.masked
             aux0 = cfdm.AuxiliaryCoordinate(
                 properties={"long_name": "Grid latitude name"},
                 data=cfdm.Data(array),