raw.set_montage(None) removes the MEG channels positions

Hi!

I am currently processing MEG empty room data for an experiment. For some of the participants, the digitization points were kept in the acquisition software while doing the empty room, which leads to issues later in the pipeline.

A solution we found was to use the method raw.set_montage(None), which as I understood would remove the digitization points from the data.

However, when running the following:

auto_noisy_channels_er, auto_flat_channels_er, auto_scores_er = mne.preprocessing.find_bad_channels_maxwell(
        raw_bids_empty_room,  # Raw data to process
        cross_talk = crosstalk_file,  # Specific to the MEG machine
        calibration = fine_cal_file,  # Specific to the MEG machine
        return_scores = True, 
        coord_frame = "meg",  # Necessary for empty room recordings
        duration = 30,  # Duration of each window (in seconds)
        min_count = 10,  # Number of window appearances to be counted as bad channel
        verbose = False,
        )

I get the following error:

ValueError                                Traceback (most recent call last)
/1_Preprocessing_batch.ipynb Cell 15 line 2
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=251'>252</a> print(\"Running the automatic detection of bad channels...\")
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=252'>253</a> t_step = dt.now()
--> <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=253'>254</a> auto_noisy_channels_er, auto_flat_channels_er, auto_scores_er = mne.preprocessing.find_bad_channels_maxwell(
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=254'>255</a>     raw_bids_empty_room,  # Raw data to process
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=255'>256</a>     cross_talk = subject.crosstalk_file, 
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=256'>257</a>     calibration = subject.fine_cal_file,
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=257'>258</a>     return_scores = True, 
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=258'>259</a>     coord_frame = \"meg\",  # Necessary for empty room recordings
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=259'>260</a>     duration = 30,  # Duration of each window (in seconds)
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=260'>261</a>     min_count = 10,  # Number of window appearances to be counted as bad channel
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=261'>262</a>     verbose = False,
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=262'>263</a>     )
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=263'>264</a> if show_timings:
    <a href='vscode-notebook-cell:/1_Preprocessing_batch.ipynb#X20sZmlsZQ%3D%3D?line=264'>265</a>     print(f\"\\tPerformed in {dt.now() - t_step}\")

File <decorator-gen-433>:10, in find_bad_channels_maxwell(raw, limit, duration, min_count, return_scores, origin, int_order, ext_order, calibration, cross_talk, coord_frame, regularize, ignore_ref, bad_condition, head_pos, mag_scale, skip_by_annotation, h_freq, extended_proj, verbose)

File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:2723, in find_bad_channels_maxwell(***failed resolving arguments***)
   2721 delta = chunk_raw.get_data(these_picks)
   2722 with use_log_level(False):
-> 2723     _run_maxwell_filter(chunk_raw, reconstruct=\"orig\", copy=False, **params)
   2725 if n_iter == 1 and len(chunk_flats):
   2726     logger.info(
   2727         \"            Flat (%2d): %s\",
   2728         len(chunk_flats),
   2729         \" \".join(chunk_flats),
   2730     )

File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:685, in _run_maxwell_filter(raw, skip_by_annotation, st_duration, st_correlation, st_only, st_when, ctc, coil_scale, this_pos_quat, meg_picks, good_mask, grad_picks, head_pos, info, _get_this_decomp_trans, S_recon, update_kwargs, ignore_ref, reconstruct, copy)
    657 def _run_maxwell_filter(
    658     raw,
    659     skip_by_annotation,
   (...)
    683     # The time it takes to recompute S and pS themselves is roughly on par
    684     # with the np.dot with the data, so not a huge gain to be made there.
--> 685     S_decomp, S_decomp_full, pS_decomp, reg_moments, n_use_in = _get_this_decomp_trans(
    686         info[\"dev_head_t\"], t=0.0
    687     )
    688     update_kwargs.update(reg_moments=reg_moments.copy())
    689     if ctc is not None:

File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:1240, in _get_decomp(trans, all_coils, cal, regularize, exp, ignore_ref, coil_scale, grad_picks, mag_picks, good_mask, mag_or_fine, bad_condition, t, mag_scale, mult)
   1235 del extended_proj
   1237 #
   1238 # Regularization
   1239 #
-> 1240 S_decomp, reg_moments, n_use_in = _regularize(
   1241     regularize, exp, S_decomp, mag_or_fine, extended_remove, t=t
   1242 )
   1243 S_decomp_full = S_decomp_full.take(reg_moments, axis=1)
   1245 #
   1246 # Pseudo-inverse of total multipolar moment basis set (Part of Eq. 37)
   1247 #

File <decorator-gen-431>:12, in _regularize(regularize, exp, S_decomp, mag_or_fine, extended_remove, t, verbose)

File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:1297, in _regularize(regularize, exp, S_decomp, mag_or_fine, extended_remove, t, verbose)
   1295 t_str = f\"{t:8.3f}\"
   1296 if regularize is not None:  # regularize='in'
-> 1297     in_removes, out_removes = _regularize_in(
   1298         int_order, ext_order, S_decomp, mag_or_fine, extended_remove
   1299     )
   1300 else:
   1301     in_removes = []

File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/preprocessing/maxwell.py:2294, in _regularize_in(int_order, ext_order, S_decomp, mag_or_fine, extended_remove)
   2292 for ii in range(n_in):
   2293     this_S = S_decomp.take(in_keepers + out_keepers, axis=1)
-> 2294     u, s, v = _safe_svd(this_S, full_matrices=False, **check_disable)
   2295     del this_S
   2296     eigs[ii] = s[[0, -1]]

File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/mne/fixes.py:92, in _safe_svd(A, **kwargs)
     90     raise ValueError(\"Cannot set overwrite_a=True with this function\")
     91 try:
---> 92     return linalg.svd(A, **kwargs)
     93 except np.linalg.LinAlgError as exp:
     94     from .utils import warn

File ~/.conda/envs/bodylingual313/lib/python3.13/site-packages/scipy/linalg/_decomp_svd.py:166, in svd(a, full_matrices, compute_uv, overwrite_a, check_finite, lapack_driver)
    164     raise LinAlgError(\"SVD did not converge\")
    165 if info < 0:
--> 166     raise ValueError('illegal value in %dth argument of internal gesdd'
    167                      % -info)
    168 if compute_uv:
    169     return u, s, v

ValueError: illegal value in 4th argument of internal gesdd

Now, I checked and it seems like using raw_bids_empty_room.set_montage(None) turned the gradiometers and magnetometers positions to nans. As an example, an excerpt of the output of raw_bids_empty_room.info["chs"] is:

[{'scanno': 1, 'logno': 113, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0113', 'coord_frame': 1 (FIFFV_COORD_DEVICE)}, 
{'scanno': 2, 'logno': 112, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0112', 'coord_frame': 1 (FIFFV_COORD_DEVICE)}, 
{'scanno': 3, 'logno': 111, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 4.1400001482561066e-11, 'coil_type': 3022 (FIFFV_COIL_VV_MAG_T1), 'loc': array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]), 'unit': 112 (FIFF_UNIT_T), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0111', 'coord_frame': 1 (FIFFV_COORD_DEVICE)}, 

I checked, and printing raw_bids_empty_room.info["chs"] before calling raw_bids_empty_room.set_montage(None) proves that the loc data was there to begin with:

[{'scanno': 1, 'logno': 113, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([-0.1066    ,  0.0464    , -0.0604    , -0.0127    ,  0.0057    ,
       -0.99990302, -0.186801  , -0.98240298, -0.0033    , -0.98232698,
        0.18674099,  0.013541  ]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0113', 'coord_frame': 1 (FIFFV_COORD_DEVICE)}, 
{'scanno': 2, 'logno': 112, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 3.250000046861601e-09, 'coil_type': 3012 (FIFFV_COIL_VV_PLANAR_T1), 'loc': array([-0.1066    ,  0.0464    , -0.0604    , -0.186801  , -0.98240298,
       -0.0033    ,  0.0127    , -0.0057    ,  0.99990302, -0.98232698,
        0.18674099,  0.013541  ]), 'unit': 201 (FIFF_UNIT_T_M), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0112', 'coord_frame': 1 (FIFFV_COORD_DEVICE)}, 
{'scanno': 3, 'logno': 111, 'kind': 1 (FIFFV_MEG_CH), 'range': 1.9073486328125e-05, 'cal': 4.1400001482561066e-11, 'coil_type': 3022 (FIFFV_COIL_VV_MAG_T1), 'loc': array([-0.1066    ,  0.0464    , -0.0604    , -0.0127    ,  0.0057    ,
       -0.99990302, -0.186801  , -0.98240298, -0.0033    , -0.98232698,
        0.18674099,  0.013541  ]), 'unit': 112 (FIFF_UNIT_T), 'unit_mul': 0 (FIFF_UNITM_NONE), 'ch_name': 'MEG0111', 'coord_frame': 1 (FIFFV_COORD_DEVICE)},

So, I tried to do something else to remove the digitization data:

montage = empty_room_fif.get_montage()  # We get the current montage
montage.remove_fiducials()  # We remove the fiducials
montage.dig = []  # We remove the digitization points

# At this point, the digitization is made of 0 points
empty_room_fif.set_montage(montage)  # We set the new montage

This worked, and allowed to remove the digitization points and coil data without removing the MEG channels positions.

I would like to know if this solution has a risk of breaking anything down the line - considering I am directly modifying an attribute instead of using a function, there might be some things I do not know about that I overlooked.

I compared the info with that of an empty room that didn’t have any digitization points to begin with. There are a few things that vary, notably dig is None in this one, instead of [] in the one we manually set (and I cannot set dig on None directly, as set_montage() then returns an error).

I am using Python 3.13.1 and MNE 1.9.0, working in a notebook on a Linux system.

Thank you!