How to reshape and interpolate filtered 2D data

I have filtered PIV data in .dat files. There are four relevant columns, the x-coordinate, (column two), y-coordinate (column 3), and two velocities (columns four and five). I've attached one for reference. There are 1,000 of these files, each representing a moment in time. I ultimately want to get these into a format that can be interpolated and then saved as two 3D matrices of (2D) space and time, one for each velocity. I need to perform a 2D interpolation to fill in the values that were filtered out, so I've been trying interp2, but to no avail. Here's one thing I tried:

 % 
clear; close all;
files = dir('D:\Dan\8_2steady\xyf*.dat'); % load files
pivgridI = importdata('grdUsedI_001.dat',' ',2); 
pivgridJ = importdata('grdUsedJ_001.dat',' ',2);
xgrid = pivgridJ.data(1,:);
ygrid = pivgridI.data(:,1);
n = length(xgrid); m = length(ygrid); % expected size of space array
L = length(files);
for i = 1:L
    Data{i} = importdata(files(i).name);
    x{i} = Data{i}(:,2);
    y{i} = Data{i}(:,3);
    u{i} = Data{i}(:,4);
    v{i} = Data{i}(:,5);  
end
xi = repmat(xgrid,1,m)'; yi = repmat(ygrid',n,1); yi = yi(:);
for i = 1:L
    ui{i} = interp2(x{i},y{i},u{i},xi,yi,'linear');
end

This tells me:

Error using griddedInterpolant

The grid vectors are not strictly monotonic increasing.

Error in interp2>makegriddedinterp (line 228)

F = griddedInterpolant(varargin{:});

And it's right, they're not monotonically increasing. My "xgrid" and "ygrid" represent pixel column and row numbers, respectively. Should I reshape first? And how when they are data filtered out? Something else?

UPDATE: I think I solved this with:

    %
for i = 1:L
    ui{i} = griddata(x{i},y{i},u{i},xgrid,ygrid,'linear');
end

This is pretty slow, though. Anyway to speed it up; I've got a lot of these to do! Also, it looks like griddata does not support extrapolation, which I definitely need. Alternatives?

FINAL UPDATE: ScatteredInterpolant is the way to go. Final code:

    % 
ui = zeros(m,n,L); vi = ui;
[Xq,Yq] = meshgrid(xgrid,ygrid);
for i = 1:L
    F = scatteredInterpolant(x{i},y{i},u{i},'linear','linear');
    G = scatteredInterpolant(x{i},y{i},v{i},'linear','linear');
    ui(:,:,i) = F(Xq,Yq);
    vi(:,:,i) = F(Xq,Yq);
end

And it reshapes it in the process, so there's a bonus!