numpy: Efficiency of approaches for 2D cut of 3D arrays?

Question

Regarding a previous thread: How can an almost arbitrary plane in a 3D dataset be plotted by matplotlib?

I found the thread very helpful, but various solutions are offered and it is not clear to me which one is the most resource efficient. As there are many related problems in physics and natural sciences generally, I would argue the general interest is there (although, a question is whether python is ideal for this heavy-lifting task?).

If a speed comparison of various methods in the quoted thread could be executed by an expert, that would be very helpful! Often, large datasets are involved in this procedure (in my case, 64x64x2000).

Thank you!

Answer 1

Here is my assessment of the solutions offered there:

• Use scipy.interpolate.griddata - this is designed to handle arbitrary scattered points in the input data. For interpolation, it needs to figure out which input points are the nearby the requested points. This is trivial for input data in a 3D array format, but this function doesn't know about that. This will be slow.

• Explicit Python code to calculate the nearest points - this is algorithmitally efficient, but it suffers from being interpreted Python code and it has some rough edges: It uses nearest-neighbor interpolation and it doesn't handle the edges very well.

• Reduced requirements - this algorithm only takes slices aligned with the axes, but it is fast.

• More explicit Python code for nearest-neighbor interpolation - this solution is quite long. I didn't analyze it in detail; it looks like it handles edge cases better but it is still interpreted Python code that doesn't use vectorization.

• Vectorized implementation with optional linear interpolation I just added that implementation. This is vectorized and therefore much, much faster.