diff --git a/fsl/wrappers/__init__.py b/fsl/wrappers/__init__.py
index 87f410cb9c5475ab4542e96c074d0c04327763e1..25f5ff9f11112426a37107fb1e8712840dfd8a03 100644
--- a/fsl/wrappers/__init__.py
+++ b/fsl/wrappers/__init__.py
@@ -107,4 +107,5 @@ from .melodic      import (melodic,         # noqa
 from .misc         import (fslreorient2std, # noqa
                            fslroi,
                            slicer,
-                           cluster)
+                           cluster,
+                           gps)
diff --git a/fsl/wrappers/misc.py b/fsl/wrappers/misc.py
index 9c21888282729dcdecf554281dec1c2a92173849..9250814830dca73288178360421b383b723e9183 100644
--- a/fsl/wrappers/misc.py
+++ b/fsl/wrappers/misc.py
@@ -85,3 +85,22 @@ def cluster(input, thresh, **kwargs):
     cmd += wutils.applyArgStyle('--=', valmap=valmap, **kwargs)
 
     return cmd
+
+@wutils.fileOrArray('out', 'init')
+@wutils.fslwrapper
+def gps(out, ndir, **kwargs):
+    """Wrapper of the ``gps`` command
+
+    Usage example to get 128 gradient orientations on the whole sphere::
+
+        from fsl.wrappers import gps, LOAD
+        bvecs = gps(LOAD, 128, optws=True)['out']
+    """
+    valmap = {name: wutils.SHOW_IF_TRUE for name in [
+        'optws', 'report', 'verbose'
+    ]}
+
+    cmd = ['gps', f'--ndir={ndir}', f'--out={out}']
+    cmd += wutils.applyArgStyle('--=', valmap=valmap, **kwargs)
+
+    return cmd
diff --git a/tests/test_wrappers/test_wrappers.py b/tests/test_wrappers/test_wrappers.py
index 695bedfb30e585a6f2dd0c966a2165b298f99d77..a75cf27346f43062e50855604690bf02e3822e36 100644
--- a/tests/test_wrappers/test_wrappers.py
+++ b/tests/test_wrappers/test_wrappers.py
@@ -336,7 +336,6 @@ def test_fast():
         assert result.stdout[0] == ' '.join(expected)
 
 
-
 def test_fsl_anat():
     with asrt.disabled(), \
          run.dryrun(), \
@@ -349,3 +348,12 @@ def test_fsl_anat():
                     '-s', '25']
 
         assert result.stdout[0] == ' '.join(expected)
+
+
+def test_gps():
+    with asrt.disabled(), run.dryrun(), mockFSLDIR(bin=('gps',)) as fsldir:
+        gps = op.join(fsldir, 'bin', 'gps')
+        result = fw.gps('bvecs', 128, optws=True, ranseed=123)
+        expected = (gps + ' --ndir=128 --out=bvecs',
+                    ('--optws', '--ranseed=123'))
+        assert checkResult(result.stdout[0], *expected)