diff --git a/requirements.conda.txt b/requirements.conda.txt
index 4ddfdf46f..686d2837c 100644
--- a/requirements.conda.txt
+++ b/requirements.conda.txt
@@ -1,4 +1,5 @@
 drmaa
+hdf5plugin
 ispyb>=11.0.1
 junit-xml>=1.9
 marshmallow-sqlalchemy
diff --git a/src/dlstbx/services/xray_centering.py b/src/dlstbx/services/xray_centering.py
index 68ec99740..735b19aa9 100644
--- a/src/dlstbx/services/xray_centering.py
+++ b/src/dlstbx/services/xray_centering.py
@@ -177,7 +177,19 @@ def garbage_collect(self):
                     del self._centering_data[dcid]
 
     def add_pia_result(self, rw, header, message):
-        """Process incoming PIA result."""
+        """Process incoming PIA result.
+
+        This is invoked repeatedly as each image in the gridscan comes through.
+
+        Recipe Step Inputs:
+            gridinfo: This is the ISPyB DataCollectionGridInfo for the DataCollection
+            parameters: Other parameters in the recipe - see the relevant per-image-analysis-xxx.json
+        Args:
+            rw: an instance of RecipeWrapper
+            header: The message header
+            message: A dict which decodes to an instance of Message
+
+        """
 
         try:
             recipe_step = RecipeStep(**rw.recipe_step)
diff --git a/src/dlstbx/util/xray_centering.py b/src/dlstbx/util/xray_centering.py
index 5caf03925..4218e341f 100644
--- a/src/dlstbx/util/xray_centering.py
+++ b/src/dlstbx/util/xray_centering.py
@@ -45,6 +45,11 @@ class GridScan2DResult(GridScanResultBase):
 def reshape_grid(
     data: np.ndarray, steps: tuple[int, int], snaked: bool, orientation: Orientation
 ) -> np.ndarray:
+    """
+    Converts the linear array to column major 2D, and unsnake.
+    NOTE, the current implementation deliberately mutates the input!
+    """
+    # Transpose the array if orientation is horizontal
     if orientation == Orientation.VERTICAL:
         data = data.reshape(steps)
     else:
diff --git a/src/dlstbx/util/xray_centering_3d.py b/src/dlstbx/util/xray_centering_3d.py
index 5cca59035..62549e8fc 100644
--- a/src/dlstbx/util/xray_centering_3d.py
+++ b/src/dlstbx/util/xray_centering_3d.py
@@ -17,6 +17,24 @@
 
 
 class GridScan3DResult(GridScanResultBase):
+    """
+        Represents a single gridscan result, corresponding to a diffracting centre.
+
+        Coordinates expressed are in terms of grid boxes, with the centre of the box
+        lying on half-integer coordinates.
+
+        Attributes:
+             centre_of_mass: The position of the centre of mass of the crystal, for a
+             crystal of size (1, 1, 1) this will be on half-integer coordinates
+             max_voxel: Position of the maximum voxel, on integer coordinates!!!
+             max_count: max count achieved in a single voxel for the crystal
+             n_voxels: Number of voxels in the diffracting centre
+             total_count: Total of above-threshold spot counts in the labelled voxels
+             bounding_box: The rectangular prism that bounds the crystal, expressed
+                as the volume of whole boxes as a half-open range i.e such that
+                p1 = (x1, y1, z1) <= p < p2 = (x2, y2, z2) and
+                p2 - p1 gives the dimensions in whole voxels.
+    """
     centre_of_mass: Tuple[float, ...]
     max_voxel: Tuple[int, ...]
     max_count: float
diff --git a/tests/util/test_xray_centering.py b/tests/util/test_xray_centering.py
index d4df29481..dc0fc041f 100644
--- a/tests/util/test_xray_centering.py
+++ b/tests/util/test_xray_centering.py
@@ -315,3 +315,63 @@ def test_single_connected_region(data, reflections_in_best_image):
     )
     assert result.centre_x == result.centre_x_box == 5
     assert result.centre_y == result.centre_y_box == 5
+
+
+EXPECTED_OUTPUT_COL_MAJOR = (np.array([
+    [1, 5, 9],
+    [2, 6, 10],
+    [3, 7, 11],
+    [4, 8, 12]
+]))
+
+GRID_INPUT_ROW_MAJOR_SNAKED_L_TO_R_FIRST = (np.array([
+    1, 2, 3, 4,
+    8, 7, 6, 5,
+    9, 10, 11, 12,
+]))
+GRID_INPUT_ROW_MAJOR_SNAKED_R_TO_L_FIRST = (np.array([
+    4, 3, 2, 1,
+    5, 6, 7, 8,
+    12, 11, 10, 9,
+]))
+GRID_INPUT_ROW_MAJOR_NOT_SNAKED_L_TO_R = (np.array([
+    1, 2, 3, 4,
+    5, 6, 7, 8,
+    9, 10, 11, 12,
+]))
+GRID_INPUT_COL_MAJOR_NOT_SNAKED = (np.array([
+    1, 5, 9,
+    2, 6, 10,
+    3, 7, 11,
+    4, 8, 12
+]))
+GRID_INPUT_COL_MAJOR_SNAKED_T_TO_B_FIRST = (np.array([
+    1, 5, 9,
+    10, 6, 2,
+    3, 7, 11,
+    12, 8, 4
+ ]))
+
+@pytest.mark.parametrize("data_in, expected_data, steps, snaked, orientation",
+                         [
+                             [GRID_INPUT_ROW_MAJOR_SNAKED_L_TO_R_FIRST, EXPECTED_OUTPUT_COL_MAJOR, (4, 3), True, dlstbx.util.xray_centering.Orientation.HORIZONTAL],
+                             [GRID_INPUT_ROW_MAJOR_NOT_SNAKED_L_TO_R, EXPECTED_OUTPUT_COL_MAJOR, (4, 3), False, 
+                              dlstbx.util.xray_centering.Orientation.HORIZONTAL],
+                             # -ve direction first snaking causes the axis to be flipped,
+                             # is there anything that uses/relies on this?
+                             # [GRID_INPUT_ROW_MAJOR_SNAKED_R_TO_L_FIRST, EXPECTED_OUTPUT_COL_MAJOR, (4, 3), True, 
+                             # dlstbx.util.xray_centering.Orientation.HORIZONTAL],
+                             [GRID_INPUT_COL_MAJOR_NOT_SNAKED, EXPECTED_OUTPUT_COL_MAJOR, (4, 3), False,
+                             dlstbx.util.xray_centering.Orientation.VERTICAL],
+                             [GRID_INPUT_COL_MAJOR_SNAKED_T_TO_B_FIRST, EXPECTED_OUTPUT_COL_MAJOR, (4, 3), True,
+                             dlstbx.util.xray_centering.Orientation.VERTICAL],
+                          ]
+                         )
+def test_reshape_grid(
+        data_in, expected_data, steps, snaked, orientation
+):
+    # old_data_in = data_in.copy()
+    data_out = dlstbx.util.xray_centering.reshape_grid(data_in, steps, snaked=snaked, orientation=orientation)
+    assert np.all(data_out == expected_data), f"{data_out} != {expected_data}"
+    # The current operation of the gridscan processing relies on this mutation of the input
+    # assert np.all(data_in == old_data_in), f"{data_in} != {old_data_in}"
diff --git a/tests/util/test_xray_centering_3d.py b/tests/util/test_xray_centering_3d.py
index 4922bbb4f..37c0843ba 100644
--- a/tests/util/test_xray_centering_3d.py
+++ b/tests/util/test_xray_centering_3d.py
@@ -50,3 +50,36 @@ def test_gridscan3d():
         "total_count": 44128.0,
         "bounding_box": ((2, 3, 2), (7, 6, 6)),
     }
+
+#  fmt: off
+@pytest.mark.parametrize("input_data, expected_results",
+                         [
+                             [  # NB x is down
+                                 (  # xy
+                                     np.array([[0, 0, 0],
+                                      [0, 0, 0],
+                                      [0, 1, 0],
+                                      [0, 0, 0],]),
+                                    # xz
+                                     np.array([[0, 0, 0],
+                                      [0, 0, 0],
+                                      [1, 0, 0],
+                                      [0, 0, 0]])
+                                 ),
+                                 [dlstbx.util.xray_centering_3d.GridScan3DResult(
+                                     centre_of_mass=(2.5, 1.5, 0.5),
+                                     max_voxel=(2.5, 1.5, 0.5),
+                                     max_count=1,
+                                     n_voxels=1,
+                                     total_count=1,
+                                     bounding_box=((2, 1, 0), (3, 2, 1))
+                                 )],
+                             ],
+                         ])
+#  fmt: on
+def test_gridscan_3d_coordinates(input_data,
+                                 expected_results):
+    results = dlstbx.util.xray_centering_3d.gridscan3d(
+        input_data
+    )
+    assert all([r == e for r, e in zip(results, expected_results, strict=True)]), f"{results} != {expected_results}"