I have an array of shape(201,201), I would like to plot some cross sections through the data, but I am having trouble accessing the relevant points. For example say I want to plot the cross section given by the line in the figure produced by,
from pylab import *
Z = randn(201,201)
x = linspace(-1,1,201)
X,Y = meshgrid(x,x)
pcolormesh(X,Y,Z)
plot(x,x*.5)
I'd like to plot these at various orientations but they will always pass through the origin if that helps...
Basically, you want to interpolate a 2D grid along a line (or an arbitrary path).
First off, you should decide if you want to interpolate the grid or just do nearest-neighbor sampling. If you'd like to do the latter, you can just use indexing.
If you'd like to interpolate, have a look at scipy.ndimage.map_coordinates. It's a bit hard to wrap your head around at first, but it's perfect for this. (It's much, much more efficient than using an interpolation routine that assumes that the data points are randomly distributed.)
I'll give an example of both. These are adapted from an answer I gave to another question. However, in those examples, everything is plotted in "pixel" (i.e. row, column) coordinates.
In your case, you're working in a different coordinate system than the "pixel" coordinates, so you'll need to convert from "world" (i.e. x, y) coordinates to "pixel" coordinates for the interpolation.
First off, here's an example of using cubic interpolation with map_coordinates:
import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt
# Generate some data...
x, y = np.mgrid[-5:5:0.1, -5:5:0.1]
z = np.sqrt(x**2 + y**2) + np.sin(x**2 + y**2)
# Coordinates of the line we'd like to sample along
line = [(-3, -1), (4, 3)]
# Convert the line to pixel/index coordinates
x_world, y_world = np.array(zip(*line))
col = z.shape[1] * (x_world - x.min()) / x.ptp()
row = z.shape[0] * (y_world - y.min()) / y.ptp()
# Interpolate the line at "num" points...
num = 1000
row, col = [np.linspace(item[0], item[1], num) for item in [row, col]]
# Extract the values along the line, using cubic interpolation
zi = scipy.ndimage.map_coordinates(z, np.vstack((row, col)))
# Plot...
fig, axes = plt.subplots(nrows=2)
axes[0].pcolormesh(x, y, z)
axes[0].plot(x_world, y_world, 'ro-')
axes[0].axis('image')
axes[1].plot(zi)
plt.show()
Alternately, we could use nearest-neighbor interpolation. One way to do this would be to pass order=0 to map_coordinates in the example above. Instead, I'll use indexing just to show another approach. If we just change the line
# Extract the values along the line, using cubic interpolation
zi = scipy.ndimage.map_coordinates(z, np.vstack((row, col)))
To:
# Extract the values along the line, using nearest-neighbor interpolation
zi = z[row.astype(int), col.astype(int)]
We'll get:
