Working with Scientific Data

Data as text

Instruments and simulations commonly save results as CSV (comma-separated values) files: plain text where each line is a row and commas separate the columns, often with a header line describing them. NumPy's np.loadtxt() reads such data straight into a 2D array.

On your own computer you would write np.loadtxt("expansion.csv", delimiter=",", skiprows=1). In these examples the data is written into the program as a string, and io.StringIO() makes that string behave like a file, so the same loadtxt call works.

import ioimport numpy as npcsv_text = """temperature_c,length_mm20,100.000100,100.185200,100.415300,100.648400,100.880500,101.115"""data = np.loadtxt(io.StringIO(csv_text), delimiter=",", skiprows=1)print(data.shape)print(data[:2])

delimiter="," says the columns are separated by commas, and skiprows=1 skips the header line, which contains text rather than numbers.

Selecting columns and converting units

Give each column a meaningful name by slicing it out of the table. Then unit conversions apply to the whole column at once.

import ioimport numpy as npcsv_text = """temperature_c,length_mm20,100.000100,100.185200,100.415300,100.648400,100.880500,101.115"""data = np.loadtxt(io.StringIO(csv_text), delimiter=",", skiprows=1)temperature_k = data[:, 0] + 273.15length_m = data[:, 1] / 1000strain = (length_m - length_m[0]) / length_m[0]print("T (K):", temperature_k)print("Strain (%):", np.round(strain * 100, 3))

Fitting a straight line

Many physical laws are linear over a useful range, and the slope of the line is the quantity you want. For thermal expansion, strain = α × ΔT, so the slope of strain against temperature is the expansion coefficient α. np.polyfit(x, y, 1) finds the best-fit straight line (a polynomial of degree 1) by least squares and returns the slope and intercept.

import ioimport numpy as npcsv_text = """temperature_c,length_mm20,100.000100,100.185200,100.415300,100.648400,100.880500,101.115"""data = np.loadtxt(io.StringIO(csv_text), delimiter=",", skiprows=1)temperature = data[:, 0]length = data[:, 1]strain = (length - length[0]) / length[0]slope, intercept = np.polyfit(temperature, strain, 1)print(f"Expansion coefficient: {slope:.2e} 1/K")print(f"That is {slope * 1e6:.1f} × 10⁻⁶ /K")

Always compare a fitted result with reference values. 23.2 × 10⁻⁶ /K matches aluminium (23.1 × 10⁻⁶ /K) and is far from copper (16.5 × 10⁻⁶ /K) or steel (about 12 × 10⁻⁶ /K), which is a useful check that the sample is what you think it is.

How good is the fit?

np.polyval() evaluates the fitted line at your x-values. The differences between the data and the line are the residuals; small, random residuals mean the straight line describes the data well.

import numpy as nptemperature = np.array([20, 100, 200, 300, 400, 500])strain = np.array([0.0, 0.00185, 0.00415, 0.00648, 0.00880, 0.01115])coefficients = np.polyfit(temperature, strain, 1)fitted = np.polyval(coefficients, temperature)residuals = strain - fittedprint("Largest residual:", f"{np.max(np.abs(residuals)):.1e}")

Cleaning data

Real datasets contain problems: a sensor drops out and records 0, a value is missing, or a reading is physically impossible. NumPy represents missing values as np.nan (not a number). Boolean masks remove invalid rows before analysis, and np.isnan() finds missing values.

import numpy as nphardness = np.array([152, 160, np.nan, 158, 0, 149, 155, 1450])   # HVvalid = ~np.isnan(hardness) & (hardness > 0) & (hardness < 1000)clean = hardness[valid]print(f"Kept {clean.size} of {hardness.size} readings: {clean}")print(f"Mean hardness: {clean.mean():.1f} HV")

~ means not for Boolean arrays, so ~np.isnan(hardness) is True for values that are present. Always report how many points you removed and why: silently discarding data is poor scientific practice.

Finding key values

np.max() gives the largest value, and np.argmax() gives its index, which you can use to look up the matching value in another column. From a tensile test, for example, this gives the ultimate tensile strength and the strain at which it occurred.

import numpy as npstrain = np.array([0.000, 0.002, 0.004, 0.010, 0.050, 0.100, 0.150, 0.180])stress = np.array([0, 140, 270, 300, 330, 345, 340, 310])   # MPapeak_index = np.argmax(stress)print(f"Ultimate tensile strength: {stress[peak_index]} MPa")print(f"Reached at a strain of {strain[peak_index]:.3f}")elastic = strain <= 0.004modulus, _ = np.polyfit(strain[elastic], stress[elastic], 1)print(f"Young's modulus: {modulus / 1000:.1f} GPa")

Assigning to _ is a convention for a value I don't need: here, the intercept of the fit.