Since mask_nans is used twice in inverted form, and never in its original form, perhaps it is worth it to invert at the original isnan location rather than inverting the original value twice?

I completely missed it and it's a good point. You mean, having this instead, correct?

   mask_nans = ~np.isnan(array_surface) #invert original

   index_array = np.arange(array_surface.shape[0])

   interpolated_array = np.interp(index_array,
                                  np.flatnonzero(mask_nans),   #not inverted
                                  array_surface[mask_nans]).   #not inverted

For arrays of floats, the docstring for assert_almost_equal suggests using i.e., assert_allclose these days for better consistency.

Thanks. I remember @orbeckst mentioned that MDAnalysis typically uses assert_almost_equal. Maybe we can have his opinion on what would be best here?

MDAnalysis uses assert_almost_equal in preference to assert_array_almost_equal. assert_allclose might be more recent than when we started using numpy tests so MDAnalysis might just be outdated and should be switching eventually. Follow @tylerjereddy 's advice :-).

And with the others too -- if you specify frame_index, you'll wind up only getting the mean of the last frame.

Do you need np.nanmean? Do you expect np.nan values in your interpolated surfaces? If not, using np.mean will reveal errors in case there ever are np.nan values.

Over which axis?

I'm not quite understanding this. mask_nans is 2D, with shape (n_x, n_y). index_array is 1D, with shape n_x. np.flatnonzero(mask_nans) and array_surface[mask_nans] will be 1D, of length between 0 to n_x * n_y. So I think that means that you are only ever operating on the first row of array_surface?

I guess my question is then, why is array_surface 2D? Why not just pass in one row at a time? In addition, as it is a 2D surface, why use a 1D interpolation function? Why not a 2D one?

Ok, I see that you are applying interpolation_by_array by row, then. In that case I suggest amending the documentation of interpolation_by_array to specify that the input is 1-dimensional with length n_x, instead of n_x, n_y. Although I think doing 2D interpolation would be less arbitrary (choosing which axis is x and which is y is basically a random, right?)

Could you please add some tests with more interesting values? e.g. in the one below, the 145 number could have come from interpolating either along the x-axis or the y-axis. In order to be diagnostic of potential current and future bugs, I'd be interested in seeing one where using 1D interpolation gives different results if it's run across the x-axis, vs. the y-axis.

(np.array([[150., 150,  150.,  150.],
               [150., np.nan, np.nan,  150.],
               [150., 130., 140.,  150.],
               [150., 150., 150.,  150.]]),
     np.array([[150., 150, 150.,  150.],
               [150., 140., 145.,  150.],
               [150., 130., 140.,  150.],
               [150., 150., 150.,  150.]])),

Here as well, tests where the expected values aren't all one number are needed to really check assumptions such as:

• is it interpolating along the axis I think it is
• is it interpolating in a linear way (vs. some random other polynomial)
• is it even interpolating, or just taking the most common number?
• are we sure it's interpolating independently for every frame, or is it interpolating across different frames?

etc.

I am, incidentally, somewhat surprised that coordinates are getting wrapped with np.nan values in them. @richardjgowers should this be happening...?