I am doing some microstate analysis and would like to compare the actual topological maps at select timepoints to the microstate topological maps at those same timepoints. I can see on evoked data it is as simple as evoked.plot_topomap([time]). But I need to do this on resting state data, so I cannot turn it into an evoked instance. Is there a way to get the topomaps for select time points on raw or epoched data? I have been playing around with mne.viz.plot_topomap(epoched, layout, ch_type=‘mag’), but this does not seem like it is going to work. I have to flatten the epochs array first and then slice out an individual time point for it to work. Any ideas?
You are on the right track, the function mne.viz.plot_topomap() accepts arbitrary data values and channel locations. Therefore, you need to extract values of interest (i.e. for a selected time point) for each channel from your raw data and pass it as the first argument. The pos argument can be an Info object with an assigned montage.
Hi, thanks for your response. I am having trouble with the pos argument due to the data I am using. I am using the HCP MEG data, which has an issue with the way the sensors are plotted in MNE, they are 45 degrees off from the head model.
I have tried figuring out how to rotate the head or the sensors but have had no luck so far. The data from the HCP has no montage in the info object. I have tried putting one in, but all I can find is the magnesWH3600 dig and layout files. The dig file ends up messing up the locations of the sensors further due to the incongruency above. the layout file does seem to work fine when I use it as an augment in some of the epochs plotting functions, but for this particular function it plots off of the head model
This is the code I am using, at some point when I get it working I will loop it so I can iteratively slice different time points from the flattened data array. Right not flat[:,0] takes just the first time point to test the code.
epoched=mne.read_epochs(path)
#%% get correct layout
layout = mne.channels.read_layout(‘magnesWH3600’)
layout.names = [meg.replace('MEG ', ‘’) for meg in layout.names]
layout.names = [z.lstrip(‘0’) for z in layout.names]
layout.names = [‘A’ + a for a in layout.names]
#%% flatten epoched data
count=0
for count in range(len(epoched.get_data())):
if count==0:
flat=epoched.get_data()[count]
else:
flat=np.concatenate((flat, epoched.get_data()[count]), axis=1)
count+=1
#%% delete ch pos for any sensor participant is missing
bads=[ ]
for x in range(len(layout.names)):
if layout.names[ x ] not in epoched.ch_names:
bads.append(x)
else: continue
bads.reverse()
for y in bads:
layout.pos=np.delete(layout.pos,y,0)
#%% plot timepoint topomap for first time point
mne.viz.plot_topomap(flat[:,0], layout.pos, ch_type=‘mag’)
Regarding your first plot, where sensors are rotated by 45°, are these channels locations already available in your MEG data? Or did you generate/import them?
The locations are already in the data. It was originally a matlab format that I changed over to a fif. There was very little in the Matlab format that goes into the fif info object. Seems like the channel locations are brought over and placed in the epoched.info[‘chs’]. The other fields such as info[‘dig’] are empty. Putting in the MNE standard dig file for the Magness sensors ends up doing nothing on the plots, but messes up things like sensor interpolation. I believe this is a known issue with the HCP MEG data. And sorry, I meant 90 degrees off, as you can see in that first plot the front of the sensors is oriented towards the right ear on the head model.
I guess trying to rotate the locations would be the best approach since this is basically almost what you want. I’m not sure if there is a way to do that with MNE functions though, but maybe someone with more experience could chime in? @mmagnuski maybe?
If the data was in eeglab format then IIRC x and y dimentions are different to mne (or was it fieldtrip?). I would first just switch x and y dimensions in the info and see if you like the outcome.
Maybe there is an easier way to do it, but I would just modify loc field of each entry in info[‘chs’].
Something like:
for ch_idx in range(len(info['chs'])):
current_pos = info['chs'][ch_idx]['loc'][:3]
info['chs'][ch_idx]['loc'][[1, 0, 2] = current_pos
Thanks @cbrnr and @mmagnuski, that loc swap did fix the first plot of the sensors being 90 degrees rotated. However, it di not help the second plot produced by mne.viz.plot_topomap(flat[:,0], layout.pos, ch_type=‘mag’). The topomap of the sensor fields is still off of the head model. I could delete the head model but since the sensor fields are off center it is messing up the plot when I try to loop this function with matplotlib so I can have a single plot with a topomap for every few milliseconds.
Also, yes you are right, it is fieldtrip the data was in originally. My schools server has the HCP data, but I was never able to get fieldtrip to work on their server so I defaulted to MNE.
I am aware of the hcp-mne package and tried using it originally, but it is not maintained and has the same issues. There is actually a whole section on their github explaining this plotting issue is a know problem with no good solution. That’s why I pulled the data directly into mne and just avoided using that package.
Awesome, got it working! Thanks for the help. It also seems like using the layout was a mistake to begin with. The maps it produces are vastly different (but are unaffected by the info[‘chs’] xy swap), not sure why this is. Prior to the xy swap, the function would not accept the info in the pos argument, which was why I defaulted to the layout originally. Here are the maps I get when I use the default MNE Magness layout.pos in the function’s pos argument.
Here is the code for anyone else trying to do the same thing, or if anyone identifies an egregious error, let me know. But I believe it does what it is supposed to do:
Hey, I just wanted to update this since I have gained a better understanding of what is being changed, and this has major implications on data analysis for any one that wants to use this code on the HCP data. The x/y swap here does not only rotate the sensors 90 degrees. It also swaps the left and right sides. Thus, the displays are flipped. This causes issues with interpolation. If you perform this swap and then interpolate, the interpolated values will be pulled from sensors that are not true neighbors. This will be obvious when you look at the timeseries or PSD for the interpolated sensors. Thus, I recommend using this swap at the very end of your analysis for viewing only. Do NOT run this before any cleaning or analysis.
There is likely a better way to rotate the sensors, but unfortunately the channel locations are in a 12 dimensional format, so without a proper dig file and the HCP’s input, it is not clear how to change these locations without unforeseen consequences, such as swapping the right and left sides, as this code does.
Alright, so the solution was simple, and seems to be a realistic fix, as the interpolation of bad electrodes gives realist results. Here is the working code:
for ch_idx in range(len(epoched.info['chs'])):
current_pos = copy.deepcopy(epoched.info['chs'][ch_idx]['loc'][:2])
epoched.info['chs'][ch_idx]['loc'][0] = current_pos[1]*-1
epoched.info['chs'][ch_idx]['loc'][1] = current_pos[0]
you can then interpolate bad sensors accurately with:
#%% set DIG with Magnes 3600wh sensor loc for
bti_path = op.abspath('/your path to/mne-python-main/mne') + '/io/bti/tests/data/'
raws = {'Magnes 3600wh': read_raw_bti(op.join(bti_path, 'test_pdf_linux'),
op.join(bti_path, 'test_config_linux'),
op.join(bti_path, 'test_hs_linux')),}
dig=mne.channels.DigMontage(dig=raws['Magnes 3600wh'].info['dig'])
epoched.set_montage(dig) #set a digitization file
