introduce assume_volumetric to interpret MRC files as densities (#149)

* fixes #150
* new MRC kwarg assume_volumetric to force parser to interpret data as volumetric, even if space group is set to 0
  (work-around for data written by programs that do not fully follow the MRC specs)
* make assume_volumetric available from Grid() and Grid.load()
* update AUTHORS
* update CHANGELOG
* add tests

---------

Co-authored-by: Oliver Beckstein <orbeckst@gmail.com>

Author: PlethoraChutney

AUTHORS

@@ -27,3 +27,4 @@ Contributors:
 * Zhiyi Wu <xiki-tempula>
 * Olivier Languin-Cattoën <ollyfutur>
 * Andrés Montoya <conradolandia> (logo)
+* Rich Waldo <plethorachutney>

CHANGELOG

@@ -28,6 +28,8 @@ The rules for this file:
   * `Grid` now accepts binary operations with any operand that can be 
     broadcasted to the grid's shape according to `numpy` broadcasting rules
     (PR #142)
+  * `Grid` now allows forcing MRC/CCP4 maps to be read as volumes even when
+    the header indicates they are stacks of 2D images. (#150, PR #149)
 
   Fixes
 

gridData/core.py

@@ -93,6 +93,15 @@ class Grid(object):
       format fails. The default is ``None`` and normally the file
       format is guessed from the file extension.
 
+    assume_volumetric : bool (optional)
+      If ``False`` (default), check the file header to determine whether
+      the data in `grid` is a 3D volume. If ``True``, assume `grid` is volumetric.
+
+      .. Note:: `assume_volumetric` only has an effect when loading
+         MRC/CCP4 files. See :class:`gridData.mrc.MRC`
+
+      .. versionadded:: 1.1.0
+
     Raises
     ------
     TypeError
@@ -196,7 +205,7 @@ class Grid(object):
 
     def __init__(self, grid=None, edges=None, origin=None, delta=None,
                  metadata=None, interpolation_spline_order=3,
-                 file_format=None):
+                 file_format=None, assume_volumetric=False):
         # file formats are guessed from extension == lower case key
         self._exporters = {
             'DX': self._export_dx,
@@ -238,7 +247,7 @@ def __init__(self, grid=None, edges=None, origin=None, delta=None,
                     filename = str(grid)
 
             if filename is not None:
-                self.load(filename, file_format=file_format)
+                self.load(filename, file_format=file_format, assume_volumetric=assume_volumetric)
             else:
                 self._load(grid, edges, metadata, origin, delta)
 
@@ -532,7 +541,8 @@ def _load(
                 "grid={0} edges={1} origin={2} delta={3}".format(
                     grid, edges, origin, delta))
 
-    def load(self, filename, file_format=None):
+    # NOTE: keep loader kwargs in sync between load() and __init__()
+    def load(self, filename, file_format=None, assume_volumetric=False):
         """Load saved grid and edges from `filename`
 
         The :meth:`load` method calls the class's constructor method and
@@ -545,32 +555,32 @@ def load(self, filename, file_format=None):
             # are not really a file
             raise IOError(errno.ENOENT, "file not found", filename)
         loader = self._get_loader(filename, file_format=file_format)
-        loader(filename)
+        loader(filename, assume_volumetric=assume_volumetric)
 
-    def _load_python(self, filename):
+    def _load_python(self, filename, **kwargs):
         with open(filename, 'rb') as f:
             saved = pickle.load(f)
         self._load(grid=saved['grid'],
                    edges=saved['edges'],
                    metadata=saved['metadata'])
 
-    def _load_mrc(self, filename):
+    def _load_mrc(self, filename, assume_volumetric=False, **kwargs):
         """Initializes Grid from a MRC/CCP4 file."""
-        mrcfile = mrc.MRC(filename)
+        mrcfile = mrc.MRC(filename, assume_volumetric=assume_volumetric)
         grid, edges = mrcfile.histogramdd()
         self._load(grid=grid, edges=edges, metadata=self.metadata)
         # Store header for access from Grid object (undocumented)
         # https://github.com/MDAnalysis/GridDataFormats/pull/100#discussion_r782604833
         self._mrc_header = mrcfile.header.copy()
 
-    def _load_dx(self, filename):
+    def _load_dx(self, filename, **kwargs):
         """Initializes Grid from a OpenDX file."""
         dx = OpenDX.field(0)
         dx.read(filename)
         grid, edges = dx.histogramdd()
         self._load(grid=grid, edges=edges, metadata=self.metadata)
 
-    def _load_plt(self, filename):
+    def _load_plt(self, filename, **kwargs):
         """Initialize Grid from gOpenMol plt file."""
         g = gOpenMol.Plt()
         g.read(filename)

gridData/mrc.py

@@ -43,6 +43,12 @@ class MRC(object):
     ----------
     filename : str (optional)
        input file (or stream), can be compressed
+    assume_volumetric : bool (optional)
+      If ``False`` (default), check the file header to determine whether
+      the data in `grid` is a 3D volume. If ``True``, assume `grid` is volumetric.
+
+      .. versionadded:: 1.1.0
+        
 
     Raises
     ------
@@ -86,17 +92,23 @@ class MRC(object):
 
     """
 
-    def __init__(self, filename=None):
+    def __init__(self, filename=None, assume_volumetric=False):
         self.filename = filename
         if filename is not None:
-            self.read(filename)
+            self.read(filename, assume_volumetric=assume_volumetric)
 
-    def read(self, filename):
+    def read(self, filename, assume_volumetric=False):
         """Populate the instance from the MRC/CCP4 file *filename*."""
         if filename is not None:
             self.filename = filename
         with mrcfile.open(filename) as mrc:
-            if not mrc.is_volume():                           #pragma: no cover
+            if assume_volumetric:
+                # non 3D volumes should always fail, regardless of assume_volumetric value
+                is_volume = mrc.data is not None and len(mrc.data.shape) == 3
+            else:
+                is_volume = mrc.is_volume()
+
+            if not is_volume:
                 raise ValueError(
                     "MRC file {} is not a volumetric density.".format(filename))
             self.header = h = mrc.header.copy()

gridData/tests/datafiles/__init__.py

@@ -7,6 +7,7 @@
 DX = importlib_resources.files(__name__) / 'test.dx'
 DXGZ = importlib_resources.files(__name__) / 'test.dx.gz'
 CCP4 = importlib_resources.files(__name__) / 'test.ccp4'
+ISPG_0 = importlib_resources.files(__name__) / "ispg_0.mrc"
 # from http://www.ebi.ac.uk/pdbe/coordinates/files/1jzv.ccp4
 # (see issue #57)
 CCP4_1JZV = importlib_resources.files(__name__) / '1jzv.ccp4'

gridData/tests/datafiles/ispg_0.mrc

@@ -113,6 +113,15 @@ def test_triclinic_ValueError():
                        "supported, not"):
         Grid(datafiles.MRC_EMD3001, file_format="MRC")
 
+def test_mrcfile_volume_check():
+    with pytest.raises(ValueError, match="is not a volumetric density"):
+        Grid(datafiles.ISPG_0)
+
+def test_mrcfile_volume_force():
+    grid = Grid(datafiles.ISPG_0, assume_volumetric=True)
+    assert_allclose(np.sum(grid.grid), 829.925)
+
+
 class TestGridMRC():
     @pytest.fixture(scope="class")
     def grid(self):