Describe the bug

Attempt to compute CSD on evoked object called "monopolar_av_ref.fif" (see attached). The evoked object contains grand-averaged EEG waveforms (64 EEG channels (averaged reference), 'biosemi64' montage + 1 EOG channel).

Steps to reproduce

GA_comp_corr_stim = mne.read_evokeds('monopolar_av_ref.fif')
GA_comp_corr_stim_CSD = mne.preprocessing.compute_current_source_density(GA_comp_corr_stim)#should be similar to the attached CSD.fif file

plots obtained:
times = np.arange(.3, 1.,.05)
fig1 = GA_comp_corr_stim.plot_topomap(times = times ,ch_type = 'eeg', cmap = 'viridis', nrows = 2, average = .03, time_unit='s')
fig3 = GA_comp_corr_stim_CSD.plot_topomap(times = times ,ch_type = 'csd', cmap = 'viridis', nrows = 2, average = .03, time_unit='s')

Actual results

plots "topo_CSD.png" and "topo_monopolar.png" in the attached. It is clear that CSD estimates are wrong, since the spatial selectivity appears lower than the spatial selectivity of the monopolar data, and differ from those obtained with a python implemenbtation of Mike Cohen's matlab code ('topo_CSD_Cohen.png')
Maybe it's a problem with the sphere= auto argument that returns:
Fitted sphere radius: 95.0 mm
Origin head coordinates: 0.0 -0.0 40.1 mm
Origin device coordinates: 0.0 -0.0 40.1 mm

Additional information

Platform: Windows-10-10.0.18362-SP0
Python: 3.7.6 (default, Jan 8 2020, 20:23:39) [MSC v.1916 64 bit (AMD64)]
Executable: C:\Users\mservant\Anaconda3\pythonw.exe
CPU: Intel64 Family 6 Model 158 Stepping 10, GenuineIntel: 12 cores
Memory: 15.8 GB

mne: 0.20.5
numpy: 1.18.1 {blas=mkl_rt, lapack=mkl_rt}
scipy: 1.4.1
matplotlib: 3.1.3 {backend=Qt5Agg}

sklearn: 0.22.1
numba: 0.48.0
nibabel: 3.0.1

cupy: Not found
pandas: 1.0.1
dipy: 1.1.1
mayavi: 4.7.1 {qt_api=pyqt5, PyQt5=5.9.2}
pyvista: Not found
vtk: 8.1.2

mne.zip