diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 6c8862394e6cc4b415c5ed8e2208bcab94674c19..f7124b6641a6c8b6613f3f6a10528f667e478f87 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -17,6 +17,16 @@ Added
 * New :func:`.affine.rescale` function, for adjusting a scaling matrix.
 
 
+Changed
+^^^^^^^
+
+
+* The :class:`.ImageWrapper` now maintains its own image data cache, rather
+  than depending on ``nibabel``.
+* Internal changes to avoid using the deprecated
+  ``nibabel.dataobj_images.DataobjImage.get_data`` method.
+
+
 Fixed
 ^^^^^
 
diff --git a/fsl/data/dicom.py b/fsl/data/dicom.py
index d38d2081d744a27c9edd71902fc3f1909dafa6d9..51c750385796b0967a25491e6e28385420a82218 100644
--- a/fsl/data/dicom.py
+++ b/fsl/data/dicom.py
@@ -36,6 +36,7 @@ import               shlex
 import               logging
 import               binascii
 
+import numpy      as np
 import nibabel    as nib
 
 import fsl.utils.tempdir as tempdir
@@ -279,9 +280,11 @@ def loadSeries(series):
 
         # copy images so nibabel no longer
         # refs to the files (as they will
-        # be deleted), and use get_data()
-        # to force-load the image data.
-        images = [nib.Nifti1Image(i.get_data(), None, i.header)
+        # be deleted), and force-load the
+        # the image data into memory (to
+        # avoid any disk accesses due to
+        # e.g. memmap)
+        images = [nib.Nifti1Image(np.asanyarray(i.dataobj), None, i.header)
                   for i in images]
 
         return [DicomImage(i, series, dcmdir, name=desc) for i in images]
diff --git a/fsl/data/image.py b/fsl/data/image.py
index 39f1f7facdd789e197340177f8e73dea77507c31..46fb47f813b3dbca90888a45a3a01245f635a350 100644
--- a/fsl/data/image.py
+++ b/fsl/data/image.py
@@ -1204,6 +1204,9 @@ class Image(Nifti):
     @property
     def nibImage(self):
         """Returns a reference to the ``nibabel`` NIFTI image instance.
+        Note that if the image data has been modified through this ``Image``,
+        it will be out of sync with what is returned by the ``nibabel`` object,
+        until a call to :meth:`save` is made.
         """
         return self.__nibImage
 
@@ -1387,6 +1390,21 @@ class Image(Nifti):
         os.close(tmphd)
 
         try:
+            # First of all, the nibabel object won't know
+            # about any image data modifications, so if
+            # any have occurred, we need to create a new
+            # nibabel image using the data managed by the
+            # imagewrapper, and the old header.
+            #
+            # Assuming here that analyze/nifti1/nifti2
+            # nibabel classes have an __init__ which
+            # expects (data, affine, header)
+            if not self.saveState:
+                self.__nibImage = type(self.__nibImage)(self[:],
+                                                        None,
+                                                        self.header)
+                self.header     = self.__nibImage.header
+
             nib.save(self.__nibImage, tmpfname)
 
             # nibabel should close any old
@@ -1395,6 +1413,8 @@ class Image(Nifti):
             self.__nibImage = None
             self.header     = None
 
+            # Copy to final destination,
+            # and reload from there
             imcp.imcp(tmpfname, filename, overwrite=True)
 
             self.__nibImage = nib.load(filename)
diff --git a/fsl/data/imagewrapper.py b/fsl/data/imagewrapper.py
index 8afcccacf7ad3e9b4e6ce035aabffa6393d4cf70..b7d5b3577c903a9bf434d14cbc8647a243c4bfd2 100644
--- a/fsl/data/imagewrapper.py
+++ b/fsl/data/imagewrapper.py
@@ -76,8 +76,7 @@ class ImageWrapper(notifier.Notifier):
      - The image data is not modified (via :meth:`__setitem__`.
 
     If any of these conditions do not hold, the image data will be loaded into
-    memory and accessed directly, via the ``nibabel.Nifti1Image.get_data``
-    method.
+    memory and accessed directly.
 
 
     *Image dimensionality*
@@ -218,7 +217,11 @@ class ImageWrapper(notifier.Notifier):
 
         self.reset(dataRange)
 
-        if loadData:
+        # We keep an internal ref to
+        # the data numpy array if/when
+        # it is loaded in memory
+        self.__data = None
+        if loadData or image.in_memory:
             self.loadData()
 
         if threaded:
@@ -232,6 +235,7 @@ class ImageWrapper(notifier.Notifier):
         the :class:`.TaskThread` is stopped.
         """
         self.__image = None
+        self.__data  = None
         if self.__taskThread is not None:
             self.__taskThread.stop()
             self.__taskThread = None
@@ -377,16 +381,11 @@ class ImageWrapper(notifier.Notifier):
         """Forces all of the image data to be loaded into memory.
 
         .. note:: This method will be called by :meth:`__init__` if its
-                  ``loadData`` parameter is ``True``.
+                  ``loadData`` parameter is ``True``. It will also be called
+                  on all write operations (see :meth:`__setitem__`).
         """
-
-        # If the data is not already loaded, this will
-        # cause nibabel to load it. By default, nibabel
-        # will cache the numpy array that contains the
-        # image data, so subsequent calls to this
-        # method will not overwrite any changes that
-        # have been made to the data array.
-        self.__image.get_data()
+        if self.__data is None:
+            self.__data = np.asanyarray(self.__image.dataobj)
 
 
     def __getData(self, sliceobj, isTuple=False):
@@ -407,14 +406,7 @@ class ImageWrapper(notifier.Notifier):
         # ArrayProxy. Otheriwse if it is in
         # memory, we can access it directly.
         #
-        # Furthermore, if it is in memory and
-        # has been modified, the ArrayProxy
-        # will give us out-of-date values (as
-        # the ArrayProxy reads from disk). So
-        # we have to read from the in-memory
-        # array to get changed values.
-        #
-        # Finally, note that if the caller has
+        # Note also that if the caller has
         # given us a 'fancy' slice object (a
         # boolean numpy array), but the image
         # data is not in memory, we can't access
@@ -422,8 +414,8 @@ class ImageWrapper(notifier.Notifier):
         # (the dataobj attribute) cannot handle
         # fancy indexing. In this case an error
         # will be raised.
-        if self.__image.in_memory: return self.__image.get_data()[sliceobj]
-        else:                      return self.__image.dataobj[   sliceobj]
+        if self.__data is not None: return self.__data[         sliceobj]
+        else:                       return self.__image.dataobj[sliceobj]
 
 
     def __imageIsCovered(self):
@@ -444,18 +436,18 @@ class ImageWrapper(notifier.Notifier):
         _, expansions = calcExpansion(slices, self.__coverage)
         expansions    = collapseExpansions(expansions, self.__numRealDims - 1)
 
-        log.debug('Updating image {} data range [slice: {}] '
-                  '(current range: [{}, {}]; '
-                  'number of expansions: {}; '
-                  'current coverage: {}; '
-                  'volume ranges: {})'.format(
-                      self.__name,
-                      slices,
-                      self.__range[0],
-                      self.__range[1],
-                      len(expansions),
-                      self.__coverage,
-                      self.__volRanges))
+        log.debug('Updating image %s data range [slice: %s] '
+                  '(current range: [%s, %s]; '
+                  'number of expansions: %s; '
+                  'current coverage: %s; '
+                  'volume ranges: %s)',
+                  self.__name,
+                  slices,
+                  self.__range[0],
+                  self.__range[1],
+                  len(expansions),
+                  self.__coverage,
+                  self.__volRanges)
 
         # As we access the data for each expansions,
         # we want it to have the same dimensionality
@@ -507,12 +499,12 @@ class ImageWrapper(notifier.Notifier):
 
         if any((oldmin is None, oldmax is None)) or \
            not np.all(np.isclose([oldmin, oldmax], [newmin, newmax])):
-            log.debug('Image {} range changed: [{}, {}] -> [{}, {}]'.format(
-                self.__name,
-                oldmin,
-                oldmax,
-                newmin,
-                newmax))
+            log.debug('Image %s range changed: [%s, %s] -> [%s, %s]',
+                      self.__name,
+                      oldmin,
+                      oldmax,
+                      newmin,
+                      newmax)
             self.notify()
 
 
@@ -591,15 +583,15 @@ class ImageWrapper(notifier.Notifier):
             slices = np.array(slices.T, dtype=np.uint32)
             slices = tuple(it.chain(map(tuple, slices), [(lowVol, highVol)]))
 
-            log.debug('Image {} data written - clearing known data '
-                      'range on volumes {} - {} (write slice: {}; '
-                      'coverage: {}; volRanges: {})'.format(
-                          self.__name,
-                          lowVol,
-                          highVol,
-                          slices,
-                          self.__coverage[:, :, lowVol:highVol],
-                          self.__volRanges[lowVol:highVol, :]))
+            log.debug('Image %s data written - clearing known data '
+                      'range on volumes %s - %s (write slice: %s; '
+                      'coverage: %s; volRanges: %s)',
+                      self.__name,
+                      lowVol,
+                      highVol,
+                      slices,
+                      self.__coverage[:, :, lowVol:highVol],
+                      self.__volRanges[lowVol:highVol, :])
 
             for vol in range(lowVol, highVol):
                 self.__coverage[:, :, vol]    = np.nan
@@ -622,7 +614,7 @@ class ImageWrapper(notifier.Notifier):
         :arg sliceobj: Something which can slice the image data.
         """
 
-        log.debug('Getting image data: {}'.format(sliceobj))
+        log.debug('Getting image data: %s', sliceobj)
 
         shape              = self.__canonicalShape
         realShape          = self.__image.shape
@@ -721,7 +713,7 @@ class ImageWrapper(notifier.Notifier):
         # have any effect.
         self.loadData()
 
-        self.__image.get_data()[sliceobj] = values
+        self.__data[sliceobj] = values
         self.__updateDataRangeOnWrite(slices, values)
 
 
diff --git a/fsl/data/mghimage.py b/fsl/data/mghimage.py
index a143e06ce1e3ce9a2367444c7e67408459e3cd36..be29b5cad5ed56c789f528b478b19d9823093bb9 100644
--- a/fsl/data/mghimage.py
+++ b/fsl/data/mghimage.py
@@ -12,6 +12,7 @@ Freesurfer ``mgh``/``mgz`` image files.
 import os.path as op
 
 import            six
+import numpy   as np
 import nibabel as nib
 
 import fsl.utils.path       as fslpath
@@ -54,7 +55,7 @@ class MGHImage(fslimage.Image):
             name     = 'MGH image'
             filename = None
 
-        data     = image.get_data()
+        data     = np.asanyarray(image.dataobj)
         xform    = image.affine
         vox2surf = image.header.get_vox2ras_tkr()
 
diff --git a/fsl/wrappers/wrapperutils.py b/fsl/wrappers/wrapperutils.py
index c89afb2aadb002e1e6838235578ce124ce4caf41..81fb86826849b6a4eb55c9012f71c0c61cf40cc6 100644
--- a/fsl/wrappers/wrapperutils.py
+++ b/fsl/wrappers/wrapperutils.py
@@ -952,16 +952,17 @@ def fileOrImage(*args, **kwargs):
 
         # create an independent in-memory
         # copy of the image file
-        img = nib.load(path, mmap=False)
+        img  = nib.load(path, mmap=False)
+        data = np.asanyarray(img.dataobj)
 
         # if any arguments were fsl images,
         # that takes precedence.
         if fslimage.Image in intypes:
-            return fslimage.Image(img.get_data(), header=img.header)
+            return fslimage.Image(data, header=img.header)
         # but if all inputs were file names,
         # nibabel takes precedence
         elif nib.nifti1.Nifti1Image in intypes or len(intypes) == 0:
-            return nib.nifti1.Nifti1Image(img.get_data(), None, img.header)
+            return nib.nifti1.Nifti1Image(data, None, img.header)
 
         # this function should not be called
         # under any other circumstances
diff --git a/tests/test_dtifit.py b/tests/test_dtifit.py
index 9e9f80ac07808fc5343ebabae4e3e974e5162225..00fc058ffca6dabc4f70a0e15459748681d0d539 100644
--- a/tests/test_dtifit.py
+++ b/tests/test_dtifit.py
@@ -169,12 +169,12 @@ def test_DTIFitTensor():
         l2file = op.join(testdir, 'dti_L2.nii')
         l3file = op.join(testdir, 'dti_L3.nii')
 
-        v1 = tests.make_random_image(v1file, (5, 5, 5, 3)).get_data()
-        v2 = tests.make_random_image(v2file, (5, 5, 5, 3)).get_data()
-        v3 = tests.make_random_image(v3file, (5, 5, 5, 3)).get_data()
-        l1 = tests.make_random_image(l1file, (5, 5, 5))   .get_data()
-        l2 = tests.make_random_image(l2file, (5, 5, 5))   .get_data()
-        l3 = tests.make_random_image(l3file, (5, 5, 5))   .get_data()
+        v1 = np.asanyarray(tests.make_random_image(v1file, (5, 5, 5, 3)).dataobj)
+        v2 = np.asanyarray(tests.make_random_image(v2file, (5, 5, 5, 3)).dataobj)
+        v3 = np.asanyarray(tests.make_random_image(v3file, (5, 5, 5, 3)).dataobj)
+        l1 = np.asanyarray(tests.make_random_image(l1file, (5, 5, 5))   .dataobj)
+        l2 = np.asanyarray(tests.make_random_image(l2file, (5, 5, 5))   .dataobj)
+        l3 = np.asanyarray(tests.make_random_image(l3file, (5, 5, 5))   .dataobj)
 
         dtiobj = dtifit.DTIFitTensor(testdir)
 
diff --git a/tests/test_ensure.py b/tests/test_ensure.py
index 961377d39b841b52cb6d804deff119021edd526c..a8ce2e68bec571c5f3dde218ade22355e51cf664 100644
--- a/tests/test_ensure.py
+++ b/tests/test_ensure.py
@@ -27,7 +27,7 @@ def test_ensureIsImage():
 
         for l in loaded:
             assert isinstance(l, nib.nifti1.Nifti1Image)
-            assert np.all(img.get_data() == l.get_data())
+            assert np.all(np.asanyarray(img.dataobj) == np.asanyarray(l.dataobj))
 
         l = None
         loaded = None
diff --git a/tests/test_image.py b/tests/test_image.py
index 3bb77d01a3d8f6e7d559fc79796db7097c04615e..e12c6acb7dd2b829a588b13551949e83d4bb9e43 100644
--- a/tests/test_image.py
+++ b/tests/test_image.py
@@ -823,7 +823,7 @@ def _test_Image_changeData(imgtype):
         img.register('name3', onDataRange, 'dataRange')
 
         # Calculate the actual data range
-        data   = img.nibImage.get_data()
+        data   = np.asanyarray(img.nibImage.dataobj)
         dmin   = data.min()
         dmax   = data.max()
         drange = dmax - dmin
@@ -1098,14 +1098,16 @@ def _test_Image_save(imgtype):
             # Load the image back in
             img2 = fslimage.Image(img.dataSource)
 
-            assert img.saveState
-            assert img.dataSource == expDataSource
+            assert img2.saveState
+            assert img2.dataSource == expDataSource
 
-            if imgtype > 0:
-                assert np.all(np.isclose(img.voxToWorldMat, xform))
+            for i in (img, img2):
+                if imgtype > 0:
+                    assert np.all(np.isclose(i.voxToWorldMat, xform))
 
-            for (x, y, z), v in zip(rvoxs, rvals):
-                assert np.isclose(img[x, y, z], v)
+                for (x, y, z), v in zip(rvoxs, rvals):
+                    assert np.isclose(i[x, y, z], v)
+            img  = None
             img2 = None
 
 
@@ -1135,7 +1137,7 @@ def _test_Image_init_xform(imgtype):
         # an image created off a
         # header should have
         # identical sform/qform
-        fimg = fslimage.Image(img.get_data(), header=img.header)
+        fimg = fslimage.Image(np.asanyarray(img.dataobj), header=img.header)
 
         fsform, fsform_code = fimg.header.get_sform(True)
         fqform, fqform_code = fimg.header.get_qform(True)
@@ -1153,7 +1155,7 @@ def _test_Image_init_xform(imgtype):
         # set to that xform,
         # qform to None, and
         # and codes set to (s2, q0)
-        fimg = fslimage.Image(img.get_data(), xform=sform)
+        fimg = fslimage.Image(np.asanyarray(img.dataobj), xform=sform)
 
         fsform, fsform_code = fimg.header.get_sform(True)
         fqform, fqform_code = fimg.header.get_qform(True)
@@ -1174,7 +1176,7 @@ def _test_Image_init_xform(imgtype):
         rxform = affine.compose(np.random.random(3),
                                 np.random.random(3),
                                 np.random.random(3))
-        fimg = fslimage.Image(img.get_data(),
+        fimg = fslimage.Image(np.asanyarray(img.dataobj),
                               header=img.header,
                               xform=rxform)
 
diff --git a/tests/test_imagewrapper.py b/tests/test_imagewrapper.py
index fc78da8460f3f1d7ed20364206d3ef12cbae0934..868245fefd24738197b47c2eab4e3214ffc28837 100644
--- a/tests/test_imagewrapper.py
+++ b/tests/test_imagewrapper.py
@@ -790,7 +790,7 @@ def _test_ImageWrapper_write_out(niters, seed, threaded):
                 wrapper[tuple(sliceobjs)] = newData
                 _ImageWraper_busy_wait(wrapper)
 
-                expLo, expHi = coverageDataRange(img.get_data(), cov, slices)
+                expLo, expHi = coverageDataRange(np.asanyarray(img.dataobj), cov, slices)
                 newLo, newHi = wrapper.dataRange
 
                 # print('Old    range: {} - {}'.format(clo,   chi))
@@ -917,10 +917,10 @@ def _test_ImageWrapper_write_in_overlap(niters, seed, threaded):
 
                 print('Expected range:    {} - {}'.format(expLo, expHi))
                 print('New range:         {} - {}'.format(newLo, newHi))
-                print('Slice min/max:     {} - {}'.format(img.get_data()[tuple(sliceobjs)].min(),
-                                                          img.get_data()[tuple(sliceobjs)].max()))
-                print('Data min/max:      {} - {}'.format(img.get_data().min(),
-                                                          img.get_data().max()))
+                print('Slice min/max:     {} - {}'.format(np.asanyarray(img.dataobj)[tuple(sliceobjs)].min(),
+                                                          np.asanyarray(img.dataobj)[tuple(sliceobjs)].max()))
+                print('Data min/max:      {} - {}'.format(np.asanyarray(img.dataobj).min(),
+                                                          np.asanyarray(img.dataobj).max()))
 
                 assert np.all(newCov == expCov)
 
diff --git a/tests/test_immv_imcp.py b/tests/test_immv_imcp.py
index dd1273d78e162a3ee96638cfba4b12850d1b080f..3e324841a118b16482d314251384188b78ad614b 100644
--- a/tests/test_immv_imcp.py
+++ b/tests/test_immv_imcp.py
@@ -14,6 +14,7 @@ import               os
 import               shutil
 import               tempfile
 
+import numpy as np
 import nibabel as nib
 
 import fsl.utils.imcp        as imcp
@@ -31,7 +32,7 @@ def print(*args, **kwargs):
 
 
 def makeImage(filename):
-    return hash(make_random_image(filename).get_data().tobytes())
+    return hash(np.asanyarray(make_random_image(filename).dataobj).tobytes())
 
 
 def checkImageHash(filename, datahash):
@@ -39,7 +40,7 @@ def checkImageHash(filename, datahash):
     """
 
     img = nib.load(filename)
-    assert hash(img.get_data().tobytes()) == datahash
+    assert hash(np.asanyarray(img.dataobj).tobytes()) == datahash
 
 
 def checkFilesToExpect(files, outdir, outputType, datahashes):
diff --git a/tests/test_mghimage.py b/tests/test_mghimage.py
index f34f635c17759eb9803d34c0405d4a2d93e4a9eb..e2ad7739f4796b41297090d09882d0ef5cac6dfe 100644
--- a/tests/test_mghimage.py
+++ b/tests/test_mghimage.py
@@ -31,7 +31,7 @@ def test_MGHImage():
     v2s      = nbimg.header.get_vox2ras_tkr()
     w2s      = affine.concat(v2s, affine.invert(nbimg.affine))
 
-    assert np.all(np.isclose(img[:],             nbimg.get_data()))
+    assert np.all(np.isclose(img[:],             np.asanyarray(nbimg.dataobj)))
     assert np.all(np.isclose(img.voxToWorldMat,  nbimg.affine))
     assert np.all(np.isclose(img.voxToSurfMat,   v2s))
     assert np.all(np.isclose(img.surfToVoxMat,   affine.invert(v2s)))
@@ -44,7 +44,7 @@ def test_MGHImage():
 
     # Load from an in-memory nibabel object
     img = fslmgh.MGHImage(nbimg)
-    assert np.all(np.isclose(img[:],            nbimg.get_data()))
+    assert np.all(np.isclose(img[:],            np.asanyarray(nbimg.dataobj)))
     assert np.all(np.isclose(img.voxToWorldMat, nbimg.affine))
     assert img.dataSource   is None
     assert img.mghImageFile is None
diff --git a/tests/test_scripts/test_fsl_apply_x5.py b/tests/test_scripts/test_fsl_apply_x5.py
index dbe2bd3a7ec6cfae76fbb07e8a7e663fc5b135ab..49b9908061f20ea678d889dcad4c5fe46e407454 100644
--- a/tests/test_scripts/test_fsl_apply_x5.py
+++ b/tests/test_scripts/test_fsl_apply_x5.py
@@ -44,7 +44,7 @@ def test_help():
 
 
 
-def test_linear():
+def test_linear(seed):
     with tempdir.tempdir():
 
         src2ref = _random_affine()
@@ -65,7 +65,7 @@ def test_linear():
         assert np.all(np.isclose(result.data, expect))
 
 
-def test_nonlinear():
+def test_nonlinear(seed):
     with tempdir.tempdir():
 
         src2ref = _random_affine()
@@ -92,7 +92,7 @@ def test_nonlinear():
         assert np.all(np.isclose(result, expect))
 
 
-def test_linear_altref():
+def test_linear_altref(seed):
     with tempdir.tempdir():
 
         src2ref = affine.scaleOffsetXform([1, 1, 1], [5,  5,  5])
@@ -119,7 +119,7 @@ def test_linear_altref():
         assert np.all(result.data == expect)
 
 
-def test_nonlinear_altref():
+def test_nonlinear_altref(seed):
     with tempdir.tempdir():
 
         src2ref = affine.scaleOffsetXform([1, 1, 1], [5,  5,  5])
@@ -149,7 +149,7 @@ def test_nonlinear_altref():
         assert np.all(result.data == expect)
 
 
-def test_linear_altsrc():
+def test_linear_altsrc(seed):
     with tempdir.tempdir():
 
         src2ref = _random_affine()
@@ -198,7 +198,7 @@ def test_linear_altsrc():
         assert np.all(np.isclose(outoff.data, expoff))
 
 
-def test_nonlinear_altsrc():
+def test_nonlinear_altsrc(seed):
     with tempdir.tempdir():
 
         src2ref = _random_affine()
@@ -252,13 +252,20 @@ def test_nonlinear_altsrc():
         assert outhi .sameSpace(ref)
         assert outoff.sameSpace(ref)
 
-        # We get boundary issues just at the first
-        # voxel, so I'm masking that voxel out
-        for img in (out, outlo, outhi, outoff,
-                    exp, explo, exphi, expoff):
-            img[0, 0, 0] = 0
-
-        assert np.all(np.isclose(out   .data, exp))
-        assert np.all(np.isclose(outlo .data, explo))
-        assert np.all(np.isclose(outhi .data, exphi))
-        assert np.all(np.isclose(outoff.data, expoff))
+        # We get boundary cropping,
+        # so ignore edge slices
+        out    = out   .data[1:-1, 1:-1, 1:-1]
+        outlo  = outlo .data[1:-1, 1:-1, 1:-1]
+        outhi  = outhi .data[1:-1, 1:-1, 1:-1]
+        outoff = outoff.data[1:-1, 1:-1, 1:-1]
+        exp    = exp[        1:-1, 1:-1, 1:-1]
+        explo  = explo[      1:-1, 1:-1, 1:-1]
+        exphi  = exphi[      1:-1, 1:-1, 1:-1]
+        expoff = expoff[     1:-1, 1:-1, 1:-1]
+
+        tol = dict(atol=1e-3, rtol=1e-3)
+
+        assert np.all(np.isclose(out,    exp,   **tol))
+        assert np.all(np.isclose(outlo,  explo, **tol))
+        assert np.all(np.isclose(outhi,  exphi, **tol))
+        assert np.all(np.isclose(outoff, expoff, **tol))
diff --git a/tests/test_scripts/test_immv_imcp.py b/tests/test_scripts/test_immv_imcp.py
index 607a7da5d761d69c824f5d85f27751b8bf5a395d..582a23d7c74f3ac53025fcafa51e579d33ff3cc4 100644
--- a/tests/test_scripts/test_immv_imcp.py
+++ b/tests/test_scripts/test_immv_imcp.py
@@ -11,6 +11,7 @@ import               tempfile
 
 import               pytest
 
+import numpy as np
 import nibabel as nib
 
 from   fsl.utils.tempdir import tempdir
@@ -269,7 +270,7 @@ def test_imcp_script_shouldPass(move=False):
                         for inf in infiles:
                             img     = nib.load(op.join(tindir, inf),
                                                mmap=False)
-                            imghash = hash(img.get_data().tobytes())
+                            imghash = hash(np.asanyarray(img.dataobj).tobytes())
                             img = None
                             imageHashes.append(imghash)
 
diff --git a/tests/test_wrapperutils.py b/tests/test_wrapperutils.py
index 601e9f8ffa25a423d05ae26360708f11ad7eab88..9e1d5d6a82c56801a60918801c28bfeb81f4e825 100644
--- a/tests/test_wrapperutils.py
+++ b/tests/test_wrapperutils.py
@@ -231,8 +231,8 @@ def test_fileOrImage():
 
     @wutils.fileOrImage('img1', 'img2', 'output')
     def func(img1, **kwargs):
-        img1   = nib.load(img1).get_data()
-        img2   = nib.load(kwargs['img2']).get_data()
+        img1   = np.asanyarray(nib.load(img1).dataobj)
+        img2   = np.asanyarray(nib.load(kwargs['img2']).dataobj)
         output = nib.nifti1.Nifti1Image(img1 * img2, np.eye(4))
         nib.save(output, kwargs['output'])
 
@@ -247,43 +247,43 @@ def test_fileOrImage():
 
         # file  file  file
         func('img1.nii', img2='img2.nii', output='output.nii')
-        assert np.all(nib.load('output.nii').get_data() == expected)
+        assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
         os.remove('output.nii')
 
         # file  file  array
         result = func('img1.nii', img2='img2.nii', output=wutils.LOAD)['output']
-        assert np.all(result.get_data() == expected)
+        assert np.all(np.asanyarray(result.dataobj) == expected)
 
         # file  array file
         func('img1.nii', img2=img2, output='output.nii')
-        assert np.all(nib.load('output.nii').get_data() == expected)
+        assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
         os.remove('output.nii')
 
         # file  array array
         result = func('img1.nii', img2=img2, output=wutils.LOAD)['output']
-        assert np.all(result.get_data() == expected)
+        assert np.all(np.asanyarray(result.dataobj) == expected)
 
         # array file  file
         func(img1, img2='img2.nii', output='output.nii')
-        assert np.all(nib.load('output.nii').get_data() == expected)
+        assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
         os.remove('output.nii')
 
         # array file  array
         result = func(img1, img2='img2.nii', output=wutils.LOAD)['output']
-        assert np.all(result.get_data() == expected)
+        assert np.all(np.asanyarray(result.dataobj) == expected)
 
         # array array file
         func(img1, img2=img2, output='output.nii')
-        assert np.all(nib.load('output.nii').get_data() == expected)
+        assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
         os.remove('output.nii')
 
         # array array array
         result = func(img1, img2=img2, output=wutils.LOAD)['output']
-        assert np.all(result.get_data() == expected)
+        assert np.all(np.asanyarray(result.dataobj) == expected)
 
         # in-memory image, file, file
         result = func(img3, img2='img2.nii', output='output.nii')
-        assert np.all(nib.load('output.nii').get_data() == expected)
+        assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
         os.remove('output.nii')
 
         # fslimage, file, load
@@ -307,7 +307,7 @@ def test_fileOrImage():
         # nib.image, nib.image, load
         result = func(img1, img2=img2, output=wutils.LOAD)['output']
         assert isinstance(result, nib.nifti1.Nifti1Image)
-        assert np.all(result.get_data()[:] == expected)
+        assert np.all(np.asanyarray(result.dataobj)[:] == expected)
 
 
 def test_fileOrThing_sequence():
@@ -357,7 +357,7 @@ def test_fileOrThing_outprefix():
 
     @wutils.fileOrImage('img', outprefix='output_base')
     def basefunc(img, output_base):
-        img = nib.load(img).get_data()
+        img = np.asanyarray(nib.load(img).dataobj)
 
         out1 = nib.nifti1.Nifti1Image(img * 5,  np.eye(4))
         out2 = nib.nifti1.Nifti1Image(img * 10, np.eye(4))
@@ -368,31 +368,31 @@ def test_fileOrThing_outprefix():
 
     with tempdir.tempdir() as td:
         img  = nib.nifti1.Nifti1Image(np.array([[1, 2], [3, 4]]), np.eye(4))
-        exp1 = img.get_data() * 5
-        exp2 = img.get_data() * 10
+        exp1 = np.asanyarray(img.dataobj) * 5
+        exp2 = np.asanyarray(img.dataobj) * 10
         nib.save(img, 'img.nii')
 
         basefunc('img.nii', 'myout')
-        assert np.all(nib.load('myout_times5.nii.gz') .get_data() == exp1)
-        assert np.all(nib.load('myout_times10.nii.gz').get_data() == exp2)
+        assert np.all(np.asanyarray(nib.load('myout_times5.nii.gz') .dataobj) == exp1)
+        assert np.all(np.asanyarray(nib.load('myout_times10.nii.gz').dataobj) == exp2)
         cleardir(td, 'myout*')
 
         basefunc(img, 'myout')
-        assert np.all(nib.load('myout_times5.nii.gz') .get_data() == exp1)
-        assert np.all(nib.load('myout_times10.nii.gz').get_data() == exp2)
+        assert np.all(np.asanyarray(nib.load('myout_times5.nii.gz') .dataobj) == exp1)
+        assert np.all(np.asanyarray(nib.load('myout_times10.nii.gz').dataobj) == exp2)
         cleardir(td, 'myout*')
 
         res = basefunc(img, 'myout', myout_times5=wutils.LOAD)
-        assert np.all(res['myout_times5'].get_data() == exp1)
+        assert np.all(np.asanyarray(res['myout_times5'].dataobj) == exp1)
         cleardir(td, 'myout*')
 
         res = basefunc(img, 'myout', myout_times10=wutils.LOAD)
-        assert np.all(res['myout_times10'].get_data() == exp2)
+        assert np.all(np.asanyarray(res['myout_times10'].dataobj) == exp2)
         cleardir(td, 'myout*')
 
         res = basefunc(img, 'myout', myout=wutils.LOAD)
-        assert np.all(res['myout_times5'] .get_data() == exp1)
-        assert np.all(res['myout_times10'].get_data() == exp2)
+        assert np.all(np.asanyarray(res['myout_times5'] .dataobj) == exp1)
+        assert np.all(np.asanyarray(res['myout_times10'].dataobj) == exp2)
         cleardir(td, 'myout*')
 
 
@@ -402,7 +402,7 @@ def test_fileOrThing_outprefix_differentTypes():
     def func(img, outpref):
 
         img  = nib.load(img)
-        img  = nib.nifti1.Nifti1Image(img.get_data() * 2, np.eye(4))
+        img  = nib.nifti1.Nifti1Image(np.asanyarray(img.dataobj) * 2, np.eye(4))
         text = '1234567890'
 
         nib.save(img, '{}_image.nii.gz' .format(outpref))
@@ -412,23 +412,23 @@ def test_fileOrThing_outprefix_differentTypes():
 
     with tempdir.tempdir() as td:
         img  = nib.nifti1.Nifti1Image(np.array([[1, 2], [3, 4]]), np.eye(4))
-        expi = img.get_data() * 2
+        expi = np.asanyarray(img.dataobj) * 2
         expt = '1234567890'
 
         func(img, 'myout')
-        assert np.all(nib.load('myout_image.nii.gz') .get_data() == expi)
+        assert np.all(np.asanyarray(nib.load('myout_image.nii.gz') .dataobj) == expi)
         with open('myout_text.txt', 'rt') as f:
             assert f.read().strip() == expt
         cleardir(td, 'myout*')
 
         res = func(img, 'myout', myout_image=wutils.LOAD)
         assert list(res.keys()) == ['myout_image']
-        assert np.all(res['myout_image'].get_data() == expi)
+        assert np.all(np.asanyarray(res['myout_image'].dataobj) == expi)
         cleardir(td, 'myout*')
 
         res = func(img, 'myout', myout=wutils.LOAD)
         assert list(res.keys()) == ['myout_image']
-        assert np.all(res['myout_image'].get_data() == expi)
+        assert np.all(np.asanyarray(res['myout_image'].dataobj) == expi)
         cleardir(td, 'myout*')
 
         res = func(img, 'myout', myout_text=wutils.LOAD)
@@ -445,8 +445,8 @@ def test_fileOrThing_outprefix_directory():
     @wutils.fileOrImage('img', outprefix='outpref')
     def func(img, outpref):
         img  = nib.load(img)
-        img2 = nib.nifti1.Nifti1Image(img.get_data() * 2, np.eye(4))
-        img4 = nib.nifti1.Nifti1Image(img.get_data() * 4, np.eye(4))
+        img2 = nib.nifti1.Nifti1Image(np.asanyarray(img.dataobj) * 2, np.eye(4))
+        img4 = nib.nifti1.Nifti1Image(np.asanyarray(img.dataobj) * 4, np.eye(4))
 
         outdir = op.abspath('{}_imgs'.format(outpref))
 
@@ -457,8 +457,8 @@ def test_fileOrThing_outprefix_directory():
 
     with tempdir.tempdir() as td:
         img  = nib.nifti1.Nifti1Image(np.array([[1, 2], [3, 4]]), np.eye(4))
-        exp2 = img.get_data() * 2
-        exp4 = img.get_data() * 4
+        exp2 = np.asanyarray(img.dataobj) * 2
+        exp4 = np.asanyarray(img.dataobj) * 4
 
         res = func(img, 'myout')
         assert len(res) == 0
@@ -469,17 +469,17 @@ def test_fileOrThing_outprefix_directory():
 
         res = func(img, 'myout', myout_imgs=wutils.LOAD)
         assert len(res) == 2
-        assert np.all(res[op.join('myout_imgs', 'img2')].get_data() == exp2)
-        assert np.all(res[op.join('myout_imgs', 'img4')].get_data() == exp4)
+        assert np.all(np.asanyarray(res[op.join('myout_imgs', 'img2')].dataobj) == exp2)
+        assert np.all(np.asanyarray(res[op.join('myout_imgs', 'img4')].dataobj) == exp4)
 
         res = func(img, 'myout', **{op.join('myout_imgs', 'img2') : wutils.LOAD})
         assert len(res) == 1
-        assert np.all(res[op.join('myout_imgs', 'img2')].get_data() == exp2)
+        assert np.all(np.asanyarray(res[op.join('myout_imgs', 'img2')].dataobj) == exp2)
 
         res = func(img, 'myout', **{op.join('myout_imgs', 'img') : wutils.LOAD})
         assert len(res) == 2
-        assert np.all(res[op.join('myout_imgs', 'img2')].get_data() == exp2)
-        assert np.all(res[op.join('myout_imgs', 'img4')].get_data() == exp4)
+        assert np.all(np.asanyarray(res[op.join('myout_imgs', 'img2')].dataobj) == exp2)
+        assert np.all(np.asanyarray(res[op.join('myout_imgs', 'img4')].dataobj) == exp4)
 
         os.mkdir('foo')
         res = func(img, op.join('foo', 'myout'))
@@ -492,8 +492,8 @@ def test_fileOrThing_outprefix_directory():
         os.mkdir('foo')
         res = func(img, op.join('foo', 'myout'), **{op.join('foo', 'myout') : wutils.LOAD})
         assert len(res) == 2
-        assert np.all(res[op.join('foo', 'myout_imgs', 'img2')].get_data() == exp2)
-        assert np.all(res[op.join('foo', 'myout_imgs', 'img4')].get_data() == exp4)
+        assert np.all(np.asanyarray(res[op.join('foo', 'myout_imgs', 'img2')].dataobj) == exp2)
+        assert np.all(np.asanyarray(res[op.join('foo', 'myout_imgs', 'img4')].dataobj) == exp4)
 
 
 def test_chained_fileOrImageAndArray():
@@ -503,7 +503,7 @@ def test_chained_fileOrImageAndArray():
         image = nib.load(image)
         array = np.loadtxt(array)
 
-        outimg = nib.nifti1.Nifti1Image(image.get_data() * 2, np.eye(4))
+        outimg = nib.nifti1.Nifti1Image(np.asanyarray(image.dataobj) * 2, np.eye(4))
 
         np.savetxt(outarray, array * 2)
         outimg.to_filename(outimage)
@@ -511,7 +511,7 @@ def test_chained_fileOrImageAndArray():
     image = nib.nifti1.Nifti1Image(np.array([[1,  2], [ 3,  4]]), np.eye(4))
     array = np.array([[5, 6, 7, 8]])
 
-    expimg = nib.nifti1.Nifti1Image(image.get_data() * 2, np.eye(4))
+    expimg = nib.nifti1.Nifti1Image(np.asanyarray(image.dataobj) * 2, np.eye(4))
     exparr = array * 2
 
     with tempdir.tempdir():
@@ -520,31 +520,31 @@ def test_chained_fileOrImageAndArray():
         np.savetxt('array.txt', array)
 
         func('image.nii', 'array.txt', 'outimg.nii', 'outarr.txt')
-        assert np.all(nib.load('outimg.nii').get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(nib.load('outimg.nii').dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(np.loadtxt('outarr.txt') == exparr)
 
         func('image.nii', array, 'outimg.nii', 'outarr.txt')
-        assert np.all(nib.load('outimg.nii').get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(nib.load('outimg.nii').dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(np.loadtxt('outarr.txt') == exparr)
 
         func( image, 'array.txt', 'outimg.nii', 'outarr.txt')
-        assert np.all(nib.load('outimg.nii').get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(nib.load('outimg.nii').dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(np.loadtxt('outarr.txt') == exparr)
 
         func( image, array, 'outimg.nii', 'outarr.txt')
-        assert np.all(nib.load('outimg.nii').get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(nib.load('outimg.nii').dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(np.loadtxt('outarr.txt') == exparr)
 
         res = func(image, array, wutils.LOAD, 'outarr.txt')
-        assert np.all(res['outimage'].get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(res['outimage'].dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(np.loadtxt('outarr.txt') == exparr)
 
         res = func(image, array, 'outimg.nii', wutils.LOAD)
-        assert np.all(nib.load('outimg.nii').get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(nib.load('outimg.nii').dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(res['outarray'] == exparr)
 
         res = func(image, array, wutils.LOAD, wutils.LOAD)
-        assert np.all(res['outimage'].get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(res['outimage'].dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(res['outarray'] == exparr)
 
 
@@ -561,7 +561,7 @@ def test_fileOrThing_chained_outprefix():
         image = nib.load(image)
         array = np.loadtxt(array)
 
-        outimg = nib.nifti1.Nifti1Image(image.get_data() * 2, np.eye(4))
+        outimg = nib.nifti1.Nifti1Image(np.asanyarray(image.dataobj) * 2, np.eye(4))
         outarr = array * 2
 
         np.savetxt('{}_array.txt'.format(out), outarr)
@@ -570,17 +570,17 @@ def test_fileOrThing_chained_outprefix():
     image = nib.nifti1.Nifti1Image(np.array([[1,  2], [ 3,  4]]), np.eye(4))
     array = np.array([[5, 6, 7, 8]])
 
-    expimg = nib.nifti1.Nifti1Image(image.get_data() * 2, np.eye(4))
+    expimg = nib.nifti1.Nifti1Image(np.asanyarray(image.dataobj) * 2, np.eye(4))
     exparr = array * 2
 
     with tempdir.tempdir():
 
         func(image, array, 'myout')
-        assert np.all(nib.load('myout_image.nii').get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(nib.load('myout_image.nii').dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(np.loadtxt('myout_array.txt') == exparr)
 
         res = func(image, array, wutils.LOAD)
-        assert np.all(res['out_image'].get_data() == expimg.get_data())
+        assert np.all(np.asanyarray(res['out_image'].dataobj) == np.asanyarray(expimg.dataobj))
         assert np.all(res['out_array'] == exparr)