How to (re)scale the x-axis to fit certain points in the graph?

Question

I want to rescale my (qualitative) x-axis so the two peaks (visible in the graph) correlate to their actual values (i. e. 500 keV and 1274 MeV). How can I do this?

import numpy as np
import matplotlib.pyplot as plt

def read_from_file(filename):
    return np.loadtxt(filename)

data = list(read_from_file("calibration.txt"))

print(data.index(max(data[:2000])))#x value 500kev
print(data.index(max(data[2000:])))#x value 1274

fig = plt.figure()
ax = fig.add_subplot(111)
x = range(len(data))
plt.plot(x, data)
plt.xlim(0, 5000)
plt.ylim(0, 7000)
plt.title("$^{22}$Na Spectrum")
plt.xlabel("Energy")
plt.ylabel("Amount of Photons")
plt.grid()
ax.annotate("500 keV", xy = (1450, 6541), xytext = (1600, 6500))
ax.annotate("1274 MeV", xy = (3500, 950), xytext = (3700, 1100))
plt.show()

Answer 1

Funny you should ask this question. I am currently trying to push an example into MatPlotLib that shows exactly how to do this. You can view the recipe here: https://github.com/madphysicist/matplotlib/blob/7b05223c85741120019b81e1248c20f9bc090c61/examples/ticks_and_spines/tick_transform_formatter.py

You do not need the entire code in the example (or the tick formatter that uses it) but the mapping function will help you create the scaled x-array (also, use np.argmax instead of index(max(...)):

ind500 = np.argmaxmax(data[:2000]))
ind1274 = np.argmax(data[2000:])) + 2000
x_scaled = (x - ind500) * (1274 - 500) / (ind1274 - ind500) + 500

You can use x_scaled to plot as usual:

plt.plot(x_scaled, data)
...

Combining it all together (and making a couple of tweaks to use OO API instead of pyplot):

import numpy as np
from matplotlib import pyplot as plt

data = np.loadtxt("calibration.txt") # Don't convert this back to a list

ind500 = np.argmaxmax(data[:2000]))
ind1274 = np.argmax(data[2000:])) + 2000
x = (np.arange(len(data)) - ind500) * (1274 - 500) / (ind1274 - ind500) + 500

fig, ax = plt.subplots()
ax.plot(x, data)
plt.title("$^{22}$Na Spectrum")
plt.xlabel("Energy")
plt.ylabel("Photons Counts")
plt.grid()
ax.annotate("500 keV", xy = (500, data[ind500]), xytext = (550, data[ind500] + 100))
ax.annotate("1274 keV", xy = (1274, data[ind1274]), xytext = (1324, data[ind1274] + 100))
plt.show()

The example I linked to would allow you to display the x-axis in entirely different units without actually modifying your x-array.

Answer 2

Using numpy, you can find the index of the two spikes (i.e. no need to convert the data to a list) using argmax.

Then, you can scale the x values using:

xnew = val1 + (x - max1) / (max2 - max1) * (val2 - val1)

where val1 and val2 are the values of your peaks, and max1 and max2 are the indices of those peaks.

Here's a bit of code that should work:

import numpy as np
import matplotlib.pyplot as plt

# Fake some data approximately in your range. You can ignore this bit!
# Random numbers for noise
data = 1000. + np.random.rand(5000) * 100.
x = np.arange(len(data))
# Add the first spike
mu1, sd1 = 1450., 300.
pdf1 = (1./(sd1*2.*np.pi) * np.exp(-(x - mu1)**2 / sd1**2)) * 1e7
data += pdf1
# Add the second spike
mu2, sd2 = 3500., 200.
pdf2 = (1./(sd2*2.*np.pi) * np.exp(-(x - mu2)**2 / sd2**2)) * 1e6
data += pdf2
# End of fake data generation

# Find the index of the first maximum (using your '2000' cutoff)
cutoff = 2000
max1 = float(np.argmax(data[:cutoff]))
# Find the index of the second cutoff
max2 = float(np.argmax(data[cutoff:]) + cutoff)

# The actual values of the two spikes
val1, val2 = 500., 1274

# Scale the xvalues
xnew = val1 + (x - max1) / (max2 - max1) * (val2 - val1)

# Plot
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xnew, data)
ax.set_ylim(0, 7000)
ax.set_title("$^{22}$Na Spectrum")
ax.set_xlabel("Energy")
ax.set_ylabel("Number of Photons")
ax.grid()

# Add some lines at the actual spikes to check scaling worked
ax.axvline(val1)
ax.axvline(val2)

plt.show()