MRG: Make use of AssociatedEmptyRoom field in *_meg.json (#795)

* Make use of AssociatedEmptyRoom field in *_meg.json

Fixes #493

* Remove new property

* flake

* Apply suggestions from code review

* Fix

* Missing whitespace

* Update changelog

* Try to fix Windows

* Another attempt to fix Windows

* Hopefully fix Windows tests

* Reformat

Author: hoechenberger

doc/whats_new.rst

@@ -41,6 +41,7 @@ Enhancements
 - :func:`mne_bids.write_raw_bids` gained a new keyword argument ``symlink``, which allows to create symbolic links to the original data files instead of copying them over. Currently works for ``FIFF`` files on macOS and Linux, by `Richard Höchenberger`_ (:gh:`778`)
 - :class:`mne_bids.BIDSPath` now has property getter and setter methods for all BIDS entities, i.e., you can now do things like ``bids_path.subject = 'foo'`` and don't have to resort to ``bids_path.update()``. This also ensures you'll get proper completion suggestions from your favorite Python IDE, by `Richard Höchenberger`_ (:gh:`786`)
 - :func:`mne_bids.write_raw_bids` now stores information about continuous head localization measurements (e.g., Elekta/Neuromag cHPI) in the MEG sidecar file, by `Richard Höchenberger`_ (:gh:`794`)
+- :func:`mne_bids.write_raw_bids` gained a new parameter `empty_room` that allows to specify an associated empty-room recording when writing an MEG data file. This information will be stored in the ``AssociatedEmptyRoom`` field of the MEG JSON sidecar file, by `Richard Höchenberger`_ (:gh:`795`)
 
 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^
@@ -49,6 +50,7 @@ API and behavior changes
 - When writing BIDS datasets, MNE-BIDS now tags them as BIDS 1.6.0 (we previously tagged them as BIDS 1.4.0), by `Richard Höchenberger`_ (:gh:`782`)
 - :func:`mne_bids.read_raw_bids` now passes ``allow_maxshield=True`` to the MNE-Python reader function by default when reading FIFF files. Previously, ``extra_params=dict(allow_maxshield=True)`` had to be passed explicitly, by `Richard Höchenberger`_ (:gh:`#787`)
 - The ``raw_to_bids`` command has lost its ``--allow_maxshield`` parameter. If writing a FIFF file, we will now always assume that writing data before applying a Maxwell filter is fine, by `Richard Höchenberger`_ (:gh:`#787`)
+- :meth:`mne_bids.BIDSPath.find_empty_room` now first looks for an ``AssociatedEmptyRoom`` field in the MEG JSON sidecar file to retrieve the empty-room recording; only if this information is missing, it will proceed to try and find the best-matching empty-room recording based on measurement date (i.e., fall back to the previous behavior), by `Richard Höchenberger`_ (:gh:`#795`)
 
 Requirements
 ^^^^^^^^^^^^

mne_bids/path.py

@@ -12,6 +12,7 @@
 from os import path as op
 from pathlib import Path
 from datetime import datetime
+import json
 from typing import Optional, Union
 
 import numpy as np
@@ -27,7 +28,7 @@
                             param_regex, _ensure_tuple)
 
 
-def _get_matched_empty_room(bids_path):
+def _find_matched_empty_room(bids_path):
     """Get matching empty-room file for an MEG recording."""
     # Check whether we have a BIDS root.
     bids_root = bids_path.root
@@ -43,11 +44,7 @@ def _get_matched_empty_room(bids_path):
     bids_fname = bids_path.update(suffix=datatype,
                                   root=bids_root).fpath
     _, ext = _parse_ext(bids_fname)
-    extra_params = None
-    if ext == '.fif':
-        extra_params = dict(allow_maxshield=True)
-
-    raw = read_raw_bids(bids_path=bids_path, extra_params=extra_params)
+    raw = read_raw_bids(bids_path=bids_path)
     if raw.info['meas_date'] is None:
         raise ValueError('The provided recording does not have a measurement '
                          'date set. Cannot get matching empty-room file.')
@@ -863,14 +860,49 @@ def find_empty_room(self):
         This will only work if the ``.root`` attribute of the
         :class:`mne_bids.BIDSPath` instance has been set.
 
+        .. note:: If the sidecar JSON file contains an ``AssociatedEmptyRoom``
+                  entry, the empty-room recording specified there will be used.
+                  Otherwise, this method will try to find the best-matching
+                  empty-room recording based on measurement date.
+
         Returns
         -------
         BIDSPath | None
             The path corresponding to the best-matching empty-room measurement.
-            Returns None if none was found.
-
+            Returns ``None`` if none was found.
         """
-        return _get_matched_empty_room(self)
+        if self.datatype not in ('meg', None):
+            raise ValueError('Empty-room data is only supported for MEG '
+                             'datasets')
+
+        if self.root is None:
+            raise ValueError('The root of the "bids_path" must be set. '
+                             'Please use `bids_path.update(root="<root>")` '
+                             'to set the root of the BIDS folder to read.')
+
+        sidecar_fname = _find_matching_sidecar(self, extension='.json')
+        with open(sidecar_fname, 'r', encoding='utf-8') as f:
+            sidecar_json = json.load(f)
+
+        if 'AssociatedEmptyRoom' in sidecar_json:
+            logger.info('Using "AssociatedEmptyRoom" entry from MEG sidecar '
+                        'file to retrieve empty-room path.')
+            emptytoom_path = sidecar_json['AssociatedEmptyRoom']
+            emptyroom_entities = get_entities_from_fname(emptytoom_path)
+            er_bids_path = BIDSPath(root=self.root, datatype='meg',
+                                    **emptyroom_entities)
+        else:
+            logger.info(
+                'The MEG sidecar file does not contain an '
+                '"AssociatedEmptyRoom" entry. Will try to find a matching '
+                'empty-room recording based on the measurement date …'
+            )
+            er_bids_path = _find_matched_empty_room(self)
+
+        if er_bids_path is not None:
+            assert er_bids_path.fpath.exists()
+
+        return er_bids_path
 
     @property
     def meg_calibration_fpath(self):

mne_bids/tests/test_path.py

@@ -821,7 +821,7 @@ def test_find_empty_room(return_bids_test_dir, tmpdir):
                        match='The root of the "bids_path" must be set'):
         bids_path.copy().update(root=None).find_empty_room()
 
-    # assert that we get error if meas_date is not available.
+    # assert that we get an error if meas_date is not available.
     raw = read_raw_bids(bids_path=bids_path)
     raw.set_meas_date(None)
     anonymize_info(raw.info)
@@ -830,6 +830,51 @@ def test_find_empty_room(return_bids_test_dir, tmpdir):
                                          'have a measurement date set'):
         bids_path.find_empty_room()
 
+    # test that the `AssociatedEmptyRoom` key in MEG sidecar is respected
+
+    bids_root = tmpdir.mkdir('associated-empty-room')
+    raw = _read_raw_fif(raw_fname)
+    meas_date = datetime(year=2020, month=1, day=10, tzinfo=timezone.utc)
+    er_date = datetime(year=2010, month=1, day=1, tzinfo=timezone.utc)
+    raw.set_meas_date(meas_date)
+
+    er_raw_matching_date = er_raw.copy().set_meas_date(meas_date)
+    er_raw_associated = er_raw.copy().set_meas_date(er_date)
+
+    # First write empty-room data
+    # We write two empty-room recordings: one with a date matching exactly the
+    # experimental measurement date, and one dated approx. 10 years earlier
+    # We will want to enforce using the older recording via
+    # `AssociatedEmptyRoom` (without AssociatedEmptyRoom, find_empty_room()
+    # would return the recording with the matching date instead)
+    er_matching_date_bids_path = BIDSPath(
+        subject='emptyroom', session='20200110', task='noise', root=bids_root,
+        datatype='meg', suffix='meg', extension='.fif')
+    write_raw_bids(er_raw_matching_date, bids_path=er_matching_date_bids_path)
+
+    er_associated_bids_path = (er_matching_date_bids_path.copy()
+                               .update(session='20100101'))
+    write_raw_bids(er_raw_associated, bids_path=er_associated_bids_path)
+
+    # Now we write experimental data and associate it with the earlier
+    # empty-room recording
+    bids_path = (er_matching_date_bids_path.copy()
+                 .update(subject='01', session=None, task='task'))
+    write_raw_bids(raw, bids_path=bids_path,
+                   empty_room=er_associated_bids_path)
+
+    # Retrieve empty-room BIDSPath
+    assert bids_path.find_empty_room() == er_associated_bids_path
+
+    # Should only work for MEG
+    with pytest.raises(ValueError, match='only supported for MEG'):
+        bids_path.copy().update(datatype='eeg').find_empty_room()
+
+    # Don't create `AssociatedEmptyRoom` entry in sidecar – we should now
+    # retrieve the empty-room recording closer in time
+    write_raw_bids(raw, bids_path=bids_path, empty_room=None, overwrite=True)
+    assert bids_path.find_empty_room() == er_matching_date_bids_path
+
 
 @pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
 def test_find_emptyroom_ties(tmpdir):

mne_bids/tests/test_write.py

@@ -493,7 +493,7 @@ def test_fif(_bids_validate, tmpdir):
     # test that an incorrect date raises an error.
     er_bids_basename_bad = BIDSPath(subject='emptyroom', session='19000101',
                                     task='noise', root=bids_root)
-    with pytest.raises(ValueError, match='Date provided'):
+    with pytest.raises(ValueError, match='The date provided'):
         write_raw_bids(raw, er_bids_basename_bad, overwrite=False)
 
     # test that the acquisition time was written properly
@@ -2736,3 +2736,38 @@ def test_symlink(tmpdir):
     p = write_raw_bids(raw=raw, bids_path=bids_path, symlink=True)
     raw = read_raw_bids(p)
     assert len(raw.filenames) == 2
+
+
+@pytest.mark.filterwarnings(warning_str['channel_unit_changed'])
+def test_write_associated_emptyroom(_bids_validate, tmpdir):
+    """Test functionality of the write_raw_bids conversion for fif."""
+    bids_root = tmpdir.mkdir('bids1')
+    data_path = testing.data_path()
+    raw_fname = op.join(data_path, 'MEG', 'sample',
+                        'sample_audvis_trunc_raw.fif')
+    raw = _read_raw_fif(raw_fname)
+    meas_date = datetime(year=2020, month=1, day=10, tzinfo=timezone.utc)
+    raw.set_meas_date(meas_date)
+
+    # First write "empty-room" data
+    bids_path_er = BIDSPath(subject='emptyroom', session='20200110',
+                            task='noise', root=bids_root, datatype='meg',
+                            suffix='meg', extension='.fif')
+    write_raw_bids(raw, bids_path=bids_path_er)
+
+    # Now we write experimental data and associate it with the empty-room
+    # recording
+    bids_path = bids_path_er.copy().update(subject='01', session=None,
+                                           task='task')
+    write_raw_bids(raw, bids_path=bids_path, empty_room=bids_path_er)
+    _bids_validate(bids_path.root)
+
+    meg_json_path = bids_path.copy().update(extension='.json')
+    with open(meg_json_path, 'r') as fin:
+        meg_json_data = json.load(fin)
+
+    assert 'AssociatedEmptyRoom' in meg_json_data
+    assert (bids_path_er.fpath
+            .as_posix()  # make test work on Windows, too
+            .endswith(meg_json_data['AssociatedEmptyRoom']))
+    assert meg_json_data['AssociatedEmptyRoom'].startswith('/')

mne_bids/write.py

@@ -564,8 +564,8 @@ def _mri_scanner_ras_to_mri_voxels(ras_landmarks, img_mgh):
     return vox_landmarks
 
 
-def _sidecar_json(raw, task, manufacturer, fname, datatype, overwrite=False,
-                  verbose=True):
+def _sidecar_json(raw, task, manufacturer, fname, datatype,
+                  emptyroom_fname=None, overwrite=False, verbose=True):
     """Create a sidecar json file depending on the suffix and save it.
 
     The sidecar json file provides meta data about the data
@@ -584,6 +584,9 @@ def _sidecar_json(raw, task, manufacturer, fname, datatype, overwrite=False,
         Filename to save the sidecar json to.
     datatype : str
         Type of the data as in ALLOWED_ELECTROPHYSIO_DATATYPE.
+    emptyroom_fname : str | mne_bids.BIDSPath
+        For MEG recordings, the path to an empty-room data file to be
+        associated with ``raw``. Only supported for MEG.
     overwrite : bool
         Whether to overwrite the existing file.
         Defaults to False.
@@ -675,6 +678,7 @@ def _sidecar_json(raw, task, manufacturer, fname, datatype, overwrite=False,
         ('SoftwareFilters', 'n/a'),
         ('RecordingDuration', raw.times[-1]),
         ('RecordingType', rec_type)]
+
     ch_info_json_meg = [
         ('DewarPosition', 'n/a'),
         ('DigitizedLandmarks', digitized_landmark),
@@ -683,15 +687,21 @@ def _sidecar_json(raw, task, manufacturer, fname, datatype, overwrite=False,
         ('MEGREFChannelCount', n_megrefchan),
         ('ContinuousHeadLocalization', chpi),
         ('HeadCoilFrequency', list(hpi_freqs))]
+
+    if emptyroom_fname is not None:
+        ch_info_json_meg.append(('AssociatedEmptyRoom', str(emptyroom_fname)))
+
     ch_info_json_eeg = [
         ('EEGReference', 'n/a'),
         ('EEGGround', 'n/a'),
         ('EEGPlacementScheme', _infer_eeg_placement_scheme(raw)),
         ('Manufacturer', manufacturer)]
+
     ch_info_json_ieeg = [
         ('iEEGReference', 'n/a'),
         ('ECOGChannelCount', n_ecogchan),
         ('SEEGChannelCount', n_seegchan)]
+
     ch_info_ch_counts = [
         ('EEGChannelCount', n_eegchan),
         ('EOGChannelCount', n_eogchan),
@@ -960,6 +970,7 @@ def make_dataset_description(path, name, data_license=None,
 def write_raw_bids(raw, bids_path, events_data=None,
                    event_id=None, anonymize=None,
                    format='auto', symlink=False,
+                   empty_room=None,
                    overwrite=False, verbose=True):
     """Save raw data to a BIDS-compliant folder structure.
 
@@ -1083,6 +1094,12 @@ def write_raw_bids(raw, bids_path, events_data=None,
            Symlinks are currently only supported on macOS and Linux. We will
            add support for Windows 10 at a later time.
 
+    empty_room : mne_bids.BIDSPath | None
+        The empty-room recording to be associated with this file. This is
+        only supported for MEG data, and only if the ``root`` attributes of
+        ``bids_path`` and ``empty_room`` are the same. Pass ``None``
+        (default) if you do not wish to specify an associated empty-room
+        recording.
     overwrite : bool
         Whether to overwrite existing files or data in files.
         Defaults to ``False``.
@@ -1189,6 +1206,9 @@ def write_raw_bids(raw, bids_path, events_data=None,
         raise RuntimeError('You passed event_id, but no events_data NumPy '
                            'array. You need to pass both, or neither.')
 
+    _validate_type(item=empty_room, item_name='empty_room',
+                   types=(BIDSPath, None))
+
     raw = raw.copy()
 
     raw_fname = raw.filenames[0]
@@ -1232,10 +1252,10 @@ def write_raw_bids(raw, bids_path, events_data=None,
 
     # check whether the info provided indicates that the data is emptyroom
     # data
-    emptyroom = False
+    data_is_emptyroom = False
     if (bids_path.datatype == 'meg' and bids_path.subject == 'emptyroom' and
             bids_path.task == 'noise'):
-        emptyroom = True
+        data_is_emptyroom = True
         # check the session date provided is consistent with the value in raw
         meas_date = raw.info.get('meas_date', None)
         if meas_date is not None:
@@ -1244,13 +1264,40 @@ def write_raw_bids(raw, bids_path, events_data=None,
                                                    tz=timezone.utc)
             er_date = meas_date.strftime('%Y%m%d')
             if er_date != bids_path.session:
-                raise ValueError("Date provided for session doesn't match "
-                                 "session date.")
+                raise ValueError(
+                    f"The date provided for the empty-room session "
+                    f"({bids_path.session}) doesn't match the empty-room "
+                    f"recording date found in the data's info structure "
+                    f"({er_date})."
+                )
             if anonymize is not None and 'daysback' in anonymize:
                 meas_date = meas_date - timedelta(anonymize['daysback'])
                 session = meas_date.strftime('%Y%m%d')
                 bids_path = bids_path.copy().update(session=session)
 
+    associated_er_path = None
+    if empty_room is not None:
+        if bids_path.datatype != 'meg':
+            raise ValueError('"empty_room" is only supported for '
+                             'MEG data.')
+        if data_is_emptyroom:
+            raise ValueError('You cannot write empty-room data and pass '
+                             '"empty_room" at the same time.')
+        if bids_path.root != empty_room.root:
+            raise ValueError('The MEG data and its associated empty-room '
+                             'recording must share the same BIDS root.')
+
+        associated_er_path = empty_room.fpath
+        if not associated_er_path.exists():
+            raise FileNotFoundError(f'Empty-room data file not found: '
+                                    f'{associated_er_path}')
+
+        # Turn it into a path relative to the BIDS root
+        associated_er_path = Path(str(associated_er_path)
+                                  .replace(str(empty_room.root), ''))
+        # Ensure it works on Windows too
+        associated_er_path = associated_er_path.as_posix()
+
     data_path = bids_path.mkdir().directory
 
     # In case of an "emptyroom" subject, BIDSPath() will raise
@@ -1317,7 +1364,7 @@ def write_raw_bids(raw, bids_path, events_data=None,
     _participants_json(participants_json_fname, True, verbose)
 
     # for MEG, we only write coordinate system
-    if bids_path.datatype == 'meg' and not emptyroom:
+    if bids_path.datatype == 'meg' and not data_is_emptyroom:
         _write_coordsystem_json(raw=raw, unit=unit, hpi_coord_system=orient,
                                 sensor_coord_system=orient,
                                 fname=coordsystem_path.fpath,
@@ -1333,7 +1380,7 @@ def write_raw_bids(raw, bids_path, events_data=None,
                        f'for data type "{bids_path.datatype}". Skipping ...')
 
     # Write events.
-    if not emptyroom:
+    if not data_is_emptyroom:
         events_array, event_dur, event_desc_id_map = _read_events(
             events_data, event_id, raw, verbose=False
         )
@@ -1351,8 +1398,10 @@ def write_raw_bids(raw, bids_path, events_data=None,
     make_dataset_description(bids_path.root, name=" ", overwrite=False,
                              verbose=verbose)
 
-    _sidecar_json(raw, bids_path.task, manufacturer, sidecar_path.fpath,
-                  bids_path.datatype, overwrite, verbose)
+    _sidecar_json(raw, task=bids_path.task, manufacturer=manufacturer,
+                  fname=sidecar_path.fpath, datatype=bids_path.datatype,
+                  emptyroom_fname=associated_er_path,
+                  overwrite=overwrite, verbose=verbose)
     _channels_tsv(raw, channels_path.fpath, overwrite, verbose)
 
     # create parent directories if needed