MRG: Make get_head_mri_trans() respect datatype and suffix attributes…

… of MEG BIDSPath (#969)

* Make get_head_mri_trans() respect datatype and suffix attributes

get_head_mri_trans() currently blindly replaces the datatype and suffix
of the provided MEG BIDSPath with 'meg'. However, when working with
derivative data, the suffix might be different, causing the function to
fail as it will create a BIDSPath that points to a non-existent file.

The change here makes it such that if the provided MEG BIDSPath has
a suffix and / or datatype attribute set, they will not be modified.
Only if the attributes are not set, we will fall back to 'meg'.

* Add a test to convince Alex

* Update changelog

* flake

* Implement changes discussed with Alex

* Update docstring

* Ensure we're working with head coords

* Update another test

* Simplify and cleanups

* Fix

* Style

doc/whats_new.rst

@@ -43,13 +43,15 @@ Enhancements
 
 - Add the ability to write multiple landmarks with :func:`mne_bids.write_anat` (e.g. to have separate landmarks for different sessions) via the new ``kind`` parameter, by `Alexandre Gramfort`_ (:gh:`955`)
 
-- :func:`mne_bids.get_head_mri_trans` and :func:`mne_bids.update_anat_landmarks` gained a new``kind`` parameter to specify which of multiple landmark sets to retrieve, by `Alexandre Gramfort`_ and `Richard Höchenberger`_ (:gh:`955`, :gh:`957`)
+- Similarly, :func:`mne_bids.get_head_mri_trans` and :func:`mne_bids.update_anat_landmarks` gained a new ``kind`` parameter to specify which of multiple landmark sets to operate on, by `Alexandre Gramfort`_ and `Richard Höchenberger`_ (:gh:`955`, :gh:`957`)
 
 API and behavior changes
 ^^^^^^^^^^^^^^^^^^^^^^^^
 
 - :func:`mne_bids.update_anat_landmarks` will now by default raise an exception if the requested MRI landmarks do not already exist. Use the new ``on_missing`` parameter to control this behavior, by `Richard Höchenberger`_ (:gh:`957`)
 
+- :func:`mne_bids.get_head_mri_trans` now raises a warning if ``datatype`` or ``suffix`` of the provided electrophysiological :class:`mne_bids.BIDSPath` are not set. In the future, this will raise an exception, by `Richard Höchenberger`_(:gh:`969`)
+
 Requirements
 ^^^^^^^^^^^^
 
@@ -70,6 +72,8 @@ Bug fixes
 
 - Avoid modifying the instance of :class:`mne_bids.BIDSPath` if validation fails when calling :meth:`mne_bids.BIDSPath.update`, by `Alexandre Gramfort`_ (:gh:`950`)
 
+- :func:`mne_bids.get_head_mri_trans` now respects ``datatype`` and ``suffix`` of the provided electrophysiological :class:`mne_bids.BIDSPath`, simplifying e.g. reading of derivaties, by `Richard Höchenberger`_ (:gh:`969`)
+
 :doc:`Find out what was new in previous releases <whats_new_previous_releases>`
 
 .. include:: authors.rst

mne_bids/dig.py

@@ -73,7 +73,6 @@ def _handle_coordsystem_reading(coordsystem_fpath, datatype):
     Handle reading the coordinate frame and coordinate unit
     of each electrode.
     """
-    # open coordinate system sidecar json
     with open(coordsystem_fpath, 'r', encoding='utf-8-sig') as fin:
         coordsystem_json = json.load(fin)
 
@@ -93,8 +92,9 @@ def _handle_coordsystem_reading(coordsystem_fpath, datatype):
         coord_frame_desc = coordsystem_json.get('iEEGCoordinateDescription',
                                                 None)
 
-    logger.info(f"Reading in coordinate system frame {coord_frame}: "
-                f"{coord_frame_desc}.")
+    msg = f'Reading coordinate system frame {coord_frame}'
+    if coord_frame_desc:
+        msg += f': {coord_frame_desc}'
 
     return coord_frame, coord_unit
 
@@ -450,15 +450,14 @@ def _read_dig_bids(electrodes_fpath, coordsystem_fpath,
     montage : mne.channels.DigMontage
         The coordinate data as MNE-Python DigMontage object.
     """
-    # read in coordinate information
     bids_coord_frame, bids_coord_unit = _handle_coordsystem_reading(
         coordsystem_fpath, datatype)
 
     if datatype == 'meg':
         if bids_coord_frame not in BIDS_MEG_COORDINATE_FRAMES:
-            warn("MEG Coordinate frame is not accepted "
-                 "BIDS keyword. The allowed keywords are: "
-                 "{}".format(BIDS_MEG_COORDINATE_FRAMES))
+            warn(f'MEG coordinate frame "{bids_coord_frame}" is not '
+                 f'supported. The supported coordinate frames are: '
+                 f'{BIDS_MEG_COORDINATE_FRAMES}')
             coord_frame = None
         elif bids_coord_frame == 'Other':
             warn("Coordinate frame of MEG data can't be determined "

mne_bids/read.py

@@ -777,7 +777,12 @@ def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None,
     Parameters
     ----------
     bids_path : BIDSPath
-        The path of the electrophysiology recording.
+        The path of the electrophysiology recording. If ``datatype`` and
+        ``suffix`` are not present, they will be set to ``'meg'``, and a
+        warning will be raised.
+
+        .. versionchanged:: 0.10
+           A warning is raised it ``datatype`` or ``suffix`` are not set.
     extra_params : None | dict
         Extra parameters to be passed to :func:`mne.io.read_raw` when reading
         the MEG file.
@@ -832,8 +837,17 @@ def get_head_mri_trans(bids_path, extra_params=None, t1_bids_path=None,
                          'Please use `bids_path.update(root="<root>")` '
                          'to set the root of the BIDS folder to read.')
 
-    # only get this for MEG data
-    meg_bids_path.update(datatype='meg', suffix='meg')
+    # if the bids_path is underspecified, only get info for MEG data
+    if meg_bids_path.datatype is None:
+        meg_bids_path.datatype = 'meg'
+        warn('bids_path did not have a datatype set. Assuming "meg". This '
+             'will raise an exception in the future.', module='mne_bids',
+             category=DeprecationWarning)
+    if meg_bids_path.suffix is None:
+        meg_bids_path.suffix = 'meg'
+        warn('bids_path did not have a suffix set. Assuming "meg". This '
+             'will raise an exception in the future.', module='mne_bids',
+             category=DeprecationWarning)
 
     # Get the sidecar file for MRI landmarks
     t1w_bids_path = (

mne_bids/tests/test_read.py

@@ -209,7 +209,9 @@ def test_get_head_mri_trans(tmp_path):
 
     # Write it to BIDS
     raw = _read_raw_fif(raw_fname)
-    bids_path = _bids_path.copy().update(root=tmp_path)
+    bids_path = _bids_path.copy().update(
+        root=tmp_path, datatype='meg', suffix='meg'
+    )
     write_raw_bids(raw, bids_path, events_data=events, event_id=event_id,
                    overwrite=False)
 
@@ -230,9 +232,13 @@ def test_get_head_mri_trans(tmp_path):
     landmarks = get_anat_landmarks(
         t1w_mgh, raw.info, trans, fs_subject='sample',
         fs_subjects_dir=subjects_dir)
+    t1w_bids_path = bids_path.copy().update(
+        datatype='anat', suffix='T1w'
+    )
     t1w_bids_path = write_anat(
-        t1w_mgh, bids_path=bids_path, landmarks=landmarks, verbose=True)
-    anat_dir = bids_path.directory
+        t1w_mgh, bids_path=t1w_bids_path, landmarks=landmarks, verbose=True
+    )
+    anat_dir = t1w_bids_path.directory
 
     # Try to get trans back through fitting points
     estimated_trans = get_head_mri_trans(
@@ -258,8 +264,9 @@ def test_get_head_mri_trans(tmp_path):
     raw = _read_raw_fif(raw_fname)
     write_raw_bids(raw, bids_path, events_data=events, event_id=event_id,
                    overwrite=True)  # overwrite with new acq
-    t1w_bids_path = write_anat(t1w_mgh, bids_path=bids_path,
-                               landmarks=landmarks, overwrite=True)
+    t1w_bids_path = write_anat(
+        t1w_mgh, bids_path=t1w_bids_path, landmarks=landmarks, overwrite=True
+    )
 
     t1w_json_fpath = t1w_bids_path.copy().update(extension='.json').fpath
     with t1w_json_fpath.open('r', encoding='utf-8') as f:
@@ -281,7 +288,9 @@ def test_get_head_mri_trans(tmp_path):
     # Test t1_bids_path parameter
     #
     # Case 1: different BIDS roots
-    meg_bids_path = _bids_path.copy().update(root=tmp_path / 'meg_root')
+    meg_bids_path = _bids_path.copy().update(
+        root=tmp_path / 'meg_root', datatype='meg', suffix='meg'
+    )
     t1_bids_path = _bids_path.copy().update(
         root=tmp_path / 'mri_root', task=None, run=None
     )
@@ -299,10 +308,12 @@ def test_get_head_mri_trans(tmp_path):
 
     # Case 2: different sessions
     raw = _read_raw_fif(raw_fname)
-    meg_bids_path = _bids_path.copy().update(root=tmp_path / 'session_test',
-                                             session='01')
+    meg_bids_path = _bids_path.copy().update(
+        root=tmp_path / 'session_test', session='01', datatype='meg',
+        suffix='meg'
+    )
     t1_bids_path = meg_bids_path.copy().update(
-        session='02', task=None, run=None
+        session='02', task=None, run=None, datatype='anat', suffix='T1w'
     )
 
     write_raw_bids(raw, bids_path=meg_bids_path)
@@ -337,6 +348,71 @@ def test_get_head_mri_trans(tmp_path):
         kind="deface")
     assert not np.allclose(trans['trans'], read_trans2['trans'])
 
+    # Test we're respecting existing suffix & data type
+    # The following path is supposed to mimic a derivative generated by the
+    # MNE-BIDS-Pipeline.
+    #
+    # XXX We MAY want to revise this once the BIDS-Pipeline produces more
+    # BIDS-compatible output, e.g. including `channels.tsv` files for written
+    # Raw data etc.
+    raw = _read_raw_fif(raw_fname)
+    deriv_root = tmp_path / 'derivatives' / 'mne-bids-pipeline'
+    electrophys_path = (
+        deriv_root / 'sub-01' / 'eeg' / 'sub-01_task-av_proc-filt_raw.fif'
+    )
+    electrophys_path.parent.mkdir(parents=True)
+    raw.save(electrophys_path)
+
+    electrophys_bids_path = BIDSPath(
+        subject='01', task='av', datatype='eeg', processing='filt',
+        suffix='raw', extension='.fif', root=deriv_root,
+        check=False
+    )
+    t1_bids_path = _bids_path.copy().update(
+        root=tmp_path / 'mri_root', task=None, run=None
+    )
+    with pytest.warns(RuntimeWarning, match='Did not find any channels.tsv'):
+        get_head_mri_trans(
+            bids_path=electrophys_bids_path,
+            t1_bids_path=t1_bids_path,
+            fs_subject='sample',
+            fs_subjects_dir=subjects_dir
+        )
+
+    # bids_path without datatype is deprecated
+    bids_path = electrophys_bids_path.copy().update(datatype=None)
+    with pytest.raises(FileNotFoundError):  # defaut location is all wrong!
+        with pytest.warns(DeprecationWarning, match='no datatype'):
+            get_head_mri_trans(
+                bids_path=bids_path,
+                t1_bids_path=t1_bids_path,
+                fs_subject='sample',
+                fs_subjects_dir=subjects_dir
+            )
+
+    # bids_path without suffix is deprecated
+    bids_path = electrophys_bids_path.copy().update(suffix=None)
+    with pytest.raises(FileNotFoundError):  # defaut location is all wrong!
+        with pytest.warns(DeprecationWarning, match='no datatype'):
+            get_head_mri_trans(
+                bids_path=bids_path,
+                t1_bids_path=t1_bids_path,
+                fs_subject='sample',
+                fs_subjects_dir=subjects_dir
+            )
+
+    # Should fail for an unsupported coordinate frame
+    raw = _read_raw_fif(raw_fname)
+    bids_root = tmp_path / 'unsupported_coord_frame'
+    bids_path = BIDSPath(
+        subject='01', task='av', datatype='meg', suffix='meg',
+        extension='.fif', root=bids_root
+    )
+    t1_bids_path = _bids_path.copy().update(
+        root=tmp_path / 'mri_root', task=None, run=None
+    )
+    write_raw_bids(raw=raw, bids_path=bids_path, verbose=False)
+
 
 def test_handle_events_reading(tmp_path):
     """Test reading events from a BIDS events.tsv file."""
@@ -849,7 +925,9 @@ def test_get_head_mri_trans_ctf(fname, tmp_path):
     ctf_data_path = op.join(testing.data_path(), 'CTF')
     raw_ctf_fname = op.join(ctf_data_path, fname)
     raw_ctf = _read_raw_ctf(raw_ctf_fname, clean_names=True)
-    bids_path = _bids_path.copy().update(root=tmp_path)
+    bids_path = _bids_path.copy().update(
+        root=tmp_path, datatype='meg', suffix='meg'
+    )
     write_raw_bids(raw_ctf, bids_path, overwrite=False)
 
     # Take a fake trans

mne_bids/utils.py

@@ -275,15 +275,26 @@ def _infer_eeg_placement_scheme(raw):
 def _extract_landmarks(dig):
     """Extract NAS, LPA, and RPA from raw.info['dig']."""
     coords = dict()
+    coord_frame = dict()
+
     landmarks = {d['ident']: d for d in dig
                  if d['kind'] == FIFF.FIFFV_POINT_CARDINAL}
     if landmarks:
         if FIFF.FIFFV_POINT_NASION in landmarks:
             coords['NAS'] = landmarks[FIFF.FIFFV_POINT_NASION]['r'].tolist()
+            coord_frame['NAS'] = (landmarks[FIFF.FIFFV_POINT_NASION]
+                                  ['coord_frame'])
         if FIFF.FIFFV_POINT_LPA in landmarks:
             coords['LPA'] = landmarks[FIFF.FIFFV_POINT_LPA]['r'].tolist()
+            coord_frame['LPA'] = landmarks[FIFF.FIFFV_POINT_LPA]['coord_frame']
         if FIFF.FIFFV_POINT_RPA in landmarks:
             coords['RPA'] = landmarks[FIFF.FIFFV_POINT_RPA]['r'].tolist()
+            coord_frame['RPA'] = landmarks[FIFF.FIFFV_POINT_RPA]['coord_frame']
+
+    # for now, we only support "head" coordinates
+    for frame in coord_frame.values():
+        assert frame == FIFF.FIFFV_COORD_HEAD
+
     return coords