Combining two gists into one practical for cathodoluminescence plots, https://gist.github.com/FilipDominec/9101258e8ba979c715a28bc527a43961 and https://gist.github.com/FilipDominec/0625e89c171e402741820dcaf11db77e

cathodoluminescence_contours.py

matplotlib.rc('font', size=12)
ys = np.array(ys)[:, :1024]
ys /= 5782000/.1517*0.02301 # np.max(ys)    # fixed normalization for all samples
x = np.poly1d([4.64708212e-15, -1.65806129e-05, 5.20397778e-01, 3.53568373e+02])(range(1024))



for yy in ys:
    yy[:] = np.convolve(yy, 2**(-np.linspace(-2,2,15)**2),mode='same')

def nGaN(energies):
    ## returns index of refraction in GaN according to [Tisch et al., JAP 89 (2001)]
    eV = [1.503, 1.655, 1.918, 2.300, 2.668, 2.757, 2.872, 3.006, 3.136, 3.229, 3.315, 3.395, 3.422] 
    n  = [2.359, 2.366, 2.383, 2.419, 2.470, 2.486, 2.511, 2.549, 2.596, 2.643, 2.711, 2.818, 2.893] 
    return np.interp(energies, eV, n)
for y in         ys:
    ## initially, x is free space wavelength (in nanometers)
    ## convert x1 to wavenumber in GaN
    x1 = 2*np.pi*nGaN(1239.8/x)/(x*1e-9)   # (1 eV photon has wavelength of 1239.8 nm in vacuum)

    ## interpolate so that the x-axis is in increasing order and equidistant
    x2 = np.linspace(np.min(x1), np.max(x1), len(x1)*4)  ## finer resample to keep all spectral details
    y2 = np.interp(x2, x1[::-1], (y[::-1]))

    ## using FFT, x2 is transformed roughly into "thicknesses of virtual Fabry-Perot resonators"
    xf = np.fft.fftfreq(len(x2), x2[1]-x2[0])
    yf = np.fft.fft(y2)
    
    ## now find and remove the spectral peak
    findmax = np.argmax(np.abs(yf[20:-20])) + 20
    yf[findmax-5:findmax+5] = (yf[findmax-5] + yf[findmax+5])/2
    yf[len(yf)-findmax-5:len(yf)-findmax+5] = (yf[findmax-5] + yf[findmax+5])/2

    # invert FFT of spectra and re-interpolate wavenumber in GaN back to wavelength
    y[:] = np.interp(2*np.pi*nGaN(1239.8/x)/(x*1e-9), x2, (np.abs(np.fft.ifft(yf))))  ## np.fft.fftshift(xf)




x = 1240./x ## crop first row, convert to eV
ys = ys / x**2 ## crop first row and normalize by dE/dW






cmaprange1, cmaprange2 = np.min(ys), np.max(ys) 
labels = [float(ll.split('k')[0]) for ll in labels]
print(labels)

levels = np.linspace(cmaprange1, cmaprange2, 50) 
contourf = ax.contourf(x, labels, ys, levels=levels, extend='both',cmap=cm.jet,alpha=.4)
ax.contour(x, labels, ys, levels=levels, extend='both',cmap=cm.jet,alpha=.4)



## Here we define the desired nonlinear dependence between dual y-axes:
def right_tick_function(p): return 10.46*(p**1.68)
ax2 = ax.twinx()
#ax2.axis['right'].major_ticklabels.set_visible(False)
ax2.set_ylim(np.array(ax.get_ylim()))
ax2.set_ylabel('electron penetration depth (nm)')
## If we wish nice round numbers on the secondary axis ticks...
right_ax_limits = sorted([right_tick_function(lim) for lim in ax.get_ylim()])
yticks2 = matplotlib.ticker.MaxNLocator(nbins=18, steps=[1,2,5]).tick_values(*right_ax_limits)
## ... we must give them correct positions. To do so, we need to numerically invert the tick values:
from scipy.optimize import brentq
def right_tick_function_inv(r, target_value=0): 
    return brentq(lambda r,q: right_tick_function(r) - target_value, ax.get_ylim()[0], ax.get_ylim()[1], 1e6)
valid_ytick2loc, valid_ytick2val = [], []
for ytick2 in yticks2:
    try:
        valid_ytick2loc.append(right_tick_function_inv(0, ytick2))
        valid_ytick2val.append(ytick2)
    except ValueError:      ## (skip tick if the ticker.MaxNLocator gave invalid target_value to brentq optimization)
        pass
## Finally, we set the positions and tick values
ax2.set_yticks(valid_ytick2loc)
from matplotlib.ticker import FixedFormatter
ax2.yaxis.set_major_formatter(FixedFormatter(["%g" % righttick for righttick in valid_ytick2val]))



## The same approach for the x-axes (relation between photon energy and vavelength)
def top_tick_function(p): return sc.h*sc.c/(p*sc.nano)/sc.eV
ax3 = ax.twiny()
#ax2.axis['top'].major_ticklabels.set_visible(False)
ax3.set_xlim(np.array(ax.get_xlim()))
ax3.set_xlabel('emission wavelength (nm)')
## If we wish nice round numbers on the secondary axis ticks...
print(ax.get_xlim())
print(top_tick_function(1.0))
print([top_tick_function(lim) for lim in ax.get_xlim()])
top_ax_limits = sorted([top_tick_function(lim) for lim in ax.get_xlim()])
print(top_ax_limits)
xticks2 = matplotlib.ticker.MaxNLocator(nbins=8, steps=[1,2,5]).tick_values(*top_ax_limits)
## ... we must give them correct positions. To do so, we need to numerically invert the tick values:
from scipy.optimize import brentq
def top_tick_function_inv(r, target_value=0): 
    return brentq(lambda r,q: top_tick_function(r) - target_value, ax.get_xlim()[0], ax.get_xlim()[1], 1e6)
valid_xtick2loc, valid_xtick2val = [], []
for xtick2 in xticks2:
    try:
        valid_xtick2loc.append(top_tick_function_inv(0, xtick2))
        valid_xtick2val.append(xtick2)
    except ValueError:      ## (skip tick if the ticker.MaxNLocator gave invalid target_value to brentq optimization)
        pass
## Finally, we set the positions and tick values
ax3.set_xticks(valid_xtick2loc)
from matplotlib.ticker import FixedFormatter
ax3.xaxis.set_major_formatter(FixedFormatter(["%g" % toptick for toptick in valid_xtick2val]))






ax.set_xlim(np.min(x), np.max(x))
ax.set_ylim(np.min(labels), np.max(labels))

ax.set_xlabel('emission energy (eV)')
ax.set_ylabel('acceleration voltage (kV)')
sample_name = labels_orig[0].split("/")[-2]
ax.set_title('Sample %s cathodoluminescence (Fabry-Pérot filtered)\n\n' % sample_name)

#axR = fig.add_subplot(122)
cax = fig.add_axes([0.8, 0.1, 0.02, 0.5])
fig.colorbar(contourf, extend='both', cax=cax, ax=[ax,ax2,ax3]) ## todo: overlaps the right axis , fraction=0.07


tosave.append(sample_name+'_CL.png')
tosave.append(sample_name+'_CL.pdf')