How to optimise/vectorize interpolant query on a multidimensional variable?

Dear all,

I have rectangular 2D-space gridded data that I want to interpolate/evaluate on a list of scattered (custom) Gridpoints.

The original (sample) grid is in 2D (longitude-latitude), and the variable to interpolate is organised in several dimensions:

%Precip Dimensions (Lon, lat, date, member, horizon) of lengths [5, 6, 8000, 4, 40]
%Query_points of lengths [820x2] -> 820 query points of 2 dimension each(lon-lat) 

As Mathworks suggests, it is faster to "create interpolant F and evaluate multiple times" than to "compute interpolation separately"

Altought preparing the 2D gridded interpolant is pretty straigthforward, I'm at a loss as to choose the best approach when the variable to interpolate is multidmensionnal.

So far, I can do this:

%preparing 2D-space grids
lon = [-75	-73,5	-72	-70,5	-69]; 
lat = [46,5	48	49,5	51	52,5	54];

[Sg1, Sg2] = ndgrid(lon,lat); %Sample grid
Query_points = [-74.1; 47,
             -74.12; 46.99,
             ...] %etc. 820x2 matrix of lon-lat pairs of points

for horizon = 1:40
    for member = 1:4
        for date = 1:8000
            
            %Extracting a single dimension at a time of the data
            temp_precip = squeeze(precip(:,:,date,member,horizon));
            
            %creating Interpolant object
            F = griddedInterpolant(Sg1,Sg2, temp_precip);
            
            %evaluate on query points
            Q_precip = F(Query_points); % a 820x1 vector of interpolated data evaluated at the requested scattered points
            
            %saving results somewhere
            Final_data(:,date) = Q_precip;
            
        end
        
        %other operations on Final_data goes here
        
    end
end

This approach works, but feels inneficcient (as well as taking a few seconds per the innermost dateloop, making the whole operation several minutes long). Is there a way to vectorize or unroll the date loop for faster execution? (the outmost and middle loop have to stay because of other operations I've ommited for the sake of simplicity)

I've tried creating a 3D interpolant by including the date dimension as a "virtual" 3D-space interpolant, such as

%preparing 3D-space grids
lon = [-75	-73.5	-72	-70.5	-69]; 
lat = [46.5	48	49.5	51	52.5	54];
date = [720000:1:728000]; %serial date numbers

[Sg1, Sg2, Sg3] = ndgrid(lon,lat,date); %Sample grid
Query_points = %some repmat of query lon-lat pair repeated 8000 times (1 for each date), ending in a 6560000*3 array

But this feels unnecessarily bloated, as I'm not evaluating "in-between" dates (the sample and query dates are identical) and the resulting output can be a very very large array.

TL;DR In other words, what is the best approach to evaluate a multidimensional array of data representing values of a variable expressed in 2D space, if not looping on every dimensions of the data?

Thank you, I've searched around but can't wrap my head around this problem,

Simon.

Answer 1

[Simon Matte]

Update:

As of Matlab2021, griddedInterpolant now accepts N-D arrays to be interpolated in a 2-D grid, and operating each queried 2D output by interpolating on each 'page' of the N-D array.